solution
首先我们有
[x^{overline n}=sum_{i=0}^nleft[egin{matrix}n\iend{matrix}
ight]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_isum_{j=0}^iinom ijx^jn^{i-j}\
=sum_{j=0}^kx^jsum_{i=j}^ka_iinom ijn^{i-j}\
=sum_{j=0}^kx^jsum_{i=0}^{k-j}a_{i+j}inom {i+j}jn^{i}\
=sum_{j=0}^kfrac 1{j!}x^jsum_{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}^kfrac 1{j!}x^jsum_{i=0}^{k-j}B_{i+j}A_{i}\
=sum_{j=0}^kfrac 1{j!}x^jsum_{i=0}^{k-j}Br_{k-i-j}A_{i}\
]
其中(Br) 满足(Br_{i}=B_{k-i}) .令(C=Br*A) ，则
[=sum_{j=0}^kfrac 1{j!}x^jC_{k-j}\
=sum_{j=0}^kfrac {C_{k-j}}{j!}x^j\
]
这样就可以求出(f(x+n)) 的各项系数，然后再和(f(x)) 卷积即可
time complexity
[T(n)=T(frac n2)+mathcal O(nlog_2n)
]
根据(master) 定理得(T(n)=mathcal O(nlog_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;
}