P5396 第二类斯特林数·列

第二类斯特林数·列

第二类斯特林数\(\begin{Bmatrix} n \\m \end{Bmatrix}\)表示把n个不同元素划分成m个相同的集合（不能有空集）的方案数。

给定\(n,k\)，对于所有的整数\(i\in[0,n]\)，你要求出\(\begin{Bmatrix} i \\k \end{Bmatrix}\)

由于答案会非常大，所以你的输出需要对\(167772161(2^{25}\cdot 5+1)\)，是一个质数）取模。

\(n,k\le 131072=2^{17}\)

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

设\(F_k=\sum a_ix^i=\sum s(i,k)x^i\)

\(F_k=\sum (s(i-1,k-1)+ks(i-1,k))x^i=xF_{k-1}+xkF_k=\frac{x}{1-kx}F_{k-1}\)

显然\(F_0=1\)故\(F_k=x^k\Pi_{i=1}^{k}\frac{1}{1-ix}\)

而\(\Pi_{i=1}^k(1-ix)=x^k\Pi_{i=1}^k(x^{-1}-i)\)

设\(y=x^{-1}\)则\(RHS=y^{-(k+1)}y^{\underline {k+1}}\)

考虑求\(y^{\underline {n}}\)

若\(n\&1==1\)，\(y^{\underline {n}}=y^{\underline {n-1}}(y-n+1)\)

若\(n\&1==0\)，设\(m=\frac{n}{2}\)，\(y^{\underline {n}}=y^{\underline {2m}}=y^{\underline {m}}(y-m)^{\underline {m}}\)

设\(g(y)=y^{\underline {m}}=\sum_{i=0}^ma_iy^i,g(y-m)=\sum_{i=0}^ma_i(y-m)^i=\sum_{i=0}^ma_i\sum_{j=0}^iC(i,j)y^j(-m)^{i-j}=\sum_{j=0}^m\sum_{i=j}^mC(i,j)a_iy^j(-m)^{i-j}=\sum_{j=0}^m\frac{y^j}{j!}\sum_{i=j}^ma_ii!\frac{(-m)^{i-j}}{(i-j)!}\)

卷积即可这部分可以复杂度为\(nlogn\)

之后将序列翻转即可得到\(\Pi_{i=1}^k(1-ix)\)

再利用多项式求逆可得\(\Pi_{i=1}^{k}\frac{1}{1-ix}\)

就可以求出\(F_k\)。

code

//#include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<ctime>
#include<cctype>
#include<queue>
#include<deque>
#include<stack>
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cstring>
#include<string>
#include<ctime>
#include<cmath>
#include<cctype>
#include<cstdlib>
#include<queue>
#include<deque>
#include<stack>
#include<vector>
#include<algorithm>
#include<utility>
#include<bitset>
#include<set>
#include<map>
#define ll long long
#define db double
#define INF 1000000000
#define inf 100000000000000000ll
#define ldb long double
#define pb push_back
#define put_(x) printf("%d ",x)
#define putl_(x) printf("%lld ",x)
#define get(x) x=read()
#define putl(x) printf("%lld\n",x)
#define rep(p,n,i) for(int i=p;i<=n;i+=1)
#define fep(n,p,i) for(int i=n;i>=p;--i)
#define go(x) for(int i=lin[x],tn=ver[i];i;tn=ver[i=nex[i]])
#define pii pair<int,int>
#define mk make_pair
#define gf(x) scanf("%lf",&x)
#define pf(x) ((x)*(x))
#define uint unsigned long long
#define ui unsigned
#define sq sqrt
#define x(w) t[w].x
#define r(w) t[w].r
#define id(w) t[w].id
#define R(w) s[w].r
#define yy p<<1|1
#define zz p<<1
#define sum(w) t[w].sum
#define mod 167772161
#define sc(A) scanf("%d",&A)
#define scs(A) scanf("%s",A);
#define put(A) printf("%d\n",A)
#define min(x,y) (x>=y?y:x)
#define max(x,y) (x>=y?x:y)
#define sub(x,y) (x-y<0?x-y+mod:x-y)
using namespace std;
const int MAXN=1<<18|1,G=3;
int n,k,m,T,cnt,lim;
int a[MAXN],b[MAXN],f[MAXN],g[MAXN],fac[MAXN],inv[MAXN],rev[MAXN];
inline int ksm(int b,int p)
{
	int cnt=1;
	while(p)
	{
		if(p&1)cnt=(ll)cnt*b%mod;
		b=(ll)b*b%mod;
		p=p>>1;
	}
	return cnt;
}
inline void NTT(int *a,int op)
{
	rep(0,lim-1,i)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int len=2;len<=lim;len=len<<1)
	{
		int mid=len>>1;
		int wn=ksm(G,op==1?(mod-1)/len:mod-1-(mod-1)/len);
		for(int j=0;j<lim;j+=len)
		{
			int d=1;
			for(int i=0;i<mid;++i)
			{
				int x=a[i+j],y=(ll)a[i+j+mid]*d%mod;
				a[i+j]=(x+y)%mod;a[i+j+mid]=(x-y+mod)%mod;
				d=(ll)d*wn%mod;
			}
		}
	}
	if(op==-1)
	{
		int IN=ksm(lim,mod-2);
		rep(0,lim-1,i)a[i]=(ll)a[i]*IN%mod;
	}
}
inline void CDQ(int n)
{
	if(n==1)
	{
		f[1]=1;
		return;
	}
	if(n&1)
	{
		CDQ(n-1);
		fep(n+1,1,i)f[i]=(f[i-1]+(ll)f[i]*(mod-n+1))%mod;
		f[0]=(ll)f[0]*(mod-n+1)%mod;
	}
	else
	{
		CDQ(n>>1);
		int m=n>>1;
		int ww=1;
		fep(m,0,i)
		{
			b[i]=(ll)ww*inv[m-i]%mod;
			a[i]=(ll)f[i]*fac[i]%mod;
			ww=(ll)ww*(mod-m)%mod;
		}
		lim=1;
		while(lim<=m+m)lim=lim<<1;
		rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
		NTT(a,1);NTT(b,1);
		rep(0,lim-1,i)a[i]=(ll)a[i]*b[i]%mod;
		NTT(a,-1);
		rep(m,m+m,j)g[j-m]=(ll)a[j]*inv[j-m]%mod;
		NTT(f,1);NTT(g,1);
		rep(0,lim-1,i)f[i]=(ll)f[i]*g[i]%mod;
		NTT(f,-1);
		rep(0,lim-1,i)
		{
			a[i]=b[i]=g[i]=0;
			if(i>n)f[i]=0;
		}
	}
}
inline void sol(int n)
{
	if(n==1)
	{
		g[0]=ksm(f[0],mod-2);
		return;
	}
	sol((n+1)>>1);
	rep(0,n-1,i)a[i]=g[i],b[i]=f[i];
	lim=1;
	while(lim<=n-1+n-1)lim=lim<<1;
	rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
	NTT(a,1);NTT(b,1);
	rep(0,lim-1,i)a[i]=(ll)a[i]*b[i]%mod*a[i]%mod;
	NTT(a,-1);
	rep(0,lim-1,i)
	{
		if(i<=n-1)g[i]=((ll)2*g[i]-a[i]+mod)%mod;
		a[i]=0;b[i]=0;
	}
}
signed main()
{
	freopen("1.in","r",stdin);
	sc(n);sc(k); 
	//put(ksm(3,(mod-1)/2));
	//put(ksm(3,(mod-1)/5));
	//由此可得mod的原根为3
	fac[0]=1;
	rep(1,k,i)fac[i]=(ll)fac[i-1]*i%mod;
	inv[k]=ksm(fac[k],mod-2);
	for(int i=k-1;i>=0;--i)inv[i]=(ll)inv[i+1]*(i+1)%mod;
	CDQ(k+1);
	rep(0,k,i)f[i]=f[i+1];
	f[k+1]=0;
	rep(0,k/2,i)swap(f[i],f[k-i]);
	if(k<=n)sol(n+1-k);
	rep(0,n,i)
	if(i>=k)
		put_(g[i-k]);
	else
		put_(0);
	return 0;
}