第一類斯特林數·行

題目

https://www.luogu.com.cn/problem/P5408

求

x n ‾ = ∑ i = 0 n [ n i ] x i x^{\overline{n}}=\sum_{i=0}^{n}\begin{bmatrix}n\\i\end{bmatrix}x^i xn=i=0∑n​[ni​]xi

思路

x 2 k ‾ = x k ‾ ( x + k ) k ‾ x^{\overline{2k}}=x^{\overline{k}}(x+k)^{\overline{k}} x2k=xk(x+k)k。假設知道 x k ‾ x^{\overline{k}} xk就能快速求出 ( x + k ) k ‾ (x+k)^{\overline{k}} (x+k)k，那麼就能用類似快速幂的方式求出。

現在考慮一個問題，假設知道一個多項式 f ( x ) = ∑ i = 0 n a i x i f(x)=\sum_{i=0}^{n}a_ix^i f(x)=∑i=0n​ai​xi，如何求

f ( x + k ) = ∑ i = 0 n b i x i f(x+k)=\sum_{i=0}^{n}b_ix^i f(x+k)=∑i=0n​bi​xi

f ( x + k ) = ∑ i = 0 n a i ( x + k ) i = ∑ i = 0 n ∑ j = 0 i ( i j ) a i x j k i − j = ∑ j = 0 n x j 1 j ! ∑ i = j n i ! a i ( i − j ) ! k i − j = ∑ j = 0 n x j 1 j ! ∑ i = 0 n − j k i i ! ( i + j ) ! a i + j = ∑ i = 0 n x i 1 i ! ∑ j = 0 n − i k j j ! ( i + j ) ! a i + j \begin{aligned} f(x+k)&=\sum_{i=0}^{n}a_i(x+k)^i\\ &=\sum_{i=0}^{n}\sum_{j=0}^{i}{i\choose j}a_ix^jk^{i-j}\\ &=\sum_{j=0}^{n}x^j\frac{1}{j!}\sum_{i=j}^{n}\frac{i!a_i}{(i-j)!}k^{i-j}\\ &=\sum_{j=0}^{n}x^j\frac{1}{j!}\sum_{i=0}^{n-j}\frac{k^{i}}{i!}(i+j)!a_{i+j}\\ &=\sum_{i=0}^{n}x^i\frac{1}{i!}\sum_{j=0}^{n-i}\frac{k^{j}}{j!}(i+j)!a_{i+j}\\ \end{aligned} f(x+k)​=i=0∑n​ai​(x+k)i=i=0∑n​j=0∑i​(ji​)ai​xjki−j=j=0∑n​xjj!1​i=j∑n​(i−j)!i!ai​​ki−j=j=0∑n​xjj!1​i=0∑n−j​i!ki​(i+j)!ai+j​=i=0∑n​xii!1​j=0∑n−i​j!kj​(i+j)!ai+j​​

是以

b i = 1 i ! ∑ j = 0 n − i k j j ! ( i + j ) ! a i + j b_i=\frac{1}{i!}\sum_{j=0}^{n-i}\frac{k^{j}}{j!}(i+j)!a_{i+j} bi​=i!1​j=0∑n−i​j!kj​(i+j)!ai+j​

令 c i = k i i ! , d i = ( n − i ) ! a n − i c_i=\frac{k^i}{i!},d_i=(n-i)!a_{n-i} ci​=i!ki​,di​=(n−i)!an−i​，則

b i = 1 i ! ∑ j = 0 n − i c j d n − i − j b_i=\frac{1}{i!}\sum_{j=0}^{n-i}c_jd_{n-i-j} bi​=i!1​j=0∑n−i​cj​dn−i−j​

這樣就有卷積的形式，再把 b i b_i bi​映射到 b n − i ′ b^{'}_{n-i} bn−i′​，則

b i = b n − i ′ = 1 i ! ∑ j = 0 n − i c j d n − i − j b i ′ = 1 ( n − i ) ! ∑ j = 0 i c j d i − j \begin{aligned} b_i=b^{'}_{n-i}&=\frac{1}{i!}\sum_{j=0}^{n-i}c_jd_{n-i-j}\\ b^{'}_i&=\frac{1}{(n-i)!}\sum_{j=0}^{i}c_jd_{i-j}\\ \end{aligned} bi​=bn−i′​bi′​​=i!1​j=0∑n−i​cj​dn−i−j​=(n−i)!1​j=0∑i​cj​di−j​​

這樣用 N T T NTT NTT求出了 b ′ b^{'} b′，在映射回去就得到 b b b，這樣就解決了問題，現在 f ( x ) = x k ‾ f(x)=x^{\overline{k}} f(x)=xk，求出 f ( x + k ) = ( x + k ) k ‾ f(x+k)=(x+k)^{\overline{k}} f(x+k)=(x+k)k

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const double Pi=acos(-1.0);
const int N=610009;
const ll mod=167772161,G=3;//G是mod的原根
ll a[N],b[N],c[N],res[N],p[N],p1[N],Gi,_inv,jie[N],jie_inv[N],f_x[N];//Gi是原根的逆元,inv是lim的逆元
int n,m,bit,lim,r[N];//lim表示目前運算的長度
ll qpow(ll a,ll b) {
    ll res=1;
    a%=mod;
    while(b) {
        if(b&1)
            res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}
void _init(int n) { //多項式可能最長的長度,隻需初始化一次
    Gi=qpow(G,mod-2);
    for(int i=1; i<n; i<<=1)
        p[i]=qpow(G,(mod-1)/(i<<1)),p1[i]=qpow(Gi,(mod-1)/(i<<1));
}
void init(int n,int m) { //n,m表示最高長度,每次運算調用一次
    lim=1,bit=0;
    while(lim<n+m-1)
        lim<<=1,bit++;
    _inv=qpow(lim,mod-2);
    for(int i=0; i<lim; i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
}
void NTT(ll *a,int type) {
    for(int i=0; i<lim; i++)
        if(i<r[i])
            swap(a[i],a[r[i]]);
    for(int mid=1; mid<lim; mid<<=1) {
        ll W=(type==1?p[mid]:p1[mid]);
        for(int r=mid<<1,j=0; j<lim; j+=r) {
            ll e=1;
            for(int k=0; k<mid; k++,e=e*W%mod) {
                ll x=a[j+k],y=e*a[j+k+mid]%mod;
                a[j+k]=(x+y)%mod,a[j+k+mid]=(x-y+mod)%mod;
            }
        }
    }
    if(type==-1)
        for(int i=0; i<lim; i++)
            a[i]=a[i]*_inv%mod;
}
void NTT(ll a[],int n,ll b[],int m,ll *c) { //n,m表示最高長度,c存答案,但是c可以和a是同一個
    init(n,m);
    NTT(a,1),NTT(b,1);
    for(int i=0; i<lim; i++)
        c[i]=a[i]*b[i]%mod;
    NTT(a,-1),NTT(b,-1);
    if((*a)!=(*c))
        NTT(c,-1);
}
void sol(ll *a,int n,ll k) {//知道f(x),求f(x+k),n為a的長度,答案存f_x,f_x用完清零
    ll tmp=1;
    for(int i=0; i<=n; i++)
        b[i]=tmp*jie_inv[i]%mod,tmp=tmp*k%mod;
    for(int i=0; i<=n; i++)
        c[i]=jie[n-i]*a[n-i]%mod;
    init(n+1,n+1);
    NTT(b,1),NTT(c,1);
    for(int i=0; i<lim; i++)
        b[i]=b[i]*c[i]%mod;
    NTT(b,-1);
    for(int i=0; i<=n; i++)
        b[i]=b[i]*jie_inv[n-i]%mod;
    for(int i=0; i<=n; i++)
        f_x[i]=b[n-i];
    for(int i=0; i<lim; i++) //清零
        b[i]=c[i]=0;
}
void cal(int ex) {
    if(ex==1) {
        res[1]=1;
        return ;
    }
    if(ex&1) {
        cal(ex-1);
        for(int i=ex; i>=1; i--)
            res[i]=(res[i]*(ex-1)%mod+res[i-1])%mod;
        res[0]=res[0]*(ex-1)%mod;
    } else {
        cal(ex/2);
        sol(res,ex/2+1,ex/2);
        NTT(res,ex/2+1,f_x,ex/2+1,res);
        for(int i=ex+1; i<lim; i++)
            res[i]=0;
        for(int i=0; i<lim; i++)
            f_x[i]=0;
    }
}
int main() {
    _init(N);
    jie[0]=1;
    for(int i=1; i<N; i++)
        jie[i]=jie[i-1]*i%mod;
    jie_inv[N-1]=qpow(jie[N-1],mod-2);
    for(int i=N-2; i>=0; i--)
        jie_inv[i]=jie_inv[i+1]*(i+1)%mod;
    scanf("%d",&n);
    cal(n);
    for(int i=0; i<=n; i++)
        printf("%lld ",res[i]);
    return 0;
}