前幾天在BZ上的考試考到有關第二類斯特林數的東西
雖然說那道題目到最后並不需要用這個東西來化簡把
不過抱着學習的態度還是學了學有關第二類斯特林數的東西
第二類斯特林數S(n,m)定義為把n個元素划分成m個無序集合的方案數
根據這個定義我們不難寫出遞推式
設狀態S(i,j)，討論第i個元素是否單獨一個集合
若單獨一個集合，則方案數等價於S(i-1,j-1)
若不是單獨一個集合，則他可以在之前任意j個集合里，方案為S(i-1,j)*j
這樣我們得到式子S(i,j)=S(i-1,j-1)+S(i-1,j)*j
遞推的時間復雜度是O(n^2)的
根據定義我們也不難寫出通項表達式
假設集合沒有非空的限制，則答案顯然是m^n
之后我們可以利用容斥原理，枚舉至少有幾個集合是空的，在套上容斥系數我們有(注意集合是無序的，所以最后要消序)
S(n,m)=1/m!*sigma((-1)^k*C(m,k)*(m-k)^n))
直接根據這個式子暴力求解某一行的斯特林數的復雜度是O(n^2)的
但顯然這是一個卷積形式，利用FFT我們可以將時間復雜度優化成O(nlogn)
我們再來討論一下貝爾數，貝爾數的定義是把n個元素划分成若干個無序集合的方案數
不難發現f(n)=sigma(S(n,m)),直接求解某一項貝爾數利用上面的通項公式我們可以做到O(nlogn)
但是我們可以做到更好，考慮第一個元素所在集合的大小，我們不難得到如下遞推式
f(i)=sigma(C(i-1,k-1)*f(i-k))
我們對於這個式子利用CDQ+FFT即可在O(nlog^2n)的時間復雜度內求解前n項的貝爾數
但是我們還可以做到更好
我們考慮每個集合是由若干個元素組成的且集合非空，集合中元素無序
則可以得到集合的生成函數為g(x)=e^x-1
而貝爾數分解成若干個集合的划分方案,則我們有f=e^g
我們構造多項式g，之后多項式求exp即可
時間復雜度O(nlogn)
有意思的是我們可以考慮把n個元素划分成若干個有序集合的方案數，設為多項式f
考慮第一個集合的大小，我們有
f(i)=sigma(C(i,k)*f(i-k))
我們一樣可以用CDQ+FFT來求解
但是我們利用之前在多項式求逆學習總結中提過的技巧進行化簡，就可以化簡成一個多項式求逆的式子
時間復雜度O(nlogn)
我們現在討論一些跟斯特林數相關的化簡技巧
其中最常用的莫過於i^k和S(k,j)之間的轉化了
我們通過第二類斯特林數的通項公式不難發現S(k,j)可以轉化成若干個帶系數的i^k的和
而我們又存在公式i^k=sigma(S(k,j)*j!*C(i,j))
不難發現這里跟i有關的項其實只有C(i,j),我們就把冪通過第二類斯特林數轉化成了組合數
更有趣的是，這個式子后面j!*C(i,j)實際上暴力計算的時間復雜度是O(j)的 QAQ
即計算i的下降階乘冪
我們來看幾道題目QAQ
BZ的題目由於某些奇怪的原因貌似不能放出來QAQ
首先是Crash的文明世界 cojs 1888
這是道很有趣的題目，首先我們如果對(i+1)^k暴力用二項式定理展開會得到一個時間復雜度為O(n*k^2)的DP方程
注意到這里狀態是O(n*k)的，而轉移是O(k)的
由於k很小，我們考慮運用第二類斯特林數將i^k轉化為C(i,j)
由於C(i,j)=C(i-1,j)+C(i-1,j-1)，我們利用這個式子可以做到轉移O(1)，這樣時間復雜度就是O(n*k)的
具體做法是設f(i,j)表示第i個點的C(dis(i,o),j)的和
然后從底向上做一遍DP，然后從上向下在做一遍DP
求出第二類斯特林數還原答案即可
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
using namespace std;
typedef long long LL;
const int maxn=50010;
const int mod=10007;
int n,k,L,u,v;
int now,A,B,Q;
int h[maxn],cnt=0;
struct edge{
int to,next;
}G[maxn<<1];
int f[maxn][151];
int jc[151],inv[151];
int S[151][151];
void add(int x,int y){
++cnt;G[cnt].to=y;G[cnt].next=h[x];h[x]=cnt;
}
void Get_read(){
scanf("%d%d%d",&n,&k,&L);
scanf("%d%d%d%d",&now,&A,&B,&Q);
for(int i=1;i<n;++i){
now=(now*A+B)%Q;
int tmp=((i<L)?i:L);
u=i-now%tmp;v=i+1;
add(u,v);add(v,u);
}
}
int pow_mod(int v,int p){
int tmp=1;
while(p){
if(p&1)tmp=tmp*v%mod;
v=v*v%mod;p>>=1;
}return tmp;
}
void Get_S(){
S[0][0]=1;
for(int i=1;i<=k;++i){
S[i][0]=0;
for(int j=1;j<=i;++j)S[i][j]=(S[i-1][j-1]+S[i-1][j]*j)%mod;
}
jc[0]=1;
for(int i=1;i<=k;++i)jc[i]=jc[i-1]*i%mod;
inv[k]=pow_mod(jc[k],mod-2);
for(int i=k-1;i>=0;--i)inv[i]=inv[i+1]*(i+1)%mod;
}
void Get_DFS(int u,int fa){
f[u][0]=1;
for(int i=h[u];i;i=G[i].next){
int v=G[i].to;
if(v==fa)continue;
Get_DFS(v,u);
f[u][0]=(f[u][0]+f[v][0])%mod;
for(int j=1;j<=k;++j)f[u][j]=(f[u][j]+f[v][j]+f[v][j-1])%mod;
}return;
}
void Get_DP(int u,int fa){
for(int i=h[u];i;i=G[i].next){
int v=G[i].to;
if(v==fa)continue;
for(int j=k;j>=2;--j){
f[v][j]=(f[u][j]+f[u][j-1]-2*f[v][j-1]-f[v][j-2])%mod;
if(f[v][j]<0)f[v][j]+=mod;
}
f[v][1]=(f[u][0]+f[u][1]-2*f[v][0])%mod;
if(f[v][1]<0)f[v][1]+=mod;
f[v][0]=f[u][0];
Get_DP(v,u);
}return;
}
int main(){
freopen("civilization.in","r",stdin);
freopen("civilization.out","w",stdout);
Get_read();Get_S();
Get_DFS(1,-1);
Get_DP(1,-1);
for(int i=1;i<=n;++i){
int ans=0;
for(int j=0;j<=k;++j){
ans=ans+S[k][j]*jc[j]%mod*f[i][j]%mod;
ans%=mod;
}printf("%d\n",ans);
}return 0;
}
cojs 2359 QAQ的圖論題
QAQ 題解在blog里，寫的非常詳細 QAQ
具體做法是考慮每個點的貢獻寫出O(n)的計算式，之后把i^k用第二類斯特林數表示
化簡之后計算的時間復雜度變成O(k)，然后利用FFT在O(klogk)的時間內求出我們所需要的斯特林數就可以了
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define G 3
using namespace std;
typedef long long LL;
const int mod=998244353;
const int maxn=300010;
int n,k,N,len,ans,now;
int jc[maxn],inv[maxn],rev[maxn];
int A[maxn],B[maxn],w[maxn];
int pow_mod(int v,int p){
int tmp=1;
while(p){
if(p&1)tmp=1LL*tmp*v%mod;
v=1LL*v*v%mod;p>>=1;
}return tmp;
}
void FFT(int *A,int n,int flag){
for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int k=2;k<=n;k<<=1){
int mk=(k>>1);
int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
w[0]=1;
for(int i=1;i<mk;++i)w[i]=1LL*w[i-1]*wn%mod;
for(int i=0;i<n;i+=k){
for(int j=0;j<mk;++j){
int a=i+j,b=i+j+mk;
int x=A[a],y=1LL*A[b]*w[j]%mod;
A[a]=(x+y)%mod;
A[b]=x-y;if(A[b]<0)A[b]+=mod;
}
}
}
if(flag==-1){
int c=pow_mod(n,mod-2);
for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
}return;
}
int main(){
freopen("QAQ_Graph.in","r",stdin);
freopen("QAQ_Graph.out","w",stdout);
scanf("%d%d",&n,&k);
jc[0]=1;
for(int i=1;i<=k;++i)jc[i]=1LL*jc[i-1]*i%mod;
inv[k]=pow_mod(jc[k],mod-2);
for(int i=k-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
for(int i=0;i<=k;++i){
if(i&1)A[i]=mod-inv[i];
else A[i]=inv[i];
}
for(int i=0;i<=k;++i)B[i]=1LL*pow_mod(i,k)*inv[i]%mod;
for(N=1,len=0;N<=k;N<<=1,len++);N<<=1,len++;
for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
FFT(A,N,1);FFT(B,N,1);
for(int i=0;i<N;++i)A[i]=1LL*A[i]*B[i]%mod;
FFT(A,N,-1);now=1;ans=0;
for(int i=0;i<=k;++i){
if(now==0)break;
ans=ans+1LL*A[i]*now%mod*pow_mod(2,n-i-1)%mod;
ans%=mod;now=1LL*now*(n-i-1)%mod;
}
ans=1LL*ans*n%mod;
ans=1LL*ans*pow_mod(2,(1LL*(n-1)*(n-2)/2)%(mod-1))%mod;
printf("%d\n",ans);
return 0;
}
cojs 2272 HEOI2016 sum
QAQ 考試的時候沒有做出來，所以要反反復復的多虐幾遍 QAQ
第一發的做法是CDQ+FFT,第二發的做法是多項式求逆
而現在會了第二類斯特林數的一些性質，我們有了第三種做法
前兩種做法都是設f(i)=sigma(S(i,j)*2^j*j!)
但是我們不妨轉換思路,設f(i)=sigma(S(j,i))
這樣我們求完之后對每一項乘以2^i*i!就可以了
考慮把S(j,i)用通項公式展開
對於S(n,m)我們有S(n,m)=sigma((-1)^k*C(m,k)*(m-k)^n)
不難發現跟n有關的項只有(m-k)^n，提出來求sigma可以用等比數列求和公式化簡掉
之后帶入原式我們發現這是一個裸的卷積形式，直接做FFT就可以了
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
#define G 3
using namespace std;
typedef long long LL;
const int maxn=300010;
const int mod=998244353;
int n,N,len,ans,now;
int jc[maxn],inv[maxn];
int A[maxn],B[maxn],rev[maxn];
int w[maxn];
int pow_mod(int v,int p){
int tmp=1;
while(p){
if(p&1)tmp=1LL*tmp*v%mod;
v=1LL*v*v%mod;p>>=1;
}return tmp;
}
void FFT(int *A,int n,int flag){
for(int i=0;i<n;++i)if(i<rev[i])swap(A[i],A[rev[i]]);
for(int k=2;k<=n;k<<=1){
int mk=(k>>1);
int wn=pow_mod(G,flag==1?(mod-1)/k:(mod-1)-(mod-1)/k);
w[0]=1;
for(int i=1;i<mk;++i)w[i]=1LL*w[i-1]*wn%mod;
for(int i=0;i<n;i+=k){
for(int j=0;j<mk;++j){
int a=i+j,b=i+j+mk;
int x=A[a],y=1LL*A[b]*w[j]%mod;
A[a]=(x+y)%mod;
A[b]=x-y;if(A[b]<0)A[b]+=mod;
}
}
}
if(flag==-1){
int c=pow_mod(n,mod-2);
for(int i=0;i<n;++i)A[i]=1LL*A[i]*c%mod;
}return;
}
int main(){
freopen("heoi2016_sum.in","r",stdin);
freopen("heoi2016_sum.out","w",stdout);
scanf("%d",&n);
jc[0]=1;
for(int i=1;i<=n;++i)jc[i]=1LL*jc[i-1]*i%mod;
inv[n]=pow_mod(jc[n],mod-2);
for(int i=n-1;i>=0;--i)inv[i]=1LL*inv[i+1]*(i+1)%mod;
for(int i=0;i<=n;++i){
if(i&1)A[i]=mod-inv[i];
else A[i]=inv[i];
}B[1]=n;
for(int i=2;i<=n;++i){
B[i]=pow_mod(i,n+1)-i;
if(B[i]<0)B[i]+=mod;
B[i]=1LL*B[i]*pow_mod(i-1,mod-2)%mod;
B[i]=1LL*B[i]*inv[i]%mod;
}
for(N=1,len=0;N<=n;N<<=1,len++);N<<=1,len++;
for(int i=0;i<N;++i)rev[i]=rev[i>>1]>>1|((i&1)<<(len-1));
FFT(A,N,1);FFT(B,N,1);
for(int i=0;i<N;++i)A[i]=1LL*A[i]*B[i]%mod;
FFT(A,N,-1);ans=1;now=2;
for(int i=1;i<=n;++i){
ans=ans+1LL*A[i]*now%mod*jc[i]%mod;
ans%=mod;now=(now<<1)%mod;
}printf("%d\n",ans);
return 0;
}
QAQ 第二類斯特林數還是很有趣的呢
准備在cojs補上一道求貝爾數的題目QAQ