LOJ 3090 「BJOI2019」勘破神机——斯特林数+递推式求通项+扩域
题目：https://loj.ac/problem/3090
题解：https://www.luogu.org/blog/rqy/solution-p5320
1.用斯特林数把下降幂化为普通的幂次求和
2.找出通项公式，使得幂次变成二项式，进而将 [ l , r ] 的部分变成等比数列求和
3.模 998244353 下没有 \( \sqrt{5} \) ，所以“扩域”，就是把数表示成 \( a+b*\sqrt{5} \) ；\( \sqrt{3} \) 也同理
注意扩域之后，不满足费马小定理，所以快速幂的指数不能对 ( mod-1 ) 取模！！！
还是不太知道怎么求的通项。为什么是 \( f[n]=A*x_{1}^{n}+B*x_{2}^{n} \) 的形式呢？如果不是二阶怎么推？
UPD:
设特征根是 x1,x2,...,xk，因为 x^n 是通解，又有线性性（？），所以通项可以写成 \( f(i)=A*x_1^i + B*x_2^i + ... \)
但是有重根的话就不是这样。 k 重根的系数是次数界为 k 的多项式。这里的次数指的是 i 的几次幂。
（k重根是针对根而言的，比如一个六次方程，x1=x2=2 , x3=5 , x4=x5=x6=1，那么 2 是2重根，5是单根，1是3重根）
比如，\( f_i = 2*f_{i-1} - f_{i-2} \)，\( f_0 = 1 , f_1 = 2 \)
解出特征根是 x1=x2=1 ，那么可以设通项公式为 \( f(i)=(A*i+B)x^i \) ，解得 A=B=1 。
又如，\( f_i = 4*f_{i-1} - f_{i-2} \) ， \( f_0 = 1 , f_1 = 7 \)
解出特征根是 x1=x2=2 ，设通项是 \( f(i)=(A*i+B)x^i \) ，解得 \( A=\frac{5}{2} , B=1 \)
^ ^
#include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; ll rdn() { ll ret=0;bool fx=1;char ch=getchar(); while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return fx?ret:-ret; } const int N=505,mod=998244353; int upt(int x){while(x>=mod)x-=mod;while(x<0)x+=mod;return x;} int pw(int x,ll k) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} int k,s1[N][N],c[N][N],bs,tlen2,ans;ll tl,tlen; struct Node{ int x,y; Node(int x=0,int y=0):x(x),y(y) {} Node operator+ (const Node &b)const { return Node(upt(x+b.x),upt(y+b.y));} Node operator- (const Node &b)const { return Node(upt(x-b.x),upt(y-b.y));} Node operator* (const Node &b)const { return Node(((ll)x*b.x+(ll)bs*y%mod*b.y)%mod,((ll)y*b.x+(ll)x*b.y)%mod);} }A[N],B[N],x1[N],x2[N],one; Node pw(Node x,ll k) { Node ret=Node(1,0); while(k){if(k&1)ret=ret*x;x=x*x;k>>=1;}return ret;} Node Inv(Node u) { int tp=upt(((ll)u.x*u.x-(ll)bs*u.y%mod*u.y)%mod); tp=pw(tp,mod-2); return Node((ll)u.x*tp%mod,upt(-(ll)u.y*tp%mod)); } void init(int lx) { s1[0][0]=1; for(int i=1;i<=k;i++) for(int j=1;j<=i;j++) s1[i][j]=(s1[i-1][j-1]+(ll)s1[i-1][j]*(i-1))%mod; for(int i=0;i<=k;i++)c[i][0]=1; for(int i=1;i<=k;i++) for(int j=1;j<=i;j++) c[i][j]=upt(c[i-1][j-1]+c[i-1][j]); one=Node(1,0); if(lx==2) { int tp=pw(5,mod-2); A[1]=Node(0,tp); B[1]=Node(0,upt(-tp)); tp=pw(2,mod-2); x1[1]=Node(tp,tp); x2[1]=Node(tp,upt(-tp)); } else { int tp=pw(6,mod-2); A[1]=Node((ll)3*tp%mod,tp); B[1]=Node((ll)3*tp%mod,upt(-tp)); x1[1]=Node(2,1); x2[1]=Node(2,upt(-1)); } A[0]=one; for(int i=2;i<=k;i++)A[i]=A[i-1]*A[1]; B[0]=one; for(int i=2;i<=k;i++)B[i]=B[i-1]*B[1]; x1[0]=one; for(int i=2;i<=k;i++)x1[i]=x1[i-1]*x1[1]; x2[0]=one; for(int i=2;i<=k;i++)x2[i]=x2[i-1]*x2[1]; } Node cal(Node x) { if(x.x==1&&x.y==0)return Node(tlen2,0);//tlen2 not tlen!!! Node d=Inv(one-x); d=d*pw(x,tl)*(one-pw(x,tlen)); return d; } void solve2() { ll l=rdn(),r=rdn(); k=rdn(); bs=5; init(2); int iv=pw((r-l+1)%mod,mod-2); l++; r++; tl=l; tlen=(r-l+1); tlen2=(r-l+1)%mod; for(int i=0,fx=((k&1)?upt(-1):1);i<=k;i++,fx=upt(-fx)) { int tp=0; for(int j=0;j<=i;j++) { Node tmp=x1[j]*x2[i-j]; Node d=cal(tmp)*A[j]*B[i-j]; tp=(tp+(ll)c[i][j]*d.x)%mod; } ans=(ans+(ll)s1[k][i]*fx%mod*tp)%mod; } ans=(ll)ans*iv%mod; int ml=1; for(int i=2;i<=k;i++)ml=(ll)ml*i%mod; ans=(ll)ans*pw(ml,mod-2)%mod; } void solve3() { ll l=rdn(),r=rdn(); k=rdn(); bs=3; init(3); int iv=pw((r-l+1)%mod,mod-2); l=(l+1)>>1; r=r>>1; tl=l; tlen=(r-l+1); tlen2=(r-l+1)%mod; for(int i=0,fx=((k&1)?upt(-1):1);i<=k;i++,fx=upt(-fx)) { int tp=0; for(int j=0;j<=i;j++) { Node tmp=x1[j]*x2[i-j]; Node d=cal(tmp)*A[j]*B[i-j]; tp=(tp+(ll)c[i][j]*d.x)%mod; } ans=(ans+(ll)s1[k][i]*fx%mod*tp)%mod; } ans=(ll)ans*iv%mod; int ml=1; for(int i=2;i<=k;i++)ml=(ll)ml*i%mod; ans=(ll)ans*pw(ml,mod-2)%mod; } int main() { int op=rdn(); op=rdn(); if(op==2)solve2(); else solve3(); printf("%d\n",ans); return 0; }
