LG P5408 第一类斯特林数·行

Description

第一类斯特林数$egin{bmatrix}n\ mend{bmatrix}$表示将$n$个**不同**元素构成$m$个圆排列的数目。

给定$n$，对于所有的整数$iin[0,n]$，你要求出$egin{bmatrix}n\ iend{bmatrix}$。

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

Solution

因为第一类斯特林数的生成函数

$$x^{overline{n}}=sum_{i=0}^negin{bmatrix}n\iend{bmatrix}x^i$$

问题转化为求上升幂

如果使用分治+FFT可以做到$O(nlog^2n)$

使用倍增，如果要求$x^{overline{2n}}$，可以先求出$x^{overline{n}}$，再考虑算出$(x+n)^{overline{n}}$

推式子可得

$$(x+n)^{overline{n}}=sum_{i=0}^n f_i(x+n)^i=sum_{j=0}^n frac{x^j}{j!} sum _{i=j}^n i!f_ifrac{n^{i-j}}{(i-j)!} $$

所以可以FFT求解$(x+n)^{overline{n}}$

当$n$为奇数时可以特判

总复杂度$O(nlog n)$

#include<iostream>
#include<cstdio>
using namespace std;
int N,rev[1050005],s,tot;
long long fac[1050005]={1},inv[1050005],f[1050005],A[1050005],B[1050005];
const long long mod=167772161;
inline int read()
{
    int w=0,f=1;
    char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') f=-1; ch=getchar();}
    while(ch>='0'&&ch<='9')w=(w<<1)+(w<<3)+ch-'0',ch=getchar();
    return w*f;
}
long long ksm(long long a,long long p)
{
    long long ret=1;
    while(p)
    {
        if(p&1) (ret*=a)%=mod;
        (a*=a)%=mod,p>>=1;
    }
    return ret;
}
void ntt(long long *a,int n,int INV)
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i<<=1)
    {
        long long wn=ksm(3,(mod-1)/i/2);
        if(INV==-1) wn=ksm(wn,mod-2);
        for(int j=0;j<n;j+=i*2)
        {
            long long w=1;
            for(int k=j;k<i+j;k++)
            {
                long long x=a[k],y=w*a[k+i]%mod;
                a[k]=(x+y)%mod,a[k+i]=(x-y+mod)%mod,(w*=wn)%=mod;
            }
        }
    }
    if(INV==-1)
    {
        long long temp=ksm(n,mod-2);
        for(int i=0;i<n;i++) (a[i]*=temp)%=mod;
    }
}
void trans(long long *f,int n,long long *g)
{
    long long base=1;
    s=2,tot=1;
    while(s<=n*2) s<<=1,tot++;
    for(int i=0;i<=n;i++) A[n-i]=f[i]*fac[i]%mod;
    for(int i=0;i<=n;i++,(base*=n)%=mod) B[i]=base*inv[i]%mod;
    for(int i=n+1;i<s;i++) A[i]=B[i]=0;
    for(int i=0;i<s;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(tot-1));
    ntt(A,s,1),ntt(B,s,1);
    for(int i=0;i<s;i++) (A[i]*=B[i])%=mod;
    ntt(A,s,-1);
    for(int i=0;i<=n;i++) g[i]=A[n-i]*inv[i]%mod;
}
void solve(int n,long long *f)
{
    if(!n) return void(f[0]=1);
    int m=n/2;
    solve(m,f);
    trans(f,m,B),s=2,tot=1;
    while(s<=n) s<<=1,++tot;
    for(int i=0;i<=m;i++) A[i]=f[i];
    for(int i=m+1;i<s;i++) A[i]=B[i]=0;
    for(int i=0;i<s;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(tot-1));
    ntt(A,s,1),ntt(B,s,1);
    for(int i=0;i<s;i++) (A[i]*=B[i])%=mod;
    ntt(A,s,-1);
    if(n&1) for(int i=0;i<=n;i++) f[i]=((i?A[i-1]:0)+(n-1)*A[i])%mod;
    else for(int i=0;i<=n;i++) f[i]=A[i];
}
int main()
{
    for(int i=1;i<=1050000;i++) fac[i]=fac[i-1]*i%mod;
    inv[1050000]=ksm(fac[1050000],mod-2);
    for(int i=1049999;~i;i--) inv[i]=inv[i+1]*(i+1)%mod;
    N=read(),solve(N,f);
    for(int i=0;i<=N;i++) printf("%lld ",f[i]);
    return 0;
}