斯特林數、歐拉數學習筆記

應用

第一類斯特林數

定義

\(\begin{bmatrix}n\\m\end{bmatrix}\) 表示將 \(n\) 個元素劃分爲 \(m\) 個圓排列的方案數，也能夠寫做 \(s(n,m)\)。

性質

\(x^{\overline n}=\sum_{i=0}^nx^is(n,i)\)

\(x^{\underline n}=\sum_{i=0}^n(-1)^{n-i}x^is(n,i)\)

\(\begin{aligned}n!=\sum_{i=0}^n s(n,i)\end{aligned}\)

一個排列惟一對應一個置換，而一個置換惟一對應一組輪換，

因此直接枚舉劃分爲 \(i\) 個輪換的方案數，求和後就是排列總數。

求法

遞推式

\(s(n,m)=s(n-1,m-1)+(n-1)s(n-1,m)\)

含義分別是把當前的元素做爲一個新的圓排列和把當前元素插入到其它元素的後面。

第一類斯特林數·行

\(x^{\overline n}=\prod\limits_{i=0}^{n-1}(x+i)=\sum\limits_{i=0}^n\begin{bmatrix} n\\i\end{bmatrix}x^i\)

分治去乘能夠作到 \(n log^2 n\)，更優秀的作法是倍增。

由於 \(x^{\overline {2n}}=x^{\overline n} (x+n)^{\overline n}\)，

因此咱們須要由 \(f(x)\) 的值快速推出 \(f(x+d)\) 的值。

\(\begin{aligned} f(x+d) &=\sum_{i=0}^d a_i (x+d)^i \\ &=\sum_{i=0}^d a_i \sum_{j=0}^i C_i^j x^j d^{i-j} \\ &=\sum_{j=0}^d \sum_{i=j}^d a_iC_i^jx^j d^{i-j} \\ &=\sum_{j=0}^d \frac{x^j}{j!} \sum_{i=j}^d a_ii! \frac{d^{i-j}}{(i-j)!} \end{aligned}\)

後面的部分就是一個減法卷積的形式。

最終的時間複雜度爲 \(O(nlogn)\)。

代碼

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=524288,G=3,mod=167772161;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
inline int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[23][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int n,jc[maxn],jcc[maxn],ny[maxn];
void pre(){
	ny[1]=1;
	for(rg int i=2;i<=n;i++) ny[i]=mulmod(mod-mod/i,ny[mod%i]);
	jc[0]=jcc[0]=1;
	for(rg int i=1;i<=n;i++) jc[i]=mulmod(jc[i-1],i),jcc[i]=mulmod(jcc[i-1],ny[i]);
}
int getC(rg int nn,rg int mm){
	return mulmod(jc[nn],mulmod(jcc[mm],jcc[nn-mm]));
}
int a[maxn],b[maxn],c[maxn],d[maxn];
void solve(rg int l,rg int r){
	if(l==r){
		a[0]=0,a[1]=1;
		return;
	}
	rg int mids=(l+r)>>1,len=r-l+1;
	if(len&1) mids--;
	solve(l,mids);
	len=mids+1;
	rg int lim=1,bit=0;
	for(;lim<=len+len;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int i=0,tmp=1;i<=len;i++,tmp=mulmod(tmp,len)) b[i]=mulmod(jcc[i],tmp);
	for(rg int i=0;i<=len;i++) c[i]=mulmod(jc[i],a[i]);
	for(rg int i=len+1;i<lim;i++) b[i]=c[i]=0;
	ntt(c,lim,-1),ntt(b,lim,1);
	for(rg int i=0;i<lim;i++) c[i]=mulmod(c[i],b[i]);
	ntt(c,lim,1);
	for(rg int i=0;i<=len;i++) c[i]=mulmod(c[i],jcc[i]);
	for(rg int i=len+1;i<lim;i++) c[i]=0;
	ntt(a,lim,1),ntt(c,lim,1);
	for(rg int i=0;i<lim;i++) a[i]=mulmod(a[i],c[i]);
	ntt(a,lim,-1);
	for(rg int i=r+2;i<lim;i++) a[i]=0;
	if((r-l+1)&1){
		for(rg int i=0;i<=r;i++) d[i+1]=a[i];
		for(rg int i=0;i<=r;i++) d[i]=addmod(d[i],mulmod(a[i],r));
		for(rg int i=0;i<=r+1;i++) a[i]=d[i];
	}
}
int main(){
	n=read();
	rg int lim=1;
	for(;lim<=n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][i-1],w[t0][1]);
	}
	pre();
	solve(0,n-1);
	for(rg int i=0;i<=n;i++) printf("%d ",a[i]);
	printf("\n");
	return 0;
}

第二類斯特林數

定義

\(\left\{\begin{matrix} n\\k\end{matrix}\right\}\)表示的是把 \(n\) 個不一樣的小球放在 \(m\) 個相同的盒子裏的方案數，也能夠寫做 \(S(n,m)\)。

性質

\(n^k=\sum_ { i=0}^k S(k,i)×i!×C(n,i)\)

左邊的含義就是 \(k\) 個球放在 \(n\) 個盒子裏，

右邊的含義是枚舉有多少個盒子放。

也能夠寫成降低冪的形式：

\(n^k=\sum_{i=0}^k S(k,i) n^{\underline i}\)

求法

遞推式

$S(n,m)=S(n-1,m-1)+mS(n-1,m) $

含義分別是新開了一個集合和在已有的集合中選擇一個放。

容斥

\(S(n,m)=\frac{1}{m!}\sum_{i=0}^m(-1)^iC_m^i(m-i)^n\)

考慮用二項式反演。

假設一共有 \(n\) 個小球，

設 \(f[i]\) 爲剛好有 \(i\) 個盒子放了小球的方案數，

\(g[i]\) 爲至多有 \(i\) 個盒子放了小球的方案數。

那麼 \(g[k]=\sum_{i=0}^k C_k^if[i]\)。

則 \(f[k]=\sum_{i=0}^k(-1)^{k-i}C_k^ig[i]=\sum_{i=0}^k(-1)^{k-i}C_k^i i^n\)，

化得好看一點就成了 \(f[k]=\sum_{i=0}^k(-1)^iC_k^i (k-i)^n\) ，

那麼 \(f[m]=\sum_{i=0}^m(-1)^iC_m^i(m-i)^n\)。

由於咱們算組合數的時候是按照盒子不一樣算的，因此還要除以一個 \(m!\)。

第二類斯特林數·行

容斥的式子能夠化成卷積的形式。

\(S(n,m)=\frac{1}{m!}\sum_{i=0}^m(-1)^iC_m^i(m-i)^n\)

\(S(n,m)=\sum_{i=0}^m\frac{(-1)^i}{i!}\frac{(m-i)^n}{(m-i)!}\)

就能夠用 \(ntt\) 在 \(nlogn\) 的時間複雜度內求出一行的第二類斯特林數

代碼

#include<cstdio>
#include<algorithm>
#include<cstring>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=524289,mod=167772161,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(A[j+k+len],w[t0][k]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int a[maxn],b[maxn],jc[maxn],jcc[maxn],ny[maxn],n;
void pre(){
	ny[1]=1;
	for(rg int i=2;i<=n;i++) ny[i]=mulmod(mod-mod/i,ny[mod%i]);
	jc[0]=jcc[0]=1;
	for(rg int i=1;i<=n;i++) jc[i]=mulmod(jc[i-1],i),jcc[i]=mulmod(jcc[i-1],ny[i]);
}
int main(){
	n=read();
	pre();
	for(rg int i=0;i<=n;i++){
		if(i&1) a[i]=mod-jcc[i];
		else a[i]=jcc[i];
		b[i]=mulmod(ksm(i,n),jcc[i]);
	}
	rg int lim=1,bit=0;
	for(;lim<=n+n;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][1],w[t0][i-1]);
	}
	ntt(a,lim,1),ntt(b,lim,1);
	for(rg int i=0;i<lim;i++) a[i]=mulmod(a[i],b[i]);
	ntt(a,lim,-1);
	for(rg int i=0;i<=n;i++) printf("%d ",a[i]);
	printf("\n");
	return 0;
}

第二類斯特林數·列

設 \(F(n,x)=\sum_{i} \begin{Bmatrix} i \\ n\end{Bmatrix} x^i\)。

則

\(\begin{aligned} F(n,x)&=\sum_{i} \begin{Bmatrix} i \\ n\end{Bmatrix}x^i \\ &=\sum_{i} (n\begin{Bmatrix} i-1 \\ n\end{Bmatrix}+\begin{Bmatrix} i-1 \\ n-1\end{Bmatrix})x^i \\ &=n\sum_{i}\begin{Bmatrix} i-1 \\ n\end{Bmatrix}x^i+\sum_{i} \begin{Bmatrix} i-1 \\ n-1\end{Bmatrix}x^i \\ &=nxF(n,x)+xF(n-1,x)\end{aligned}\)

因此 \(F(n,x)=\frac{x}{1-nx}F(n-1,x)\)。

又由於 \(F(0,x)=1\)，

因此 \(F(n,x)=\frac{x^n}{\prod_{i=1}^n (1-ix)}\)。

把下面的部分分治乘起來再求逆便可。

代碼

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
#include<vector>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=1<<19,mod=167772161,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[25][maxn],wz[maxn];
void ntt(std::vector<int>& A,rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A.begin()+1,A.end());
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int n,k;
std::vector<int> g[maxn];
void dfs(rg int da,rg int l,rg int r){
	if(l==r){
		g[da].resize(2);
		g[da][0]=1,g[da][1]=mod-l;
		return;
	}
	rg int mids=(l+r)>>1;
	dfs(da<<1,l,mids);
	dfs(da<<1|1,mids+1,r);
	rg int lim=1,bit=0,len=r-l+1;
	for(;lim<=len+len;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	g[da].resize(lim),g[da<<1].resize(lim),g[da<<1|1].resize(lim);
	ntt(g[da<<1],lim,1),ntt(g[da<<1|1],lim,1);
	for(rg int i=0;i<lim;i++) g[da][i]=mulmod(g[da<<1][i],g[da<<1|1][i]);
	ntt(g[da],lim,-1);
	for(rg int i=len+1;i<lim;i++) g[da][i]=0;
}
std::vector<int> finv,ans;
void getinv(std::vector<int>&A,std::vector<int>&B,rg int deg){
	if(deg==1){
		B.resize(1);
		B[0]=ksm(A[0],mod-2);
		return;
	}
	getinv(A,B,(deg+1)>>1);
	rg int lim=1,bit=0;
	for(;lim<deg+deg;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	finv.resize(lim);
	if(A.size()<lim) A.resize(lim);
	for(rg int i=0;i<deg;i++) finv[i]=A[i];
	for(rg int i=deg;i<lim;i++) finv[i]=0;
	B.resize(lim);
	ntt(finv,lim,1),ntt(B,lim,1);
	for(rg int i=0;i<lim;i++) B[i]=delmod(mulmod(2,B[i]),mulmod(B[i],mulmod(B[i],finv[i])));
	ntt(B,lim,-1);
	for(rg int i=deg;i<lim;i++) B[i]=0;
}
int main(){
	n=read(),k=read();
	rg int lim=1;
	for(;lim<=n+n;lim<<=1);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][i-1],w[t0][1]);
	}
	dfs(1,1,k);
	if(n<k){
		for(rg int i=0;i<=n;i++) printf("0 ");
		printf("\n");
		return 0;
	}
	getinv(g[1],ans,n-k+2);
	for(rg int i=0;i<k && i<=n;i++) printf("0 ");
	for(rg int i=k;i<=n;i++) printf("%d ",ans[i-k]);
	printf("\n");
	return 0;
}

斯特林反演

\(F_i=\sum_{j=0}^{i}\begin{Bmatrix}i\\ j\end{Bmatrix}G_j\Leftrightarrow G_i=\sum_{j=0}^{i}(-1)^{i-j}\begin{bmatrix}i\\ j\end{bmatrix}F_j\)

\(F_i=\sum_{j=0}^{i}\begin{bmatrix}i\\ j\end{bmatrix}G_j\Leftrightarrow G_i=\sum_{j=0}^{i}(-1)^{i-j}\begin{Bmatrix}i\\ j\end{Bmatrix}F_j\)

歐拉數

定義

歐拉數 \(\langle\begin{smallmatrix}n\\ k\end{smallmatrix}\rangle\) 是 \(\{1,2,...,n\}\) 的有 \(k\) 個升高的排列 \(\pi_1 \pi_2 \cdots \pi_n\) 的個數，也就是說，在其中 \(k\) 個地方有 \(\pi_j<\pi_{j+1}\)

性質

\(\left\langle\begin{matrix}n\\m\end{matrix}\right\rangle=\sum_{k=0}^{m}\binom{n+1}{k}(-1)^k(m+1-k)^n\)

求法

遞推式

\(\left\langle\begin{matrix}n\\m\end{matrix}\right\rangle=(n-m) \left\langle\begin{matrix}n-1 \\ m-1 \end{matrix}\right\rangle+(m+1)\left\langle\begin{matrix}n-1\\m\end{matrix}\right\rangle\)

強制從小到大放數。

若是放在最前面，升高的排列數不會增長，若是放在以前某個升高的排列中間，那麼會減小一個排列增長一個排列，總數不變，不然放在其它的位置升高的排列數都會增長一個。

求一行歐拉數

構造兩個多項式，令

\(f(x)=\sum_k \binom{n+1}{k}(-1)^kx^k\)

\(g(x)=\sum_k (k+1)^nx^k\)

這兩個多項式卷積獲得的結果就是第 \(n\) 行的歐拉數。

代碼

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=1<<19,mod=998244353,G=3;
inline int addmod(rg int now1,rg int now2){
	return now1+=now2,now1>=mod?now1-mod:now1;
}
inline int delmod(rg int now1,rg int now2){
	return now1-=now2,now1<0?now1+mod:now1;
}
inline int mulmod(rg long long now1,rg int now2){
	return now1*=now2,now1>=mod?now1%mod:now1;
}
int ksm(rg int ds,rg int zs){
	rg int nans=1;
	while(zs){
		if(zs&1) nans=mulmod(nans,ds);
		ds=mulmod(ds,ds);
		zs>>=1;
	}
	return nans;
}
int w[21][maxn],wz[maxn];
void ntt(rg int A[],rg int lim,rg int typ){
	for(rg int i=0;i<lim;i++) if(i<wz[i]) std::swap(A[i],A[wz[i]]);
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		for(rg int j=0,now=len<<1;j<lim;j+=now){
			for(rg int k=0;k<len;k++){
				rg int x=A[j+k],y=mulmod(w[t0][k],A[j+k+len]);
				A[j+k]=addmod(x,y),A[j+k+len]=delmod(x,y);
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+lim);
		rg int ny=ksm(lim,mod-2);
		for(rg int i=0;i<lim;i++) A[i]=mulmod(A[i],ny);
	}
}
int ny[maxn],jc[maxn],jcc[maxn];
int getC(rg int nn,rg int mm){
	return mulmod(jc[nn],mulmod(jcc[nn-mm],jcc[mm]));
}
int n,a[maxn],b[maxn];
int main(){
	n=read();
	rg int lim=1,bit=0;
	for(;lim<=n+n;lim<<=1) bit++;
	for(rg int i=0;i<lim;i++) wz[i]=(wz[i>>1]>>1)|((i&1)<<(bit-1));
	for(rg int len=1,t0=0;len<lim;len<<=1,t0++){
		w[t0][0]=1,w[t0][1]=ksm(G,(mod-1)/(len<<1));
		for(rg int i=2;i<len;i++) w[t0][i]=mulmod(w[t0][i-1],w[t0][1]);
	}
	ny[1]=1;
	for(rg int i=2;i<=n+1;i++) ny[i]=mulmod(mod-mod/i,ny[mod%i]);
	jc[0]=jcc[0]=1;
	for(rg int i=1;i<=n+1;i++) jc[i]=mulmod(jc[i-1],i),jcc[i]=mulmod(jcc[i-1],ny[i]);
	for(rg int i=0;i<=n;i++) a[i]=ksm(i+1,n),b[i]=i&1?mod-getC(n+1,i):getC(n+1,i);
	ntt(a,lim,1),ntt(b,lim,1);
	for(rg int i=0;i<lim;i++) a[i]=mulmod(a[i],b[i]);
	ntt(a,lim,-1);
	for(rg int i=0;i<=n;i++) printf("%d ",a[i]);
	printf("\n");
	return 0;
}