总的思路就是把Fn 转化成用Fj 表示的线性递推的形式。
然后用矩阵快速幂加速。
具体分析参考这个博主(分析的很棒!)
不过所有能用矩阵快速幂的线性递推。。
都能用杜教BM开挂秒。。。
https://blog.csdn.net/The___Flash/article/details/106066577
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define ls (o<<1)
#define rs (o<<1|1)
#define pb push_back
const double PI= acos(-1.0);
const int M = 1e5+7;
/*
int head[M],cnt;
void init(){cnt=0,memset(head,0,sizeof(head));}
struct EDGE{int to,nxt,val;}ee[M*2];
void add(int x,int y){ee[++cnt].nxt=head[x],ee[cnt].to=y,head[x]=cnt;}
*/
const int mod = 998244353;
int sz;
struct mat{
ll a[110][110];
};
mat operator *(mat &a, mat &b){
mat ans={0};
for(int i=0;i<sz;i++)for(int j=0;j<sz;j++){
for(int k=0;k<sz;k++) {
ans.a[i][j] += a.a[i][k] * b.a[k][j] % mod;
}
ans.a[i][j]%=mod;
}
return ans;
}
ll C[210][210],f[210],pown[210],g[210];
ll qpow(ll a,ll b)
{
ll ans=1;
while(b)
{
if(b&1)ans=ans*a%mod;
a=a*a%mod;
b/=2;
}
return ans;
}
int main()
{
ll n,k;
cin>>n>>k;
C[0][0]=pown[0]=1;
ll mn=n%mod;
for(int i=1;i<=200;i++)
{
pown[i]=pown[i-1]*mn%mod;
C[i][0]=1;
for(int j=1;j<=200;j++)
C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;
}
mat tp,ans;
tp.a[0][0]=0;tp.a[0][1]=1;
tp.a[1][1]=1;tp.a[1][0]=1;
tp.a[1][2]=1;tp.a[2][2]=1;
for(int i=1;i<=k;i++)
for(int j=0;j<=i;j++)
tp.a[j+2][i+2]=C[i][j];
//G(j,1)={1,1,1,1,1}
ll b=n-1;//需要乘tp的次数
sz=k+3;
for(int i=0;i<sz;i++)ans.a[i][i]=1;
while(b)
{
if(b&1)ans=ans*tp;
tp=tp*tp;
b/=2;
}
for(int i=0;i<=k;i++)
for(int j=0;j<=2;j++)
g[i]=(g[i]+ans.a[j][i+2])%mod;
// for(int i=0;i<=k;i++)cout<<i<<" "<<g[i]<<endl;
f[0]=g[0];
for(int i=1;i<=k;i++)
{
ll nw=g[i];
for(int j=0;j<i;j++)
{
int z=1;
if(j&1)z=-1;
nw=(nw-C[i][j]*z%mod*pown[i-j]%mod*f[j]%mod+mod)%mod;
}
f[i]=nw;
if(i&1)f[i]*=-1;
f[i]=(f[i]+mod)%mod;
// cout<<"- "<<i<<" "<<f[i]<<endl;
}
cout<<f[k]<<endl;
return 0;
}