斯特林數求解

Deep_Kevin發表於2020-10-08

正題

第二類斯特林數.行

根據通項公式有: $S(n,m)=\frac{1}{m!}\sum_{i=0}^mC_m^i(-1)^i(m-i)^n$

拆開來就可以變成卷積的形式.

第一類斯特林數.行:

第一類斯特林數實際上是上升冪展開成普通冪的係數,那麼只需要多項式平移就可以倍增求出 $s(n,0)-s(n,n)$

#include<bits/stdc++.h>
#define vi vector<int>
using namespace std;

const int mod=167772161;
int ksm(int x,int t){
	int tot=1;
	while(t){
		if(t&1) tot=1ll*tot*x%mod;
		x=1ll*x*x%mod;
		t/=2;
	}
	return tot;
}
const int N=1<<19;//1<<17
int where[N],limit;
int poor[N*5],inv[N],n,m,fac[N],finv[N],h[N];
int*w[2][19],*now=poor;
void ad(int&x,int y){x=(x+y>=mod)?(x+y-mod):(x+y);}
void pre(){
	inv[1]=1;for(int i=2;i<=n;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
	fac[0]=1;for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%mod;
	finv[n]=ksm(fac[n],mod-2);for(int i=n-1;i>=0;i--) finv[i]=1ll*finv[i+1]*(i+1)%mod;
	for(int u=0,l=2;u<19;u++,l<<=1){
		w[1][u]=now;now+=l+1;
		w[0][u]=now;now+=l+1;
		w[1][u][0]=1;w[1][u][1]=ksm(3,(mod-1)/l);
		for(int i=2;i<=l;i++) w[1][u][i]=1ll*w[1][u][i-1]*w[1][u][1]%mod;
		for(int i=0;i<=l;i++) w[0][u][i]=w[1][u][l-i];
	}
}
void prepare(int n){
	int l=0;limit=1;
	while(limit<=n) limit<<=1,l++;
	for(int i=0;i<limit;i++) where[i]=(where[i>>1]>>1)|((i&1)<<(l-1));
}
void DFT(vi&now,int op){
	int**w1=w[op];
	for(int i=0;i<limit;i++) if(i<where[i]) swap(now[i],now[where[i]]);
	for(int l=1,u=0,r=2;r<=limit;l=r,r<<=1,u++){
		int*w2=w1[u];
		for(int i=0;i<limit;i+=r){
			for(int v=0;v<l;v++){
				int&x=now[i+v],&y=now[i+l+v],b=1ll*y*w2[v]%mod;
				y=(x>=b)?(x-b):(x+mod-b);
				x=(x+b>=mod)?(x+b-mod):(x+b);
			}
		}
	}
	if(!op){
		int tmp=ksm(limit,mod-2);
		for(int i=0;i<limit;i++) now[i]=1ll*now[i]*tmp%mod;
	}
}
vi operator*(vi f,vi g){
	int tot=f.size()+g.size()-2;
	prepare(tot);f.resize(limit);g.resize(limit);
	DFT(f,1);DFT(g,1);
	for(int i=0;i<limit;i++) f[i]=1ll*f[i]*g[i]%mod;
	DFT(f,0);f.resize(tot+1);
	return f;
}
vi A,B;

void gas(vi&f){
	int n=f.size()-1,m=n/2;
	if(n==1){f[1]=1;return ;}
	vi f0;f0.resize(m+1);gas(f0);
	vi g;g.resize(m+1);
	int t=1;f=f0;
	for(int i=0;i<=m;i++) 
		f0[i]=1ll*fac[i]*f0[i]%mod,g[i]=1ll*t*finv[i]%mod,t=1ll*t*m%mod;
	reverse(g.begin(),g.end());g=f0*g;
	for(int i=0;i<=m;i++) g[i]=1ll*g[m+i]*finv[i]%mod;
	g.resize(m+1);f=f*g;
	if(n%2==1) {
		f.push_back(0);
		for(int i=n-1;i>=0;i--) ad(f[i+1],f[i]),f[i]=1ll*f[i]*(n-1)%mod;
	}
}

int main(){
	scanf("%d",&n);pre();
	A.resize(n+1);gas(A);
	for(int i=0;i<=n;i++) printf("%d ",A[i]);
}

第二類斯特林數.列

將遞推公式寫出來,設 $G_k=\sum_{i=0}S(i,k)x^i$

顯然有 $\\G_k(x)=kxG_k(x)+xG_{k-1}(x) \\G_k(x)=\frac{x}{1-kx}G_{k-1}(x)$

那麼就是要求一個字首的(1-kx)卷積,多項式平移即可,最後討論一下係數.

#include<bits/stdc++.h>
#define vi vector<int>
using namespace std;

const int mod=167772161;
int ksm(int x,int t){
	int tot=1;
	while(t){
		if(t&1) tot=1ll*tot*x%mod;
		x=1ll*x*x%mod;
		t/=2;
	}
	return tot;
}
const int N=1<<19;//1<<17
int where[N],limit;
int poor[N*5],inv[N],n,m,fac[N],finv[N],h[N];
int*w[2][19],*now=poor;
void ad(int&x,int y){x=(x+y>=mod)?(x+y-mod):(x+y);}
void pre(){
	inv[1]=1;for(int i=2;i<=n;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
	fac[0]=1;for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%mod;
	finv[n]=ksm(fac[n],mod-2);for(int i=n-1;i>=0;i--) finv[i]=1ll*finv[i+1]*(i+1)%mod;
	for(int u=0,l=2;u<19;u++,l<<=1){
		w[1][u]=now;now+=l+1;
		w[0][u]=now;now+=l+1;
		w[1][u][0]=1;w[1][u][1]=ksm(3,(mod-1)/l);
		for(int i=2;i<=l;i++) w[1][u][i]=1ll*w[1][u][i-1]*w[1][u][1]%mod;
		for(int i=0;i<=l;i++) w[0][u][i]=w[1][u][l-i];
	}
}
void prepare(int n){
	int l=0;limit=1;
	while(limit<=n) limit<<=1,l++;
	for(int i=0;i<limit;i++) where[i]=(where[i>>1]>>1)|((i&1)<<(l-1));
}
void DFT(vi&now,int op){
	int**w1=w[op];
	for(int i=0;i<limit;i++) if(i<where[i]) swap(now[i],now[where[i]]);
	for(int l=1,u=0,r=2;r<=limit;l=r,r<<=1,u++){
		int*w2=w1[u];
		for(int i=0;i<limit;i+=r){
			for(int v=0;v<l;v++){
				int&x=now[i+v],&y=now[i+l+v],b=1ll*y*w2[v]%mod;
				y=(x>=b)?(x-b):(x+mod-b);
				x=(x+b>=mod)?(x+b-mod):(x+b);
			}
		}
	}
	if(!op){
		int tmp=ksm(limit,mod-2);
		for(int i=0;i<limit;i++) now[i]=1ll*now[i]*tmp%mod;
	}
}
vi operator*(vi f,vi g){
	int tot=f.size()+g.size()-2;
	prepare(tot);f.resize(limit);g.resize(limit);
	DFT(f,1);DFT(g,1);
	for(int i=0;i<limit;i++) f[i]=1ll*f[i]*g[i]%mod;
	DFT(f,0);f.resize(tot+1);
	return f;
}

void INV(vi&f){
	int n=f.size();
	if(n==1) {f[0]=ksm(f[0],mod-2);return ;}
	vi f0=f;f0.resize((n+1)/2);INV(f0);
	prepare(n+2*(n+1)/2-3);f.resize(limit);f0.resize(limit);
	DFT(f,1);DFT(f0,1);
	for(int i=0;i<limit;i++) f[i]=1ll*f0[i]*(mod+2-1ll*f[i]*f0[i]%mod)%mod;
	DFT(f,0);f.resize(n);
}
vi A,B;

void gas(vi&f){
	int n=f.size()-1,m=n/2;
	if(n==1){f[1]=1;return ;}
	vi f0;f0.resize(m+1);gas(f0);
	vi g;g.resize(m+1);
	int t=1;f=f0;
	for(int i=0;i<=m;i++) 
		f0[i]=1ll*fac[i]*f0[i]%mod,g[i]=1ll*t*finv[i]%mod,t=1ll*t*m%mod;
	reverse(g.begin(),g.end());g=f0*g;
	for(int i=0;i<=m;i++) g[i]=1ll*g[m+i]*finv[i]%mod;
	g.resize(m+1);f=f*g;
	if(n%2==1) {
		f.push_back(0);
		for(int i=n-1;i>=0;i--) ad(f[i+1],f[i]),f[i]=1ll*f[i]*(n-1)%mod;
	}
}

int main(){
	scanf("%d %d",&n,&m);pre();
	if(n<m){for(int i=0;i<=n;i++) printf("0 ");return 0;}
	A.resize(m+2);gas(A);
	int t=(m&1)?(mod-1):1;
	for(int i=0;i<=m;i++) A[i]=1ll*A[i+1]*t%mod,t=(t==1)?(mod-1):(1);
	A.pop_back();reverse(A.begin(),A.end());
	A.resize(n-m+1);INV(A);
	for(int i=0;i<m;i++) printf("0 ");
	for(int i=0;i<n-m+1;i++) printf("%d ",A[i]);
}

第一類斯特林數.列

實際上就是求把若干個個元素拆成k個置換環的方案數,顯然可以構造單個置換環的EGF,然後求這個EGF的k次,再除以k的階乘,然後就可以的得到i個元素拆成k個置換環的EGF

#include<bits/stdc++.h>
#define vi vector<int>
using namespace std;
 
const int mod=167772161;
int ksm(int x,int t){
	int tot=1;
	while(t){
		if(t&1) tot=1ll*tot*x%mod;
		x=1ll*x*x%mod;
		t/=2;
	}
	return tot;
}
const int N=3000010,I=ksm(3,(mod-1)/4),Iinv=ksm(I,mod-2);
int limit,inv[N],where[N],fac[N],finv[N],n,m;
int poor[16778000];
int*w[2][22],*now=poor;
void pre(){
	inv[1]=1;for(int i=2;i<=3000000;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
	fac[0]=1;for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%mod;
	finv[n]=ksm(fac[n],mod-2);for(int i=n-1;i>=0;i--) finv[i]=1ll*finv[i+1]*(i+1)%mod;
	for(int l=2,u=0;u<=21;u++,l<<=1){
		w[0][u]=now;now+=l+1;
		w[1][u]=now;now+=l+1;
		w[1][u][0]=1;w[1][u][1]=ksm(3,(mod-1)/l);
		for(int i=2;i<=l;i++) w[1][u][i]=1ll*w[1][u][i-1]*w[1][u][1]%mod;
		for(int i=0;i<=l;i++) w[0][u][i]=w[1][u][l-i];
	}
}
void prepare(int n){
	int l=0;limit=1;
	while(limit<=n) limit<<=1,l++;
	for(int i=0;i<limit;i++) where[i]=(where[i>>1]>>1)|((i&1)<<(l-1));
}
void DFT(vi&now,int op){
	for(int i=0;i<limit;i++) if(i<where[i]) swap(now[i],now[where[i]]);
	for(int l=1,u=0;l<limit;l<<=1,u++){
		for(int i=0;i<limit;i+=l*2)
			for(int v=0;v<l;v++){
				int x=i+v,y=i+l+v,a=now[x],b=1ll*now[y]*w[op][u][v]%mod;
				now[x]=a+b,now[x]>=mod?now[x]-=mod:0;
				now[y]=a+mod-b;now[y]>=mod?now[y]-=mod:0;
			}
	}
	if(!op){
		int tmp=ksm(limit,mod-2);
		for(int i=0;i<limit;i++) now[i]=1ll*now[i]*tmp%mod;
	}
}
vi operator*(vi f,vi g){
	int tot=f.size()+g.size()-2;
	prepare(tot);f.resize(limit);g.resize(limit);
	DFT(f,1);DFT(g,1);
	for(int i=0;i<limit;i++) f[i]=1ll*f[i]*g[i]%mod;
	DFT(f,0);f.resize(tot+1);
	return f;
}
vi operator+(vi f,vi g){
	int n=max(f.size(),g.size());
	f.resize(n);g.resize(n);
	for(int i=0;i<n;i++) f[i]=f[i]+g[i],f[i]>=mod?f[i]-=mod:0;
	return f;
}
vi operator-(vi f,vi g){
	int n=max(f.size(),g.size());
	f.resize(n);g.resize(n);
	for(int i=0;i<n;i++) f[i]=f[i]+mod-g[i],f[i]>=mod?f[i]-=mod:0;
	return f;
}
void INV(vi&f){
	int n=f.size();
	if(n==1) {f[0]=ksm(f[0],mod-2);return ;}
	vi f0=f;f0.resize((n+1)/2);INV(f0);
	prepare(n+2*(n+1)/2-3);f.resize(limit);f0.resize(limit);
	DFT(f,1);DFT(f0,1);
	for(int i=0;i<limit;i++) f[i]=1ll*f0[i]*(mod+2-1ll*f[i]*f0[i]%mod)%mod;
	DFT(f,0);f.resize(n);
}
void DERI(vi&f){for(int i=0;i<f.size()-1;i++) f[i]=1ll*f[i+1]*(i+1)%mod;f.pop_back();}
void INTER(vi&f){f.push_back(0);for(int i=f.size()-1;i>=1;i--) f[i]=1ll*f[i-1]*inv[i]%mod;f[0]=0;}
void LN(vi&f){int n=f.size();vi g=f;DERI(g);INV(f);f=f*g;INTER(f);f.resize(n);}
void EXP(vi&f){
	int n=f.size();
	if(n==1) {f[0]=1;return ;}
	vi f0=f;f0.resize((n+1)/2);EXP(f0);
	vi g=f0;g.resize(n);LN(g);
	f=f0-(g-f)*f0;f.resize(n);
}
void KSM(vi&f,int a,int b,bool full){//k mod , k mod-1 , k >= n
	int n=f.size(),c,d=n;
	for(int i=0;i<n;i++) if(f[i]) {c=f[i],d=i;break;}
	if(d && full || 1ll*d*a>=n){for(int i=0;i<n;i++) f[i]=0;return ;}
	int tmp=ksm(c,mod-2);c=ksm(c,b);
	for(int i=0;i<n-d;i++) f[i]=1ll*f[i+d]*tmp%mod;
	f.resize(n-d);LN(f);
	for(int i=0;i<n-d;i++) f[i]=1ll*f[i]*a%mod;
	EXP(f);f.resize(n);
	for(int i=n-1;i>=d*a;i--) f[i]=1ll*f[i-d*a]*c%mod;
	for(int i=0;i<d*a;i++) f[i]=0;
}
 
int main(){
	scanf("%d %d",&n,&m);pre();
	vi A;A.resize(n+1);A[0]=0;
	for(int i=1;i<=n;i++) A[i]=1ll*fac[i-1]*finv[i]%mod;
	KSM(A,m,m,0);
	for(int i=0;i<=n;i++) printf("%lld ",1ll*A[i]*fac[i]%mod*finv[m]%mod);
}

【3】斯特林數
2024-08-30
斯特林
2024-11-07
卡特蘭數和斯特林數
2024-07-26
第二類斯特林數小記
2024-12-09
下降冪及斯特林數雜談
2024-06-10
路徑規劃演算法 - 求解最短路徑 - Dijkstra(迪傑斯特拉)演算法
2023-12-07
演算法
求解數獨
2020-10-28
cf932E. Team Work(第二類斯特靈數組合數)
2018-09-29
奈奎斯特取樣定理
2024-04-12
迪傑斯特拉與spfa
2020-10-14
利用瓊斯矩陣求解一般偏振問題
2024-08-04
矩陣
大數翻倍法求解CRT
2021-05-05
最短路-迪傑斯特拉(dijkstra)
2024-05-23
Python中的布林值用途及特點介紹！
2023-10-13
Python
面試題：布林變數
2020-03-12
面試題變數
路特斯財報：2024年上半年路特斯總收入達到3.98億美元同比增長206%
2024-08-29
梯度下降法中導數的求解
2021-09-09
梯度
表示式求解（非一位數）
2019-05-14
馬斯克起訴OpenAI和奧特曼
2024-03-02
馬斯克OpenAI奧特曼
離散數學（格與布林代數）
2021-10-31
求區間不同數的個數【主席樹求解】
2018-09-12
奈奎斯特穩定判據（Nyquist stability criterion）
2024-09-04
UI
440億美元，馬斯克成功收購推特
2022-04-26
馬斯克
蘭切斯特方程的分析和案例解釋
2018-12-17
Cplex混合整數規劃求解(Python API)
2023-10-01
PythonAPI
回溯演算法求解橋本分數式
2018-07-02
演算法
伯努利數,求解自然數冪和的關鍵係數
2020-10-08
打造數字化時代的“場景+管理”方案，特斯聯為什麼選擇SAP？
2020-07-01
弗林德斯大學：研究發現近視人群睡眠質量更差
2021-05-28
復旦DMG：2021中國開放數林指數
2021-07-24
求區間不同數的個數【樹狀陣列求解】
2018-09-12
陣列
每據內帶制萬斯入二特七dtb
2022-03-08
心斯採公滿八列合特維安kba
2022-02-28
求最短路徑-----迪傑斯特拉演算法
2020-12-13
演算法
最短路徑—Dijkstra(迪傑斯特拉)演算法
2018-08-01
演算法
更相減損術求解最大公約數
2020-11-07
運籌優化（十）--整數規劃求解
2019-01-17
優化
弗羅斯特：自動駕駛每年將為英國經濟貢獻數百億英鎊
2019-04-07
自動駕駛

斯特林數求解

正題

相關文章