斯特林数、欧拉数学习笔记

时间  2021-03-12
标签 html ios spa code htm blog get string 栏目 HTML 繁體版
原文   原文链接

应用html

第一类斯特林数

定义

\(\begin{bmatrix}n\\m\end{bmatrix}\) 表示将 \(n\) 个元素划分为 \(m\) 个圆排列的方案数，也能够写做 \(s(n,m)\)。ios

性质

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

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

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

一个排列惟一对应一个置换，而一个置换惟一对应一组轮换，blog

因此直接枚举划分为 \(i\) 个轮换的方案数，求和后就是排列总数。get

求法

递推式

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

含义分别是把当前的元素做为一个新的圆排列和把当前元素插入到其它元素的后面。it

第一类斯特林数·行

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

分治去乘能够作到 \(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;
}
相关文章
相关标签/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公众号
   欢迎关注本站公众号,获取更多信息
联系我们 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程