洛谷P5396 【模板】第二类斯特林数·列（生成函数+分治ntt/倍增ntt）

最新推荐文章于 2021-12-01 11:20:41 发布

SC.ldxcaicai

最新推荐文章于 2021-12-01 11:20:41 发布

阅读量562

点赞数 1

CC 4.0 BY-SA版权

分类专栏： # 生成函数 # ntt # 斯特林数

本文链接：https://blog.youkuaiyun.com/dreaming__ldx/article/details/98383385

ntt 同时被 3 个专栏收录

28 篇文章

订阅专栏

生成函数

27 篇文章

订阅专栏

斯特林数

7 篇文章

订阅专栏

博客介绍了如何在O(klog²k)的时间复杂度内计算第二类斯特林数的第k列，通过生成函数和递推关系式展示解题思路。提供了使用分治和倍增NTT算法实现斯特林数计算的方法。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

传送门
给出同阶的 $n, k$
要求在 $O(k\log^2k)$ 时间内处理出第二类斯特林数的第 $k$ 列。
即求出 $S_{i}^{k},0\le i\le n$

思路：
考虑第 $k$ 列斯特林数的生成函数：
$S_k(x)=\sum_{i=0}^{\infin}S_{i}^kx^i=\sum_{i=0}^{\infin}(S_{i-1}{k-1}+S_{i-1}{k}k)x^i=kxS_k(x)+S_{k-1}(x)$
$\Rightarrow S_k(x)=\frac{x^k}{\prod_{i=1}^k(1-ix)}=\frac{x^k}{\prod_{i=1}^k(x-i)}$
那么可以用 $O(n\log^2n)$ 的分治 $n t t$ 或者 $O(n\log n)$ 的倍增 $n t t$ 算分母
分治 $n t t$ ：

#include<bits/stdc++.h>
#define ri register int
using namespace std;
typedef long long ll;
const int mod=167772161;
inline int add(const int&a,const int&b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return a>=b?a-b:a-b+mod;}
inline int mul(const int&a,const int&b){return (ll)a*b%mod;}
inline void Add(int&a,const int&b){a=a+b>=mod?a+b-mod:a+b;}
inline void Dec(int&a,const int&b){a=a>=b?a-b:a-b+mod;}
inline void Mul(int&a,const int&b){a=(ll)a*b%mod;}
inline int ksm(int a,int p){int ret=1;for(;p;p>>=1,a=mul(a,a))if(p&1)Mul(ret,a);return ret;}
vector<int>w[20],pos[20];
int invv[20];
const int Lim=1<<20;
inline void ntt_init(){
	for(ri mt=1,inv=(mod+1)>>1,W,i=1,t=0;i<Lim;i<<=1,++t,Mul(mt,inv)){
		invv[t]=mt;
		pos[t].resize(i),w[t].resize(i),w[t][0]=1,W=ksm(3,(mod-1)>>(t+1));
		for(ri j=1;j<i;++j)w[t][j]=mul(w[t][j-1],W),pos[t][j]=(pos[t][j>>1]>>1)|((j&1)<<(t-1));
	}
}
int lim,tim;
inline void init(const int&up){lim=1,tim=0;while(lim<=up)lim<<=1,++tim;}
inline void ntt(vector<int>&a,int type){
	for(ri i=0;i<lim;++i)if(i<pos[tim][i])swap(a[i],a[pos[tim][i]]);
	for(ri i=1,a0,a1,t=0;i<lim;i<<=1,++t)for(ri j=0,len=i<<1;j<lim;j+=len)
	for(ri k=0;k<i;++k)a0=a[j+k],a1=mul(a[j+k+i],w[t][k]),a[j+k]=add(a0,a1),a[j+k+i]=dec(a0,a1);
	if(~type)return;
	reverse(++a.begin(),a.end());
	for(ri i=0;i<lim;++i)Mul(a[i],invv[tim]);
}
struct poly{
	vector<int>a;
	poly(int k=0,int x=0){a.resize(k+1),a[k]=x;}
	inline int&operator[](const int&k){return a[k];}
	inline const int&operator[](const int&k)const{return a[k];}
	inline int deg()const{return a.size()-1;}
	inline poly extend(const int&k){poly ret=*this;return ret.a.resize(k+1),ret;}
	inline poly rev(){poly ret=*this;return reverse(ret.a.begin(),ret.a.end()),ret;}
};
inline poly operator*(poly a,poly b){
	int n=a.deg(),m=b.deg();
	init(n+m);
	a.a.resize(lim),ntt(a.a,1);
	b.a.resize(lim),ntt(b.a,1);
	for(ri i=0;i<lim;++i)Mul(a[i],b[i]);
	return ntt(a.a,-1),a;
}
inline poly operator-(poly a,poly b){
	int n=a.deg(),m=b.deg();
	poly ret(max(n,m));
	for(ri i=0;i<=n;++i)Add(ret[i],a[i]);
	for(ri i=0;i<=m;++i)Dec(ret[i],b[i]);
	return ret;
}
inline poly operator*(poly a,int b){for(ri i=0,n=a.deg();i<=n;++i)Mul(a[i],b);return a;}
inline poly poly_inv(poly a,int k){
	if(k==1)return poly(0,ksm(a[0],mod-2));
	a.extend(k);
	poly f0=poly_inv(a,k>>1);
	return f0*2-((f0*f0)*a).extend(k);
}
inline poly Poly_inv(poly a,int k){
	int up=1;
	while(up<k)up<<=1;
	return poly_inv(a.extend(up),up).extend(k);
}
inline poly ntt_dc(int l,int r){
	if(l==r){poly ret(1,mod-l);return ret[0]=1,ret;}
	int mid=l+r>>1;
	return (ntt_dc(l,mid)*ntt_dc(mid+1,r)).extend(r-l+1);
}
int n,k;
int main(){
	ntt_init();
	cin>>n>>k;
	poly t=(Poly_inv(ntt_dc(1,k).extend(max(k,n-k)),max(k,n-k)))*poly(k,1).extend(n);
	for(ri i=0;i<=n;++i)cout<<t[i]<<' ';
	return 0;
}

倍增 $n t t$ ：

#include<bits/stdc++.h>
#define ri register int
using namespace std;
typedef long long ll;
const int mod=167772161;
inline int add(const int&a,const int&b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(const int&a,const int&b){return a>=b?a-b:a-b+mod;}
inline int mul(const int&a,const int&b){return (ll)a*b%mod;}
inline void Add(int&a,const int&b){a=a+b>=mod?a+b-mod:a+b;}
inline void Dec(int&a,const int&b){a=a>=b?a-b:a-b+mod;}
inline void Mul(int&a,const int&b){a=(ll)a*b%mod;}
inline int ksm(int a,int p){int ret=1;for(;p;p>>=1,a=mul(a,a))if(p&1)Mul(ret,a);return ret;}
vector<int>w[21],pos[21];
int invv[21];
const int Lim=1<<21;
inline void ntt_init(){
	for(ri mt=1,inv=(mod+1)>>1,W,i=1,t=0;i<Lim;i<<=1,++t,Mul(mt,inv)){
		invv[t]=mt;
		pos[t].resize(i),w[t].resize(i),w[t][0]=1,W=ksm(3,(mod-1)>>(t+1));
		for(ri j=1;j<i;++j)w[t][j]=mul(w[t][j-1],W),pos[t][j]=(pos[t][j>>1]>>1)|((j&1)<<(t-1));
	}
}
int lim,tim;
const int N=1<<20|5;
int fac[N],ifac[N];
inline void fac_init(const int&up){
	fac[0]=fac[1]=ifac[0]=ifac[1]=1;
	for(ri i=2;i<=up;++i)fac[i]=mul(fac[i-1],i),ifac[i]=mul(ifac[mod-mod/i*i],mod-mod/i);
	for(ri i=2;i<=up;++i)Mul(ifac[i],ifac[i-1]);
}
inline void init(const int&up){lim=1,tim=0;while(lim<=up)lim<<=1,++tim;}
inline void ntt(vector<int>&a,int type){
	for(ri i=0;i<lim;++i)if(i<pos[tim][i])swap(a[i],a[pos[tim][i]]);
	for(ri i=1,a0,a1,t=0;i<lim;i<<=1,++t)for(ri j=0,len=i<<1;j<lim;j+=len)
	for(ri k=0;k<i;++k)a0=a[j+k],a1=mul(a[j+k+i],w[t][k]),a[j+k]=add(a0,a1),a[j+k+i]=dec(a0,a1);
	if(~type)return;
	reverse(++a.begin(),a.end());
	for(ri i=0;i<lim;++i)Mul(a[i],invv[tim]);
}
struct poly{
	vector<int>a;
	poly(int k=0,int x=0){a.resize(k+1),a[k]=x;}
	inline int&operator[](const int&k){return a[k];}
	inline const int&operator[](const int&k)const{return a[k];}
	inline int deg()const{return a.size()-1;}
	inline poly extend(const int&k){poly ret=*this;return ret.a.resize(k+1),ret;}
	inline poly rev(){poly ret=*this;return reverse(ret.a.begin(),ret.a.end()),ret;}
};
inline poly operator*(poly a,poly b){
	int n=a.deg(),m=b.deg();
	init(n+m);
	a.a.resize(lim),ntt(a.a,1);
	b.a.resize(lim),ntt(b.a,1);
	for(ri i=0;i<lim;++i)Mul(a[i],b[i]);
	ntt(a.a,-1);
	return a;
}
inline poly operator-(poly a,poly b){
	int n=a.deg(),m=b.deg();
	poly ret(max(n,m));
	for(ri i=0;i<=n;++i)Add(ret[i],a[i]);
	for(ri i=0;i<=m;++i)Dec(ret[i],b[i]);
	return ret;
}
inline poly operator*(poly a,int b){for(ri i=0,n=a.deg();i<=n;++i)Mul(a[i],b);return a;}
inline poly poly_inv(poly a,int k){
	if(k==1)return poly(0,ksm(a[0],mod-2));
	a.extend(k);
	poly f0=poly_inv(a,k>>1);
	return (f0*(poly(0,2)-(f0*a).extend(k))).extend(k);
}
inline poly Poly_inv(poly a,int k){
	int up=1;
	while(up<k)up<<=1;
	return poly_inv(a,up).extend(k);
}
inline poly getffp(int k){
	if(k==1)return poly(1,1);
	if(k&1){
		poly a=getffp(k-1),b=a;
		a.a.push_back(0);
		for(ri i=k;i;--i)a[i]=a[i-1];
		a[0]=0;
		for(ri i=0;i<k;++i)Dec(a[i],mul(b[i],k-1));
		return a;
	}
	poly L=getffp(k>>1),t(k>>1),R(k>>1);
	for(ri i=0,n=k>>1;i<=n;++i)Mul(L[i],fac[i]);
	for(ri i=0,mt=1,n=k>>1;i<=n;++i,Mul(mt,n))t[i]=mul(mt,i&1?mod-ifac[i]:ifac[i]);
	t=L*t.rev();
	for(ri i=0,n=k>>1;i<=n;++i)R[i]=mul(t[n+i],ifac[i]),Mul(L[i],ifac[i]);
	return (L*R).extend(k);
}
inline poly ntt_dc(int l,int r){
	if(l==r){poly ret(1,1);return ret[0]=mod-l,ret;}
	int mid=l+r>>1;
	return (ntt_dc(l,mid)*ntt_dc(mid+1,r)).extend(r-l+1);
}
inline poly fix(poly a){
	int n=a.deg();
	for(ri i=0;i<n;++i)a[i]=a[i+1];
	return a.a.pop_back(),a;
}
int n,k;
int main(){
	ntt_init();
	fac_init(1<<20);
	cin>>n>>k;
	poly t=((Poly_inv(fix(getffp(k+1)).rev(),max(k,n-k))).extend(max(n-k,0)))*poly(k,1);
	for(ri i=0;i<=n;++i)cout<<t[i]<<' ';
	return 0;
}