[Luogu5320][BJOI2019]堪破神机(DP+斯特林数)

https://www.cnblogs.com/cjyyb/p/10747543.html

特征方程+斯特林反演化简式子，要注意在模998244353意义下5没有二次剩余，所以每个数都要用$a+b\sqrt{5}$的形式表示，运算类似复数。

斯特林反演的几个用法：

1.下降幂转幂：连续求和时可以通过等比数列求和公式加速。

2.幂转下降幂：类似自然数幂和地用有限微积分加速，或帮助设计DP状态。

 #include<cstdio>

 #include<algorithm>

 #define rep(i,l,r) for (int i=(l); i<=(r); i++)

 typedef long long ll;

 using namespace std;

 const int N=,mod=,i2=,i5=,i6=;

 ll V,L,R;

 int T,m,K,S[N][N],C[N][N];

 struct num{ int a,b; }z;

 num operator +(const num &a,const num &b){ return (num){(a.a+b.a)%mod,(a.b+b.b)%mod}; }

 num operator -(const num &a,const num &b){ return (num){(a.a-b.a+mod)%mod,(a.b-b.b+mod)%mod}; }

 num operator *(const num &a,const num &b){ return (num){int((1ll*a.a*b.a+V*a.b*b.b)%mod),int((1ll*a.b*b.a+1ll*a.a*b.b)%mod)}; }

 num operator *(const num &a,int b){ return (num){int(1ll*a.a*b%mod),int(1ll*a.b*b%mod)}; }

 int ksm(int a,ll b){

     int res=;

     for (; b; a=1ll*a*a%mod,b>>=)

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

     return res;

 }

 num ksm(num a,ll b){

     num res=(num){,};

     for (; b; a=a*a,b>>=)

         if (b & ) res=res*a;

     return res;

 }

 num Inv(num a){ return (num){a.a,(mod-a.b)%mod}*ksm((1ll*a.a*a.a-V*a.b*a.b%mod+mod)%mod,mod-); }

 int main(){

     freopen("calc.in","r",stdin);

     freopen("calc.out","w",stdout);

     scanf("%d%d",&T,&m);

     if (m==) V=; else V=;

     while (T--){

         scanf("%lld%lld%d",&L,&R,&K); S[][]=;

         rep(i,,K) rep(j,,i) S[i][j]=(S[i-][j-]-1ll*S[i-][j]*(i-)%mod+mod)%mod;

         rep(i,,K) C[i][]=;

         rep(i,,K) rep(j,,i) C[i][j]=(C[i-][j]+C[i-][j-])%mod;

         num a,b,A,B; ll len=R-L+;

         if (m==) L++,R++,a=(num){i2,i2},b=(num){i2,mod-i2},A=(num){,i5},B=(num){,mod-i5};

         else R>>=,L=(L+)>>,a=(num){,},b=(num){,mod-},A=(num){i2,i6},B=(num){i2,mod-i6};

         int ans=; ll t=R-L;

         rep(i,,K){

             num res=z;

             rep(j,,i){

                 num s=ksm(A,j)*ksm(B,i-j)*C[i][j];

                 num p=ksm(ksm(a,L),j)*ksm(ksm(b,L),(i-j));

                 num q=ksm(a,j)*ksm(b,i-j);

                 num w=p*((ksm(q,t+)-(num){,})*Inv(q-(num){,}));

                 if (q.a== && q.b==) w=p*((R-L+)%mod);

                 res=res+s*w;

             }

             ans=(ans+1ll*res.a*S[K][i])%mod;

         }

         int Ans=1ll*ans*ksm(len%mod,mod-)%mod;

         rep(i,,K) Ans=1ll*Ans*ksm(i,mod-)%mod;

         printf("%d\n",Ans);

     }

     return ;

 }