继续浏览精彩内容
慕课网APP
程序员的梦工厂
打开
继续
感谢您的支持,我会继续努力的
赞赏金额会直接到老师账户
将二维码发送给自己后长按识别
微信支付
支付宝支付

BM

holdtom
关注TA
已关注
手记 1703
粉丝 240
获赞 991


#include<bits/stdc++.h>

using namespace std;

typedef long long ll;

typedef unsigned long long ull;

#define rep(i,a,n) for(ll i=a;i<n;i++)

#define per(i,a,n) for(ll i=n-1;i>=a;i--)

#define pb push_back

#define mp make_pair

#define all(x) (x).begin(),(x).end()

#define fi first

#define se second

#define SZ(x) (ll)(x).size()

typedef vector<ll>VI;

typedef pair<ll,ll>PII;

const int mod=1e9+7;

ll powmod(ll a,ll b){

 ll res=1;

 a%=mod;

 assert(b>=0);

 for(;b;b>>=1){

  if(b&1)res=res*a%mod;

  a=a*a%mod;

 }return res;

}

ll _,n;

namespace gauss{

 const ll N=10010;

 ll res[N],base[N],_c[N],_md[N];

 vector<ll>Md;

 void mul(ll *a,ll *b,ll k){

  rep(i,0,k+k)_c[i]=0;

  rep(i,0,k)if(a[i])rep(j,0,k)

  _c[i+j]=(_c[i+j]+a[i]*b[j])%mod;

  for(ll i=k+k-1;i>=k;i--)

   if(_c[i])

    rep(j,0,SZ(Md))

     _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;

  rep(i,0,k)a[i]=_c[i];

 }

 ll solve(ll n,VI a,VI b){

  ll ans=0,pnt=0;

  ll k=SZ(a);

  assert(SZ(a)==SZ(b));

  rep(i,0,k)_md[k-1-i]=-a[i];_md[k]=1;

  Md.clear();

  rep(i,0,k)if(_md[i]!=0)Md.pb(i);

  rep(i,0,k)res[i]=base[i]=0;

  res[0]=1;

  while((1ll<<pnt)<=n)pnt++;

  for(ll p=pnt;p>=0;p--){

   mul(res,res,k);

   if((n>>p)&1){

    for(ll i=k-1;i>=0;i--)res[i+1]=res[i];res[0]=0;

     rep(j,0,SZ(Md))res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod;

   }

  }

  rep(i,0,k)ans=(ans+res[i]*b[i])%mod;

  if(ans<0)ans+=mod;

  return ans;

 }

 VI BM(VI s){

  VI C(1,1),B(1,1);

  ll L=0,m=1,b=1;

  rep(n,0,SZ(s)){

   ll d=0;

   rep(i,0,L+1)d=(d+1ll*C[i]*s[n-i])%mod;

   if(d==0)++m;

   else if(2*L<=n){

    VI T=C;

    ll c=mod-d*powmod(b,mod-2)%mod;

    while(SZ(C)<SZ(B)+m)C.pb(0);

    rep(i,0,SZ(B))C[i+m]=(C[i+m]+c*B[i])%mod;

    L=n+1-L;B=T;b=d;m=1;

   }

   else {

    ll c=mod-d*powmod(b,mod-2)%mod;

    while(SZ(C)<SZ(B)+m)C.pb(0);

    rep(i,0,SZ(B))C[i+m]=(C[i+m]+c*B[i])%mod;

    ++m;

   }

  }return C;

 }

 ll gao(VI a,ll n){

  VI c=BM(a);

  c.erase(c.begin());

  rep(i,0,SZ(c))c[i]=(mod-c[i])%mod;

  return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));

 }

};

VI g;

ll f[1100];

ll sum[1100];

int main() 

{

 ll n,k;

 f[1]=1;

 f[2]=2;

 cin>>n>>k;

 for(int i=3;i<=200;i++)f[i]=(f[i-1]+f[i-2])%mod;

 g.clear();

 sum[0]=0;

 for(int i=1;i<=200;i++){

  sum[i]=(f[i]*powmod(i,k)+sum[i-1])%mod;

  g.pb(sum[i]);

 }

 printf("%lld\n",gauss::gao(g,n-1));

}

©著作权归作者所有:来自51CTO博客作者qinXpeng的原创作品,如需转载,请注明出处,否则将追究法律责任


打开App,阅读手记
0人推荐
发表评论
随时随地看视频慕课网APP