斯特林数学习笔记
好像全忘了，那就再开一个博客吧。。
第一类斯特林数
\(\begin{bmatrix} n\\m\end{bmatrix}\)表示将\(n\)个不同元素构成\(m\)个圆排列的方案数。
递推式为\(\begin{bmatrix} n\\m\end{bmatrix}=\begin{bmatrix} n-1\\m-1\end{bmatrix}+\begin{bmatrix} n-1\\m\end{bmatrix}(n-1)\)
\(\begin{bmatrix} n-1\\m-1\end{bmatrix}\)就是第\(n\)个元素新开一个圆排列，否则将第\(n\)个元素插进前面任意一个元素的后面。
有两个式子：
\(x^{\overline n}=\sum_{i=0}^nx^i\begin{bmatrix} n\\i\end{bmatrix}\)
\(x^{\underline n}=\sum_{i=0}^n(-1)^{n-i}x^i\begin{bmatrix} n\\i\end{bmatrix}\)
求一行
https://www.luogu.org/problemnew/show/P5408
套上面的上升幂式子，可以直接分治ntt，\(O(n\log^2n)\)
还有一种类似快速幂的做法：
用类似快速幂来推这玩意儿，要支持快速做两种操作：
\(x^{\overline n}\rightarrow x^{\overline{n+1}}\)
\(x^{\overline n}\rightarrow x^{\overline{2n}}\)
第一个直接暴力乘上去就好了。
第二个列一下式子：
\(x^{\overline{2n}}=x^{\overline{n}}(x+n)^{\overline{n}}\)，那么要求\((x+n)^{\overline{n}}\)
套上面的式子：\((x+n)^{\overline{n}}=\sum_{i=0}^n(x+n)^i\begin{bmatrix} n\\i\end{bmatrix}\)
\((x+n)^{\overline{n}}=\sum_{i=0}^n\begin{bmatrix} n\\i\end{bmatrix}\sum_{j=0}^ix^jn^{i-j}C_{i}^{j}\)
化成NTT可以接受的形式：
\(=\sum_{j=0}^n\frac{x^j}{n^jj!}\sum_{i=j}^n\begin{bmatrix} n\\i\end{bmatrix}n^ii! \cdot \frac{1}{(i-j)!}\)
因为你已经求出来了\(n\)以内的斯特林数，所以可以求出\((x+n)^{\overline{n}}\)，然后乘起来可以得到\(x^{\overline{2n}}\)。
#include<bits/stdc++.h>
#define il inline
#define vd void
#define mod 167772161
#define poly std::vector<int>
typedef long long ll;
il ll gi(){
ll x=0,f=1;
char ch=getchar();
while(!isdigit(ch))f^=ch=='-',ch=getchar();
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
il int pow(int x,int y){
int ret=1;
while(y){
if(y&1)ret=1ll*ret*x%mod;
x=1ll*x*x%mod;y>>=1;
}
return ret;
}
#define maxn 524289
poly pA,pB;
int rev[maxn],_lstN,P[maxn],iP[maxn];
il vd ntt(int*A,int N,int t){
for(int i=0;i<N;++i)if(rev[i]>i)std::swap(A[i],A[rev[i]]);
for(int o=1;o<N;o<<=1){
int W=t?P[o]:iP[o];
for(int*p=A;p!=A+N;p+=o<<1)
for(int i=0,w=1;i<o;++i,w=1ll*w*W%mod){
int t=1ll*w*p[i+o]%mod;
p[i+o]=(p[i]-t+mod)%mod;p[i]=(p[i]+t)%mod;
}
}
if(!t){
int inv=pow(N,mod-2);
for(int i=0;i<N;++i)A[i]=1ll*A[i]*inv%mod;
}
}
int N,lg;
il vd setN(int n){
N=1,lg=0;
while(N<n)N<<=1,++lg;
if(N!=_lstN)for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
}
il vd ntt(poly&a,int t){
static int A[maxn];
for(int i=0;i<a.size();++i)A[i]=a[i];memset(A+a.size(),0,4*(N-a.size()));
ntt(A,N,t);
a.resize(N);
for(int i=0;i<N;++i)a[i]=A[i];
int s=a.size();while(s&&!a[s-1])--s;
a.resize(s);
}
il poly mul(poly a,poly b,int newn=-1){
if(newn==-1)newn=a.size()+b.size()-1;
setN(a.size()+b.size()-1);
ntt(a,1),ntt(b,1);
for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod;
ntt(a,0);a.resize(newn);
return a;
}
il poly operator+(poly a,const poly&b){
if(a.size()<b.size())a.resize(b.size());
for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]+b[i])%mod;
return a;
}
il poly operator-(poly a,const poly&b){
if(a.size()<b.size())a.resize(b.size());
for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]-b[i]+mod)%mod;
return a;
}
il poly operator*(poly a,int b){
for(auto&i:a)i=1ll*i*b%mod;
return a;
}
il poly qiudao(poly a){
for(int i=0;i<a.size()-1;++i)a[i]=1ll*a[i+1]*(i+1)%mod;
a.erase(a.end()-1);
return a;
}
il poly jifen(poly a){
a.insert(a.begin(),0);
for(int i=1;i<a.size();++i)a[i]=1ll*a[i]*pow(i,mod-2)%mod;
return a;
}
il poly getinv(poly a){
if(a.size()==1)return poly(1,pow(a[0],mod-2));
int n=a.size(),m=a.size()+1>>1;
poly _a(m);
for(int i=0;i<m;++i)_a[i]=a[i];
poly b=getinv(_a);
setN(n+m*2-2);
ntt(a,1);ntt(b,1);
for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod*b[i]%mod;
ntt(a,0),ntt(b,0);
a.resize(n);
return b*2-a;
}
il poly getln(poly a,int n=-1){
if(n==-1)n=a.size();
a.resize(n);
return jifen(mul(qiudao(a),getinv(a),n));
}
il poly getexp(poly a){
if(a.size()==1)return a[0]=1,a;
int n=a.size(),m=a.size()+1>>1;
poly _a(m);
for(int i=0;i<m;++i)_a[i]=a[i];
poly b=getexp(_a);
return mul(b,poly(1,1)-getln(b,a.size())+a,a.size());
}
il poly operator^(poly a,int b){
int n=a.size();
a=getexp(getln(a)*b);a.resize(n);
return a;
}
il poly sqrt(poly a){
if(a.size()==1)return a;
int n=a.size(),m=a.size()+1>>1;
poly _a(m);
for(int i=0;i<m;++i)_a[i]=a[i];
poly b=sqrt(_a);b.resize(n);
return (b+mul(a,getinv(b),n))*(mod+1>>1);
}
il vd poly_init(){
int G=3,iG=pow(G,mod-2);
for(int i=1;i<maxn;i<<=1)P[i]=pow(G,(mod-1)/(i<<1)),iP[i]=pow(iG,(mod-1)/(i<<1));
}
int fact[262147],ifact[262147];
il poly get_stelin(int n){
if(n==1){poly r(2);r[1]=1;return r;}
int m=n>>1;
poly g=get_stelin(m),_g(g),_h(m+1),h(m+1);
for(int i=0,pm=1;i<=m;++i,pm=1ll*pm*m%mod)_g[i]=1ll*_g[i]*fact[i]%mod*pm%mod;
for(int i=0;i<=m;++i)_h[i]=ifact[m-i];
_h=mul(_h,_g);
int invm=pow(m,mod-2);
for(int i=0,pm=1;i<=m;++i,pm=1ll*pm*invm%mod)h[i]=1ll*_h[m+i]*pm%mod*ifact[i]%mod;
g=mul(h,g,n+1);
if(n&1){
g.resize(n+1);
for(int i=n;i;--i)g[i]=(g[i-1]+1ll*(n-1)*g[i])%mod;
g[0]=1ll*g[0]*(n-1)%mod;
}
return g;
}
int main(){
#ifdef XZZSB
freopen("in.in","r",stdin);
freopen("out.out","w",stdout);
#endif
poly_init();
int n=gi();
fact[0]=1;for(int i=1;i<=n;++i)fact[i]=1ll*fact[i-1]*i%mod;
ifact[n]=pow(fact[n],mod-2);for(int i=n-1;~i;--i)ifact[i]=1ll*ifact[i+1]*(i+1)%mod;
poly f=get_stelin(n);
for(int i=0;i<=n;++i)printf("%d ",f[i]);
return 0;
}求一列
https://www.luogu.org/problemnew/show/P5409
这东西太神仙了
考虑用两种方式拆开\((x+1)^t\)。
第一种就是快速幂的方法，\((x+1)^t=\exp(t\ln(x+1))=\sum_{i=0}^{\infty}\frac{t^i\ln(x+1)^i}{i!}\)
第二种是拆二项式，\((x+1)^t=\sum_{i=0}^{\infty}x^iC_{t}^i\)
\((x+1)^t=\sum_{i=0}^{\infty}\frac{x^i}{i!}t^{\underline i}\)
拆下降幂
\((x+1)^t=\sum_{i=0}^{\infty}\frac{x^i}{i!}\sum_{j=0}^{\infty}(-1)^{i-j}\begin{bmatrix}i\\j\end{bmatrix}t^j\)
\((x+1)^t=\sum_{j=0}^{\infty}t^j\sum_{i=0}^{\infty}\frac{x^i}{i!}(-1)^{i-j}\begin{bmatrix}i\\j\end{bmatrix}\)
将两种拆法放一起
\(\sum_{i=0}^{\infty}\frac{t^i\ln(x+1)^i}{i!}=\sum_{i=0}^{\infty}t^i\sum_{j=0}^{\infty}\frac{x^j}{j!}(-1)^{j-i}\begin{bmatrix}j\\i\end{bmatrix}\)
\(\frac{t^k\ln(x+1)^k}{k!}=\sum_{j=0}^{\infty}\frac{x^j}{j!}(-1)^{j-k}\begin{bmatrix}j\\i\end{bmatrix}\)
因为\(t\)是变量，所以任意取一个\(i\)两边都是相等的，\(i\)取\(k\)就可以求得第\(k\)列的生成函数，\(\mod x^{n+1}\)可以求得0到n的\(\begin{bmatrix}i\\k\end{bmatrix}\)
#pragma GCC optimize("O3,Ofast,no-stack-protector,unroll-loops,fast-math")
#include<bits/stdc++.h>
#define il inline
#define vd void
typedef long long ll;
il ll gi(){
ll x=0,f=1;
char ch=getchar();
while(!isdigit(ch))f^=ch=='-',ch=getchar();
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
//================多项式板子部分===================
#define mod 167772161
#define maxn 524289
#define Gmod 3
#define poly std::vector<int>
il int pow(int x,int y){
int ret=1;
while(y){
if(y&1)ret=1ll*ret*x%mod;
x=1ll*x*x%mod;y>>=1;
}
return ret;
}
int rev[maxn],_lstN,P[maxn],iP[maxn];
il vd ntt(int*A,int N,int t){
for(int i=0;i<N;++i)if(rev[i]>i)std::swap(A[i],A[rev[i]]);
for(int o=1;o<N;o<<=1){
int W=t?P[o]:iP[o];
for(int*p=A;p!=A+N;p+=o<<1)
for(int i=0,w=1;i<o;++i,w=1ll*w*W%mod){
int t=1ll*w*p[i+o]%mod;
p[i+o]=(p[i]-t+mod)%mod;p[i]=(p[i]+t)%mod;
}
}
if(!t){
int inv=pow(N,mod-2);
for(int i=0;i<N;++i)A[i]=1ll*A[i]*inv%mod;
}
}
int N,lg;
il vd setN(int n){
N=1,lg=0;
while(N<n)N<<=1,++lg;
if(N!=_lstN){_lstN=N;for(int i=0;i<N;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);}
}
il vd ntt(poly&a,int t){
static int A[maxn];
for(int i=0;i<a.size();++i)A[i]=a[i];memset(A+a.size(),0,4*(N-a.size()));
ntt(A,N,t);
a.resize(N);
for(int i=0;i<N;++i)a[i]=A[i];
int s=a.size();while(s&&!a[s-1])--s;
a.resize(s);
}
il poly mul(poly a,poly b,int newn=-1){
if(newn==-1)newn=a.size()+b.size()-1;
setN(a.size()+b.size()-1);
ntt(a,1),ntt(b,1);
for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod;
ntt(a,0);a.resize(newn);
return a;
}
il poly operator+(poly a,const poly&b){
if(a.size()<b.size())a.resize(b.size());
for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]+b[i])%mod;
return a;
}
il poly operator-(poly a,const poly&b){
if(a.size()<b.size())a.resize(b.size());
for(int i=0;i<a.size();++i)if(i<b.size())a[i]=(a[i]-b[i]+mod)%mod;
return a;
}
il poly operator*(poly a,int b){
for(auto&i:a)i=1ll*i*b%mod;
return a;
}
il poly qiudao(poly a){
for(int i=0;i<a.size()-1;++i)a[i]=1ll*a[i+1]*(i+1)%mod;
a.erase(a.end()-1);
return a;
}
il poly jifen(poly a){
a.insert(a.begin(),0);
for(int i=1;i<a.size();++i)a[i]=1ll*a[i]*pow(i,mod-2)%mod;
return a;
}
il poly getinv(poly a){
if(a.size()==1)return poly(1,pow(a[0],mod-2));
int n=a.size(),m=a.size()+1>>1;
poly _a(m);
for(int i=0;i<m;++i)_a[i]=a[i];
poly b=getinv(_a);
setN(n+m*2-2);
ntt(a,1);ntt(b,1);
for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%mod*b[i]%mod;
ntt(a,0),ntt(b,0);
a.resize(n);
return b*2-a;
}
il poly getln(poly a,int n=-1){
if(n==-1)n=a.size();
a.resize(n);
return jifen(mul(qiudao(a),getinv(a),n));
}
il poly getexp(poly a){
if(a.size()==1)return a[0]=1,a;
int n=a.size(),m=a.size()+1>>1;
poly _a(m);
for(int i=0;i<m;++i)_a[i]=a[i];
poly b=getexp(_a);
return mul(b,poly(1,1)-getln(b,a.size())+a,n);
}
il poly Pow(poly a,int b,int n){
int cnt0=0;while(!a[cnt0]&&cnt0<a.size())++cnt0;
a.erase(a.begin(),a.begin()+cnt0);
if(1ll*b*cnt0>=n)return poly(n,0);
int m=a.size(),a0=a[0],ia0=pow(a[0],mod-2);
for(int i=0;i<m;++i)a[i]=1ll*ia0*a[i]%mod;
a=getexp(getln(a)*b);a.resize(n);
poly ret(n,0);int p=1ll*cnt0*b;
for(int i=0;i<m;++i){
ret[p++]=1ll*a[i]*a0%mod;
if(p==n)break;
}
return ret;
}
il poly operator%(poly a,poly b){
int n=a.size(),m=b.size();
if(n<m)return a;
std::reverse(a.begin(),a.end());
std::reverse(b.begin(),b.end());
b.resize(n);
poly c=mul(a,getinv(b),n-m+1);
std::reverse(a.begin(),a.end());
b.resize(m);std::reverse(b.begin(),b.end());
std::reverse(c.begin(),c.end());
a=a-mul(b,c);
int s=a.size();while(s&&!a[s-1])--s;
a.resize(s);
return a;
}
namespace sqrt_mod{
struct numb{int real,imag;};
int a,w;
il numb operator*(const numb&a,const numb&b){return(numb){(1ll*a.real*b.real+1ll*a.imag*b.imag%mod*w)%mod,(1ll*a.real*b.imag+1ll*a.imag*b.real)%mod};}
il int sqrt(int x){
while(1){
a=rand()%(mod-1)+1;w=(1ll*a*a-x+mod)%mod;
if(pow(w,mod-1>>1)!=1)break;
}
numb s={a,1},ret={1,0};int y=mod+1>>1;
while(y){
if(y&1)ret=ret*s;
s=s*s;y>>=1;
}
return std::min(mod-ret.real,ret.real);
}
}
il poly sqrt(poly a){
if(a.size()==1)return a[0]=sqrt_mod::sqrt(a[0]),a;
int n=a.size(),m=a.size()+1>>1;
poly _a(m);
for(int i=0;i<m;++i)_a[i]=a[i];
poly b=sqrt(_a);b.resize(n);
return (b+mul(a,getinv(b),n))*(mod+1>>1);
}
il vd poly_init(){
int G=Gmod,iG=pow(G,mod-2);
for(int i=1;i<maxn;i<<=1)P[i]=pow(G,(mod-1)/(i<<1)),iP[i]=pow(iG,(mod-1)/(i<<1));
}
struct _poly_auto_init{_poly_auto_init(){poly_init();}}_auto_init;
//End==============多项式板子部分===================
int fact[131113],ifact[131113];
int main(){
#ifdef XZZSB
freopen("in.in","r",stdin);
freopen("out.out","w",stdout);
#endif
int n=gi(),k=gi();
if(k>n){++n;while(n--)printf("0 ");return 0;}
fact[0]=1;for(int i=1;i<=n;++i)fact[i]=1ll*fact[i-1]*i%mod;
ifact[n]=pow(fact[n],mod-2);for(int i=n;i;--i)ifact[i-1]=1ll*ifact[i]*i%mod;
poly f(2,1);
f=Pow(getln(f,n-k+1),k,n+1);
for(int i=0;i<k;++i)printf("0 ");
for(int i=k;i<=n;++i)printf("%d ",1ll*f[i]*ifact[k]%mod*fact[i]%mod*(i-k&1?mod-1:1)%mod);
return 0;
}第二类斯特林数
\(\begin{Bmatrix} n\\m\end{Bmatrix}\)表示将\(n\)个不同元素放进\(m\)个相同盒子的方案数。
这个东西好像比第一类水一点，题也多一点？
递推式：\(\begin{Bmatrix} n\\m\end{Bmatrix}=\begin{Bmatrix} n-1\\m\end{Bmatrix}m+\begin{Bmatrix} n-1\\m-1\end{Bmatrix}\)
通项：容斥枚举放了\(i\)个空盒子。\(\sum_{i=0}^m\frac{(-1)^i}{i!}\frac{(m-i)^n}{(m-i)!}\)
性质：分解\(m^n\)，枚举有一个盒子非空。\(m^n=\sum_{i=0}^mC_m^i\begin{Bmatrix} n\\i\end{Bmatrix}\)
求一行
直接用通项ntt
求一列
https://www.luogu.org/problemnew/show/P5396
构造生成函数\(F_m(x)=\sum_{i=m}^{\infty}\begin{Bmatrix} i\\m\end{Bmatrix}x^i\)
根据递推式可得\(F_m(x)=F_{m-1}(x)\cdot x+F_m(x)\cdot mx\)
移项得\(F_m(x)=\frac{F_{m-1}(x)}{1-mx}x\)
显然\(F(0)=1\)，所以可以一路转移过去，最后写出来是这样的
\(F_m(x)=\frac{x^m}{(1-x)(1-2x)\cdots(1-mx)}\)
可以分治ntt解决。
具体数学上说\(\begin{bmatrix} n\\m\end{bmatrix}=\begin{Bmatrix} -m\\-n\end{Bmatrix}\)，那么应该第2类的列和第1类的行有点关系，然而好像不可以用那种倍增方法优化，因为\((1-x)\cdots(1-nx)\)并不能修改x推到\([1-(n+1)x]\cdots(1-2nx)\)
好像只能分治ntt了233