斯特林数、欧拉数学习笔记
第一类斯特林数
定义
\(\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_n = \sum_{i=0}^{n}\begin{Bmatrix} n\\ i\end{Bmatrix}g_i\Leftrightarrow g_n = \sum_{i=0}^{n}\begin{bmatrix}n\\ i\end{bmatrix}(-1)^{n-i}f_i\)
\(f_n = \sum_{i=0}^{n}\begin{Bmatrix}n\\ i\end{Bmatrix}(-1)^{n-i}g_i\Leftrightarrow g_n = \sum_{i=0}^{n}\begin{bmatrix}n\\ i\end{bmatrix}f_i\)
\(f_n = \sum_{i=n}^{\infty}\begin{Bmatrix}i\\ n\end{Bmatrix}(-1)^{i-n}g_i\Leftrightarrow g_n = \sum_{i=n}^{\infty}\begin{bmatrix}i\\ n\end{bmatrix}f_i\)
\(f_n = \sum_{i=n}^{\infty}\begin{Bmatrix}i\\ n\end{Bmatrix}g_i\Leftrightarrow g_n = \sum_{i=n}^{\infty}\begin{bmatrix}i\\ n\end{bmatrix}(-1)^{i-n}f_i\)
欧拉数
定义
欧拉数 \(\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;
}
