最主要的两个式子：
套路1：
$$egin{array}{l}x^k=sum_{i=0}^kegin{pmatrix}x\iend{pmatrix}egin{Bmatrix}k\iend{Bmatrix}i!\end{array}$$
左边的式子可以看成将k个球放到x个有标号的盒子中，总共的方案数。
右边的式子可以看成先选出i个非空的盒子，将k个球放到这i个盒子中且不空的方案数。但由于第二类斯特林数是无序的，还要乘上i阶乘。
套路2：
$$egin{array}{l}egin{Bmatrix}n\mend{Bmatrix}=frac1{m!};sum_{i=0}^m;(-1)^i;;egin{pmatrix}m\iend{pmatrix}(m-i)^n\end{array}$$
这是第二类斯特林数的容斥形式。
右边相当于枚举有多少个空的盒子，然后剩下的盒子随便放。由于这是有标号的结果，最后要除以m阶乘。
有了套路2，就能NTT在O(nlogn)时间内算出$$egin{array}{l}\egin{Bmatrix}n\0end{Bmatrix}simegin{Bmatrix}n\mend{Bmatrix}end{array}$$>。
1.求第二类斯特林数S(n,0)~S(n,m)，n,m<=5E5，模998244353。
朴素NTT即可。
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long int ll; 4 const ll mod=998244353; 5 const ll G=3; 6 const ll Gi=332748118; 7 const ll maxn=1E6+5; 8 ll a[maxn],b[maxn],fac[maxn],n,m,len,limit,sum[maxn],ans; 9 int r[maxn]; 10 int re(int x) 11 { 12 int ans=0; 13 for(int i=0;i<len;++i)ans=ans*2+((x&(1<<i))>0); 14 return ans; 15 } 16 ll qpow(ll x,ll y) 17 { 18 ll ans=1,base=x; 19 while(y) 20 { 21 if(y&1)ans=ans*base%mod; 22 base=base*base%mod; 23 y>>=1; 24 } 25 return ans; 26 } 27 void init() 28 { 29 fac[0]=1; 30 for(int i=1;i<=n;++i)fac[i]=fac[i-1]*i%mod; 31 for(int i=0;i<=n;++i)b[i]=qpow(fac[i],mod-2)*qpow(i,n)%mod; 32 for(int i=0;i<=n;++i)a[i]=(qpow(-1,i)*qpow(fac[i],mod-2)%mod+mod)%mod; 33 } 34 void NTT(ll*A,int g) 35 { 36 for(int i=0;i<limit;++i) 37 if(r[i]<i)swap(A[i],A[r[i]]); 38 for(int i=2;i<=limit;i*=2) 39 { 40 ll w; 41 if(g==1)w=qpow(G,(mod-1)/i); 42 else w=qpow(Gi,(mod-1)/i); 43 for(int j=0;j<limit/i;++j) 44 { 45 ll d=1; 46 for(int k=0;k<i/2;++k) 47 { 48 ll a=A[i*j+k],b=A[i*j+i/2+k]*d%mod; 49 A[i*j+k]=(a+b)%mod; 50 A[i*j+i/2+k]=(a-b+mod)%mod; 51 d=d*w%mod; 52 } 53 } 54 } 55 } 56 int main() 57 { 58 ios::sync_with_stdio(false); 59 cin>>n; 60 init(); 61 len=0,limit=1; 62 while(limit<=n*2) 63 { 64 ++len; 65 limit<<=1; 66 } 67 for(int i=n+1;i<=limit;++i)a[i]=0; 68 for(int i=n+1;i<=limit;++i)b[i]=0; 69 for(int i=0;i<limit;++i)r[i]=re(i); 70 NTT(a,1); 71 NTT(b,1); 72 for(int i=0;i<limit;++i)a[i]=a[i]*b[i]%mod; 73 NTT(a,-1); 74 ll g=qpow(limit,mod-2); 75 for(int i=0;i<limit;++i)a[i]=a[i]*g%mod; 76 for(int i=0;i<=n;++i)cout<<a[i]<<" ";cout<<endl; 77 return 0; 78 }
2.n<=1E5,求$$sum_{i=0}^n;sum_{j=0}^i;egin{Bmatrix}i\jend{Bmatrix} imes2^j imes(j!);$$
改变求和顺序，将斯特林数展开即可。
$$egin{array}{l}sum_{i=0}^n;sum_{j=0}^i;egin{Bmatrix}i\jend{Bmatrix} imes2^j imes(j!);\=sum_{j=0}^n;2^j imes(j!)sum_{i=j}^n;egin{Bmatrix}i\jend{Bmatrix}\=sum_{j=0}^n;2^j imes(j!)sum_{i=0}^n;egin{Bmatrix}i\jend{Bmatrix}end{array}$$
令$$egin{array}{l}\f(m)=sum_{i=0}^n;egin{Bmatrix}i\mend{Bmatrix}\f(m)=sum_{i=0}^n;frac{displaystyle1}{displaystyle m!};sum_{j=0}^n;(-1)^jegin{pmatrix}m\jend{pmatrix}(m-j)^i\=sum_{i=0}^nsum_{j=0}^n;frac{displaystyle1}{displaystyle m!}frac{displaystyle m!}{displaystyle j!(m-j)!}(m-j)^i(-1)^j\=sum_{i=0}^nsum_{j=0}^n;frac{(m-j)^i}{(m-j)!} imesfrac{(-1)^j}{j!}\=sum_{j=0}^n;frac{displaystyle(-1)^j}{displaystyle j!}sum_{i=0}^n;frac{displaystyle(m-j)^i}{displaystyle(m-j)!}\end{array}$$
右边只是一个等比数列求和。求出所有的值后，显然可以卷积。
最后答案为$$sum_{i=0}^n;2^i imes(i!) imes f(i)$$
注意等比数列求和时若x=1，需特判。
1 // luogu-judger-enable-o2 2 #include<bits/stdc++.h> 3 using namespace std; 4 typedef long long int ll; 5 const ll mod=998244353; 6 const ll G=3; 7 const ll Gi=332748118; 8 const ll maxn=1E6+5; 9 ll a[maxn],b[maxn],fac[maxn],n,m,len,limit,sum[maxn],ans; 10 int r[maxn]; 11 int re(int x) 12 { 13 int ans=0; 14 for(int i=0;i<len;++i)ans=ans*2+((x&(1<<i))>0); 15 return ans; 16 } 17 ll qpow(ll x,ll y) 18 { 19 ll ans=1,base=x; 20 while(y) 21 { 22 if(y&1)ans=ans*base%mod; 23 base=base*base%mod; 24 y>>=1; 25 } 26 return ans; 27 } 28 ll get(ll x){return (qpow(x,n+1)-1)*qpow(x-1,mod-2)%mod;} 29 void init() 30 { 31 fac[0]=1; 32 for(int i=1;i<=n;++i)fac[i]=fac[i-1]*i%mod; 33 for(int i=0;i<=n;++i)a[i]=(qpow(-1,i)*qpow(fac[i],mod-2)%mod+mod)%mod; 34 for(int i=0;i<=n;++i)b[i]=get(i)*qpow(fac[i],mod-2)%mod; 35 } 36 void NTT(ll*A,int g) 37 { 38 for(int i=0;i<limit;++i) 39 if(r[i]<i)swap(A[i],A[r[i]]); 40 for(int i=2;i<=limit;i*=2) 41 { 42 ll w; 43 if(g==1)w=qpow(G,(mod-1)/i); 44 else w=qpow(Gi,(mod-1)/i); 45 for(int j=0;j<limit/i;++j) 46 { 47 ll d=1; 48 for(int k=0;k<i/2;++k) 49 { 50 ll a=A[i*j+k],b=A[i*j+i/2+k]*d%mod; 51 A[i*j+k]=(a+b)%mod; 52 A[i*j+i/2+k]=(a-b+mod)%mod; 53 d=d*w%mod; 54 } 55 } 56 } 57 } 58 int main() 59 { 60 ios::sync_with_stdio(false); 61 cin>>n; 62 init(); 63 len=0,limit=1; 64 while(limit<=n*2) 65 { 66 ++len; 67 limit<<=1; 68 } 69 for(int i=0;i<limit;++i)r[i]=re(i); 70 b[1]=n+1; 71 NTT(a,1); 72 NTT(b,1); 73 for(int i=0;i<limit;++i)a[i]=a[i]*b[i]%mod; 74 NTT(a,-1); 75 ll g=qpow(limit,mod-2); 76 for(int i=0;i<limit;++i)a[i]=a[i]*g%mod; 77 for(int i=0;i<=n;++i)ans=(ans+qpow(2,i)*fac[i]%mod*a[i]%mod)%mod; 78 cout<<ans<<endl; 79 return 0; 80 }
https://www.luogu.org/problemnew/show/P4091
3.n个点，求（（所有有编号无向图（每个点度数的k次方的和））的和）。
考虑每个点的贡献，答案为：$$n;sum_{i=0}^{n-1};i^k imesegin{pmatrix}n-1\iend{pmatrix} imes2^{egin{pmatrix}n-1\2end{pmatrix}}$$
即枚举每个点的度数，n-1个点中要有i个点连向他，剩下的点随便连边。
方便起见，只需要求
$$egin{array}{l}sum_{i=0}^n;i^k imesegin{pmatrix}n\iend{pmatrix};\=sum_{i=0}^n;egin{pmatrix}n\iend{pmatrix};sum_{j=0}^k;egin{Bmatrix}k\jend{Bmatrix}egin{pmatrix}i\jend{pmatrix}j!\=sum_{j=0}^kegin{Bmatrix}k\jend{Bmatrix}j!;;sum_{i=0}^n;egin{pmatrix}n\iend{pmatrix}egin{pmatrix}i\jend{pmatrix}\=sum_{j=0}^kegin{Bmatrix}k\jend{Bmatrix}j!;;sum_{i=0}^n;egin{pmatrix}n\jend{pmatrix}egin{pmatrix}n-j\i-jend{pmatrix}\=sum_{j=0}^kegin{Bmatrix}k\jend{Bmatrix}j!egin{pmatrix}n\jend{pmatrix};;sum_{i=0}^n;egin{pmatrix}n-j\i-jend{pmatrix};\=sum_{j=0}^kegin{Bmatrix}k\jend{Bmatrix}j!egin{pmatrix}n\jend{pmatrix};;2^{n-j}end{array}$$
最后几步看不懂翻翻具体数学。P143
1 #pragma GCC optimize(2) 2 #include<bits/stdc++.h> 3 using namespace std; 4 typedef long long int ll; 5 const ll mod=998244353; 6 const ll G=3; 7 const ll Gi=332748118; 8 const ll maxn=2E6+5; 9 ll min(ll x,ll y){return x<y?x:y;} 10 ll max(ll x,ll y){return x>y?x:y;} 11 int r[maxn],len,limit; 12 ll fac[maxn],a[maxn],b[maxn],k,n,ans,inv[maxn]; 13 int re(int x) 14 { 15 int ans=0; 16 for(int i=0;i<len;++i)ans=ans*2+((x&(1<<i))>0); 17 return ans; 18 } 19 ll qpow(ll x,ll y) 20 { 21 ll ans=1,base=x; 22 while(y) 23 { 24 if(y&1)ans=ans*base%mod; 25 base=base*base%mod; 26 y>>=1; 27 } 28 return ans; 29 } 30 void init() 31 { 32 fac[0]=1; 33 for(int i=1;i<=k;++i)fac[i]=fac[i-1]*i%mod; 34 for(int i=0;i<=k;++i)inv[i]=qpow(fac[i],mod-2); 35 for(int i=0;i<=k;++i)a[i]=(qpow(-1,i)*inv[i]%mod+mod)%mod; 36 for(int i=0;i<=k;++i)b[i]=qpow(i,k)*inv[i]%mod; 37 limit=1; 38 while(limit<=2*k) 39 { 40 ++len; 41 limit<<=1; 42 } 43 for(int i=0;i<limit;++i)r[i]=re(i); 44 } 45 void NTT(ll*A,int g) 46 { 47 for(int i=0;i<limit;++i) 48 if(i<r[i])swap(A[i],A[r[i]]); 49 for(int i=2;i<=limit;i*=2) 50 { 51 ll w; 52 if(g==1)w=qpow(G,(mod-1)/i); 53 else w=qpow(Gi,(mod-1)/i); 54 for(int j=0;j<limit/i;++j) 55 { 56 ll d=1; 57 for(int k=0;k<i/2;++k) 58 { 59 ll a=A[i*j+k],b=d*A[i*j+i/2+k]%mod; 60 A[i*j+k]=(a+b)%mod; 61 A[i*j+i/2+k]=(a-b+mod)%mod; 62 d=d*w%mod; 63 } 64 } 65 } 66 } 67 int main() 68 { 69 ios::sync_with_stdio(false); 70 cin>>n>>k; 71 --n; 72 init(); 73 NTT(a,1); 74 NTT(b,1); 75 for(int i=0;i<=limit;++i)a[i]=a[i]*b[i]%mod; 76 NTT(a,-1); 77 ll g=qpow(limit,mod-2); 78 for(int i=0;i<=limit;++i)a[i]=a[i]*g%mod; 79 ll w=qpow(2,n),GG=qpow(2,mod-2); 80 for(int i=0;i<=min(k,n);++i) 81 { 82 ans=(ans+a[i]*w)%mod; 83 w=w*(n-i)%mod*GG%mod; 84 } 85 ans=ans*(n+1)%mod*qpow(2,n*(n-1)/2)%mod; 86 cout<<ans<<endl; 87 return 0; 88 }
https://www.lydsy.com/JudgeOnline/problem.php?id=5093
4.求树上所有点对距离的k次方和，点对中两个点不一样。n<=5E4,k<=150。
还是展开。
1 // luogu-judger-enable-o2 2 #include<bits/stdc++.h> 3 using namespace std; 4 typedef long long int ll; 5 const ll mod=10007; 6 ll n,f1[50005][255],f2[50005][255],head[50005*2],size,x,y,s[555][555],fac[555],k,tmp[555]; 7 struct edge{int to,next;}E[50005*2]; 8 void add(int u,int v) 9 { 10 E[++size].to=v; 11 E[size].next=head[u]; 12 head[u]=size; 13 } 14 void init() 15 { 16 s[0][0]=1; 17 for(int i=1;i<=200;++i) 18 for(int j=1;j<=200;++j) 19 s[i][j]=(j*s[i-1][j]%mod+s[i-1][j-1])%mod; 20 fac[0]=1; 21 for(int i=1;i<=200;++i)fac[i]=fac[i-1]*i%mod; 22 } 23 void dfs1(int u,int F) 24 { 25 f1[u][0]=1; 26 for(int e=head[u];e;e=E[e].next) 27 { 28 int v=E[e].to; 29 if(v==F)continue; 30 dfs1(v,u); 31 for(int i=k;i>=1;--i)f1[u][i]=(f1[u][i]+f1[v][i]+f1[v][i-1])%mod; 32 f1[u][0]=(f1[u][0]+f1[v][0])%mod; 33 } 34 } 35 void dfs2(int u,int F) 36 { 37 for(int i=0;i<=k;++i)f2[u][i]=f1[u][i]; 38 if(F!=0) 39 { 40 for(int i=1;i<=k;++i) 41 tmp[i]=(f2[F][i]-f1[u][i]-f1[u][i-1]+mod*2)%mod; 42 tmp[0]=(f2[F][0]-f1[u][0]+mod)%mod; 43 for(int i=1;i<=k;++i) 44 f2[u][i]=(f2[u][i]+tmp[i]+tmp[i-1])%mod; 45 f2[u][0]=(f2[u][0]+tmp[0])%mod; 46 } 47 for(int e=head[u];e;e=E[e].next) 48 { 49 int v=E[e].to; 50 if(v==F)continue; 51 dfs2(v,u); 52 } 53 } 54 int main() 55 { 56 // freopen("a.in","r",stdin); 57 ios::sync_with_stdio(false); 58 cin>>n>>k; 59 for(int i=1;i<=n-1;++i) 60 { 61 cin>>x>>y; 62 add(x,y); 63 add(y,x); 64 } 65 init(); 66 dfs1(1,0); 67 dfs2(1,0); 68 for(int u=1;u<=n;++u) 69 { 70 int ans=0; 71 for(int i=0;i<=k;++i)ans=(ans+s[k][i]*fac[i]%mod*f2[u][i]%mod)%mod; 72 cout<<ans<<endl; 73 } 74 return 0; 75 }
https://www.luogu.org/problemnew/show/P4827