第一类斯特林数求行

solution

首先我们有

\[x^{\overline n}=\sum_{i=0}^n\left[\begin{matrix}n\\i\end{matrix}\right]x^i \]

因此只用求出\(x^{\overline n}\) 的各项系数就可以算出一行的第一类斯特林数了

如何快速计算？考虑倍增。假设已经计算出\(x^{\overline n}\) ,要求\(x^{\overline{2n}}\) ,首先我们有

\[x^{\overline {2n}}=x^{\overline n}\cdot (x+n)^{\overline n} \]

由此问题转化为已知\(f(x)\) ,如何求出\(f(x+n)\) 。推导如下(以下记\(f_i=[x^i]f(x)\))

\[f(x+n)=\sum_{i=0}^ka_i(x+n)^i\\ =\sum_{i=0}^ka_i\sum_{j=0}^i\binom ijx^jn^{i-j}\\ =\sum_{j=0}^kx^j\sum_{i=j}^ka_i\binom ijn^{i-j}\\ =\sum_{j=0}^kx^j\sum_{i=0}^{k-j}a_{i+j}\binom {i+j}jn^{i}\\ =\sum_{j=0}^k\frac 1{j!}x^j\sum_{i=0}^{k-j}a_{i+j} {(i+j)!}\frac {n^{i}}{i!}\\ \]

令

\[B_i=a_ii!\\ A_i=\frac {n^i}{i!} \]

则

\[=\sum_{j=0}^k\frac 1{j!}x^j\sum_{i=0}^{k-j}B_{i+j}A_{i}\\ =\sum_{j=0}^k\frac 1{j!}x^j\sum_{i=0}^{k-j}Br_{k-i-j}A_{i}\\ \]

其中\(Br\) 满足\(Br_{i}=B_{k-i}\) .令\(C=Br*A\) ，则

\[=\sum_{j=0}^k\frac 1{j!}x^jC_{k-j}\\ =\sum_{j=0}^k\frac {C_{k-j}}{j!}x^j\\ \]

这样就可以求出\(f(x+n)\) 的各项系数，然后再和\(f(x)\) 卷积即可

time complexity

\[T(n)=T(\frac n2)+\mathcal O(n\log_2n) \]

根据\(master\) 定理得\(T(n)=\mathcal O(n\log_2n)\)

code

#include<bits/stdc++.h>
using namespace std;
const int N=262150,G=3,ivG=55924054,mod=167772161;
struct poly
{
	vector<int>v;
	inline int&operator[](int x)
	{
		while(v.size()<=x)v.push_back(0);
		return v[x];
	}
	inline void pre(int p,int lim)
	{
		int t=min((int)v.size(),lim);
		for(int i=p;i<t;++i)v[i]=0;
	}
};
namespace Math
{
	int inv[N],fac[N],fiv[N];
	inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
	inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
	inline int qpow(int x,int y)
	{
		int ans=1;
		for(;y;y>>=1,x=1ll*x*x%mod)
			if(y&1)ans=1ll*ans*x%mod;
		return ans;
	}
	inline void init(int n)
	{
		inv[1]=1;for(int i=2;i<=n;++i)inv[i]=mod-1ll*mod/i*inv[mod%i]%mod;
		fac[0]=1;for(int i=1;i<=n;++i)fac[i]=1ll*fac[i-1]*i%mod;
		fiv[n]=qpow(fac[n],mod-2);for(int i=n-1;~i;--i)fiv[i]=1ll*fiv[i+1]*(i+1)%mod;
	}
}
using namespace Math;
namespace Basis
{
	int rev[N<<2],fr[N<<2],wn[N<<2];
	inline void pre(int lim,int l)
	{
		for(int i=0;i<lim;++i)
			rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	}
	inline void NTT(int lim,poly&f,bool tp)
	{
		for(int i=0;i<lim;++i)fr[i]=f[rev[i]];
		for(int mid=1;mid<lim;mid<<=1)
		{
			int len=mid<<1,t=qpow(tp?ivG:G,(mod-1)/len);
			wn[0]=1;for(int i=1;i<mid;++i)wn[i]=1ll*wn[i-1]*t%mod;
			for(int j=0;j+len-1<lim;j+=len)
				for(int k=0;k<mid;++k)
				{
					int x=fr[j+k],y=1ll*wn[k]*fr[j+mid+k]%mod;
					fr[j+k]=add(x,y),fr[j+mid+k]=dec(x,y);
				}
		}
		for(int i=0;i<lim;++i)f[i]=fr[i];
	}
	inline poly mul(int n,poly f,int m,poly g)
	{
		int lim=1,l=0;while(lim<n+m)lim<<=1,++l;
		pre(lim,l);f.pre(n,lim),g.pre(m,lim);
		NTT(lim,f,0);NTT(lim,g,0);
		for(int i=0;i<lim;++i)f[i]=1ll*f[i]*g[i]%mod;
		NTT(lim,f,1);int iv=qpow(lim,mod-2);
		for(int i=0;i<lim;++i)f[i]=1ll*f[i]*iv%mod;
		return f;
	}
}
poly solve(int n)
{
	if(n==1){poly f;f[0]=1;return f;}
	if(n&1)
	{
		int m=n/2+1;
		poly f=solve(m),A,Br,fm;
		A[0]=1;for(int i=1;i<m;++i)A[i]=1ll*A[i-1]*(m-1)%mod*inv[i]%mod;
		for(int i=0;i<m;++i)Br[m-1-i]=1ll*f[i]*fac[i]%mod;
		poly C=Basis::mul(m,A,m,Br);
		for(int i=0;i<m;++i)fm[i]=1ll*C[m-1-i]*fiv[i]%mod;
		return Basis::mul(m,fm,m,f);
	}
	else
	{
		poly f=solve(n-1),g;
		for(int i=1;i<n;++i)g[i]=add(f[i-1],1ll*(n-2)*f[i]%mod);
		return g;
	}
}
int main()
{
	int n;scanf("%d",&n);init(n);
	poly f=solve(n+1);
	for(int i=0;i<=n;++i)printf("%d ",f[i]);
	return 0;
}