多项式

· · 算法·理论

多项式全家桶

多项式

定义

形如 \sum a_nx^n 的求和式若为有限项相加，那么称作多项式 记作f(x)=\sum_{n=0}^m a_nx^n

对于如 \sum_{n=0}^{\infty}a_nx^n 有一个非负整数次幂乘以一个常数，称作幂级数

运算

加减运算

两个多项式 F(x)=\sum_{n\geq0}a_nx^n 与 G(x)=\sum_{n\geq0}b_nx^n 有

F(x)\pm G(x)=\sum_{n\geq0}(a_n\pm b_n)x^n
卷积（乘法）

两个多项式 F(x)=\sum_{n\geq0}a_nx^n 与 G(x)=\sum_{n\geq0}b_nx^n 有

F(x)\times G(x)=\sum_{i\geq0}a_ix^i \sum_{j\geq0}b_jx^j\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\sum_{n\geq0}x^n\sum_{i\geq0}a_ib_{n-i}

FFT

单位根

对于 n 次单位根为 \omega_n 表示在复平面上以原点为圆心的半径为1的园，并从实部正半轴开始逆时针将园 n 等分后的原点与第一个 n 等分点构成的复数。

在弧度制下，该角弧度为 \frac{2\pi}{n} ，那么

\omega_n^k=\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n})

那么可以推导出：

(\omega_n^n)=(\omega_n^0)=1\\ \omega_{2n}^{2k}=\cos(2k\times\frac{2\pi}{2n})+i\sin(2k\times\frac{2\pi}{2n})=\omega_n^k\\ \omega_n^{\frac{n}{2}}=\cos(\frac{n}{2}\times\frac{2\pi}{n})+i\sin(\frac{n}{2}\times\frac{2\pi}{n})=\cos(\pi)+i\sin(\pi)=-1\\ \omega_n^{k+\frac n2}可以理解为在选择\frac{2k\pi}n 的基础上继续旋转\frac{2\pi \frac{n}2}{n}=\pi,即多旋转180^\circ,\\从向量的角度看，可以感性理解为其相乘，则\omega_n^{k+\frac n2}=-\omega_n^{k} (\omega_n^k)^2=(\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}))\times(\cos (k\times \frac{2\pi}{n}) +i\sin (k\times \frac{2\pi}{n}))\\ =\cos^2(k\times \frac{2\pi}{n})-\sin^2(k\times \frac{2\pi}{n}) +i2\cos(k\times \frac{2\pi}{n})\sin(k\times \frac{2\pi}{n})\\ \\\because \begin{bmatrix} \cos \alpha&-\sin\alpha\\ \sin \alpha&\cos\alpha\\ \end{bmatrix} \begin{bmatrix} \cos \alpha\\ \sin \alpha \end{bmatrix}= \begin{bmatrix} \cos^2\alpha-\sin^2\alpha\\ 2\cos\alpha\sin\alpha\\ \end{bmatrix}\\ \therefore \cos 2\alpha=\cos^2\alpha-\sin^2\alpha\\ \sin 2\alpha=2\cos\alpha\sin\alpha\\ \therefore(\omega_n^k)^2=\cos(2k\times\frac{2\pi}{n})+i\sin(2k\times\frac{2\pi}{n})=\omega_n^{2k} 可以感性理解一下(\omega_n^i)(\omega_n^j)=(\omega_n^{i+j})，因为其可以看作向量旋转

DFT

考虑将这 n 个 n 次单位根带入一个 n 次多项式，获得其点值表示法，复杂度肯定是 O(n^2) ，但是如果对下标进行奇偶性分离：

A(x)=\sum_{i=0}^{n-1} a_ix^i=a_0+a_1x^1+a_2x^2+\dots +a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\\ =a_0+a_2x^2+\dots+a_{n-2}x^{n-2} \\+a_1x^1+a_3x^3+\dots+a_{n-1}x^{n-1}\\ 令A_1(x)=a_0+a_2x^1+\dots+a_{n-2}x^{\frac n2-1}\\A_2(x)=a_1+a_3x^1+\dots+a_{n-1}x^{\frac n2-1}\\ 那么A(x)=A_1(x^2)+xA_2(x^2)\\ 将其中一个单位根 \omega_n^k \ (k<\frac n2) 带入：\\ A(\omega_n^k)=A_1((\omega_n^k)^2)+\omega_n^kA_2((\omega_n^k)^2)\\ =A_1(\omega_n^2k)+\omega_n^kA_2(\omega_n^2k)\\ 若将\omega_n^{k+\frac n2}带入：\\ A`(-\omega_n^k)=A_1(\omega_n^2k)-\omega_n^kA_2(\omega_n^2k)\\

我们可以发现，A 与 A` 只有一个常数项不同，所以在求 A 时可以 O(1) 算出第二个，所以 k 的取值范围时 [0,\frac n 2-1] ，然后发现 A_1,A_2,A 的性质相同，所以可以递归转化，复杂度为 O(n\log_2 n)。

IDFT

现在已知 (y_0,y_1,y_2,\dots,y_{n-1}) 是 (a_0,a_1,a_2,\dots,a_{n-1}) 的点值表示法。

假设有个向量 (c_0,c_1,c_2,\dots,c_{n-1}) ，即多项式 B(x)=\sum\limits_{i=0}^{n-1}y_ix^i 在 (\omega_n^0,\omega_n^{-1},\omega_n^{-2},\dots,\omega_n^{-(n-1)}) 上的值。

那么：

\begin{aligned} c_k&=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j)(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i)\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omega_n^{-k})^i\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i\\ &=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ \end{aligned} 设S(x,n)=\sum_{i=0}^{n-1}x^i\\ 当x=\omega_n^k时：\\ S(\omega_n^k,n)=1+(\omega_n^k)+(\omega_n^k)^2+\dots+(\omega_n^k)^{n-1}\\ (\omega_n^k)S(\omega_n^k,n)=(\omega_n^k)+(\omega_n^k)^2+(\omega_n^k)^3+\dots+(\omega_n^k)^n\\ 两式相减：\\ (\omega_n^k-1)S(\omega_n^k,n)=(\omega_n^k)^n-1\\ S(\omega_n^k,n)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}\\ \therefore S(\omega_n^k,n)=\frac{1-1}{\omega_n^k-1}\\ 但是在k=0时S(\omega_n^k,n)=n\\ 所以在原式c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j-k})^i\\ 当i\neq k时，后一个\sum中值为0\\ 当i= k时，后一个\sum中值为n\\ 故c_k=na_k\\ a_k=\frac{c_k}{n}

总结

那么两个多项式相乘可以先将其分别转化为点值表示法，然后对于位置相乘得到结果的点值表示法，然后在跑一遍那个将多项式转化为点值的函数，将点值转化为多项式，最后依据 IDFT 除以即可。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const double PI=acos(-1.0);

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

struct comple{
    double a,b;
    comple(double a=0,double b=0):
        a(a),b(b){;}
    comple operator +(const comple &x)const{return comple(a+x.a,b+x.b);}
    comple operator -(const comple &x)const{return comple(a-x.a,b-x.b);}
    comple operator *(const comple &x)const{return comple(a*x.a-b*x.b,a*x.b+b*x.a);}
}; 

void Fast(int now,comple *a,int tag){
    if(now==1)return;
    comple a1[(now>>1)+10],a2[(now>>1)+10];
    for(int i=0;i<(now>>1);++i){
        a1[i]=a[i<<1];
        a2[i]=a[i<<1|1];
    }
    Fast(now>>1,a1,tag);
    Fast(now>>1,a2,tag);
    comple W=comple(cos(2.0*PI/now),tag*(sin(2.0*PI/now))),w=comple(1,0);
    for(int i=0;i<(now>>1);++i,w=w*W){
        comple z=w*a2[i];
        a[i]=a1[i]+z;
        a[i+(now>>1)]=a1[i]-z;
    }
}

comple a[Max],b[Max];

bool Med;
signed main(){
    int n=read();
    int m=read();
    for(int i=0;i<=n;++i)a[i].a=raed();
    for(int i=0;i<=m;++i)b[i].a=read();
    int pos=1;
    while(pos<=n+m)pos<<=1;
    Fast(pos,a,1);
    Fast(pos,b,1);
    for(int i=0;i<=pos;++i)a[i]=a[i]*b[i];
    Fast(pos,a,-1);
    for(int i=0;i<=n+m;++i)cout << (int)(a[i].a/pos+0.5)<< ' ';

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}

实际上对于上述算法还可以进行优化，就是不递归的写法。

因为我们发现一个数的位置是其元下标二进制翻转后的位置，所有可以预处理然后枚举长度合并。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const double PI=acos(-1.0);

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

struct comple {

    double a,b;
    comple (double a=0,double b=0):
        a(a),b(b){;}
    comple operator +(const comple &x)const{return comple(a+x.a,b+x.b);}
    comple operator -(const comple &x)const{return comple(a-x.a,b-x.b);}
    comple operator *(const comple &x)const{return comple(a*x.a-b*x.b,a*x.b+b*x.a);}
};

int rev[Max];

int limit;

void Fast(comple *a,int tag){
    for(int i=0;i<limit;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        comple W=comple(cos(PI/mid),tag*sin(PI/mid));
        for(int R=mid<<1,j=0;j<limit;j+=R){
            comple w=comple(1,0);
            for(int l=0;l<mid;++l,w=w*W){
                comple x=a[l+j],y=a[l+j+mid]*w;
                a[l+j]=x+y;
                a[l+j+mid]=x-y;
            }
        }
    }
}

comple a[Max],b[Max];

bool Med;
signed main(){
    int n=read();
    int m=read();
    for(int i=0;i<=n;++i)a[i].a=read();
    for(int i=0;i<=m;++i)b[i].a=read();
    limit=1;
    int len=0;
    while(limit<=n+m)limit<<=1,++len;
    for(int i=0;i<=limit;++i)rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
    Fast(a,1);
    Fast(b,1);
    for(int i=0;i<=limit;++i)a[i]=a[i]*b[i];
    Fast(a,-1);
    for(int i=0;i<=n+m;++i)cout << (int)(a[i].a/limit+0.5) << ' ';
    cout << "\n";

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}

NTT

原根

将最小的满足 a^l\equiv1\pmod p 的 l 称作 a 在模 p 意义下的阶。

如果 g 在模 p 意义下的阶是 \varphi(p)，那么称 g 是 p 的原根。

求原根可以这样做：如果 \gcd(g,p)=1，设 p_1,p_2 \dots p_k 为 p 的所有质因子，如果 g 是 p 的原根就满足对于任意 i\in[1,k],g^{\frac{\varphi(p)}{p_i}}\not\equiv1\pmod p

可以发现原根有一个非常优秀的性质就是对于在 1 到 \varphi(p) 中的任意一个幂次都是不一样的，而且是以 \varphi(p) 为长度的一个循环节，这与 \omega 的性质很像。而且当 p 为质数时，如果 \varphi(p) 有很多 2 这个因子就十分优秀，像 998244353=2^{23}\times 7\times 17+1 它的原根是 3 。

那么现在我们可以用 g^{\frac{\varphi(p)}{n}} 代替 n 次单位根就行。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=3e6+10;
const int mod=998244353;
const int inf=1e9+10;
const int g=3;
const int invg=332748118;

template <int mod>
struct modint{

    int val;

    static int norm(const int &x){return x<0?x+mod:x;}

    modint inv()const{
        int a=val,b=mod,u=1,v=0,t;
        while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
        return modint(u);
    }
    modint():val(0){}
    modint(const int &m):val(norm(m)){}
    modint(const ll &m):val(norm(m%mod)){}
    modint operator -()const{return modint(norm(-val));}
    bool operator ==(const modint &x){return val==x.val;}
    bool operator !=(const modint &x){return val!=x.val;}
    bool operator <=(const modint &x){return val<=x.val;}
    bool operator >=(const modint &x){return val>=x.val;}
    bool operator >(const modint &x){return val>x.val;}
    bool operator <(const modint &x){return val<x.val;}
    modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
    modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
    modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
    modint& operator >>=(const modint &x){return val>>=x.val,*this;}
    modint& operator <<=(const modint &x){return val<<=x.val,*this;}
    modint& operator ^=(const modint &x){return val^=x.val,*this;}
    modint operator >>(const modint &x){return modint(*this)>>=x;}
    modint operator <<(const modint &x){return modint(*this)<<=x;}
    modint& operator /=(const modint &x){return *this*=x.inv();}
    modint operator +(const modint &x){return modint(*this)+=x;}
    modint operator -(const modint &x){return modint(*this)-=x;}
    modint operator *(const modint &x){return modint(*this)*=x;}
    modint operator /(const modint &x){return modint(*this)/=x;}
    modint operator ^(const modint &x){return modint(*this)^=x;}
    friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
    friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}

};

typedef modint<998244353>m98;

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

m98 ksm(m98 a,int b){
    m98 ans=1;
    for(;b;b>>=1){
        if(b&1)ans=ans*a;
        a=a*a;
    }
    return ans;
}

m98 a[Max],b[Max];

int rev[Max],limit;

void Fast(m98 *a,int tag){
    for(int i=0;i<limit;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
        for(int R=mid<<1,j=0;j<limit;j+=R){
            m98 w=1;
            for(int p=0;p<mid;++p,w=w*W){
                m98 x=a[p+j],y=w*a[p+j+mid];
                a[p+j]=x+y;
                a[p+j+mid]=x-y;
            }
        }
    }
    if(tag==-1){
        m98 inv=m98(limit).inv();
        for(int i=0;i<limit;++i)a[i]=a[i]*inv;
    }
}

bool Med;
signed main(){
    int n=read();
    int m=raed();
    for(int i=0;i<=n;++i)a[i].val=read();
    for(int i=0;i<=m;++i)b[i].val=raed();
    limit=1;
    int len=0;
    while(limit<=n+m)limit<<=1,++len;
    for(int i=0;i<=limit;++i)rev[i]=((rev[i>>1]>>1)|((i&1)<<(len-1)));
    Fast(a,1);
    Fast(b,1);
    for(int i=0;i<=limit;++i)a[i]=a[i]*b[i];
    Fast(a,-1);
    for(int i=0;i<=n+m;++i)cout << a[i] << ' ';
    cout << "\n";

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}
/*

*/

多项式牛顿迭代

给定多项式 g(x) 求在模 x^n 下的 f(x) 使其满足：

g(f(x))\equiv0\pmod {x^n}

通过倍增和泰勒展开(不会)得到：

f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n}

考虑证明：

假设我们已经知道 f(x) 在 (\mod x^{\lceil\frac n2\rceil}) 下的函数 f_0(x) 先要求其在 \mod x^n 发函数：

考虑将 g(f(x)) (\mod x^n) 在 f_0(x) 位置下泰勒展开：

\sum_{i\geq 0} \frac{g^{(i)}(f_0(x))}{i!}(f(x)-f_0(x))^i \equiv0\pmod {x^n}

可以注意到 f(x)-f_0(x) 的最低非零项次最低是 \lceil \frac n2\rceil 所以当 i\geq 2 时 f(x)-f_0(x)\equiv 0\pmod{x^n} 。

所以有

g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\equiv 0\pmod{x^n}\\ \therefore f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}

多项式求逆

设给定的多项式是 h(x)，那么有：

g(f(x))=\frac{1}{f(x)}-h(x)\equiv 0\pmod{x^n}

于是有：

f(x)\equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}}\pmod{x^n}\\ \equiv 2f_0(x)-f_0^2(x)h(x)\pmod{x^n}\\

而当 n=1 时，[x^0]f_0(x)=[x^0]h(x)^{-1}

多项式开根

设给定的多项式是 h(x)，那么有：

g(f(x))\equiv f^2(x)-h(x)\pmod{x^n}\\

于是有：

f(x)\equiv f_0(x)-\frac{f^2_0(x)-h(x)}{2f_0(x)}\pmod{x^n}\\ \equiv\frac{f_0^2(x)+h(x)}{2f_0(x)}

多项式对数函数

假设我们要求 G(x)\equiv\ln(F(x))\pmod{x^n} ，设函数 f(x)=\ln(x)，那么我们要求的是：

G(x)\equiv f(F(x))\pmod{x^n}\\ G'(x)=f(F(x))'\\ G'(x)=f'(F(x))F'(x)\\

注意到 \ln(x) 的导数是 \frac{1}{x}，所以有：

G'\equiv \frac{F'(x)}{F(x)} \pmod{x^n}

因为我们只关心余数，所以可以求逆就行。

​

多项式指数函数

设给定的多项式是 h(x)，那么有：

g(f(x))\equiv\ln(f(x))-h(x)\pmod{x^n}

于是有：

f(x)\equiv f_0(x)-\frac{\ln(f_0(x))-h(x)}{\frac{1}{f_0(x)}}\pmod{x^n}\\ \equiv(f_0(x))(1-\ln(f_0(x))+h(x))\pmod{x^n}

多项式快速幂

设给定的多项式是 h(x)，那么求：

G(x)\equiv h(x)^k \pmod{x^n}\\ \ln(G(x))\equiv k\ln(h(x))\pmod{x^n}\\ G(x)\equiv e^{k\ln(h(x))}\pmod{x^n}

注意到要求 \ln 所以要保证第一位位 1 ，所以对于部分不合法多项式应该左移 z 位到第一位不为 0，如果是别的数要先除以 a_0，但是最后要每一项乘以 a_0^{k\mod \varphi(p)} ，然后还要乘以一个 x^{zk\mod p}

Code

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=6e5+10;
const int mod=998244353;
const int inf=1e9+10;

template <int mod>
struct modint{

    int val;

    static int norm(const int &x){return x<0?x+mod:x;}

    modint inv()const{
        int a=val,b=mod,u=1,v=0,t;
        while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
        return modint(u);
    }

    modint():val(0){}
    modint(const int &m):val(norm(m)){}
    modint(const ll &m):val(norm(m%mod)){}
    modint operator -()const{return modint(norm(-val));}
    bool operator ==(const modint &x){return val==x.val;}
    bool operator !=(const modint &x){return val!=x.val;}
    bool operator <=(const modint &x){return val<=x.val;}
    bool operator >=(const modint &x){return val>=x.val;}
    bool operator >(const modint &x){return val>x.val;}
    bool operator <(const modint &x){return val<x.val;}
    modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
    modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
    modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
    modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
    modint& operator >>=(const modint &x){return val>>=x.val,*this;}
    modint& operator ^=(const modint &x){return val^=x.val,*this;}
    modint operator >>(const modint &x){return modint(*this)>>=x;}
    modint operator <<(const modint &x){return modint(*this)<<=x;}
    modint& operator /=(const modint &x){return *this*=x.inv();}
    modint operator +(const modint &x){return modint(*this)+=x;}
    modint operator -(const modint &x){return modint(*this)-=x;}
    modint operator *(const modint &x){return modint(*this)*=x;}
    modint operator /(const modint &x){return modint(*this)/=x;}
    modint operator ^(const modint &x){return modint(*this)^=x;}
    friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
    friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<998244353>m98;

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

m98 ksm(m98 a,int b){
    m98 ans=1;
    for(;b;b>>=1){
        if(b&1)ans=ans*a;
        a=a*a;
    }
    return ans;
}

m98 inv[Max],frac[Max];

void init(int len) {
    frac[0]=inv[0]=1; 
    for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
    inv[len]=frac[len].inv();for(int i=len-1;i;--i)inv[i]=inv[i+1]*(i+1);
    for(int i=1;i<=len;++i)inv[i]*=frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
    const int g=3;
    const int invg=332748118;

    m98 a[Max];
    void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
    void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
    void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
    void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
    m98& operator [](const int x){return a[x];}

    void NTT(int tag,int lim){
        for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            m98 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
            for(int R=mid<<1,j=0;j<lim;j+=R){
                m98 w=1;
                for(int p=0;p<mid;++p,w=w*W){
                    m98 x=a[j+p],y=a[j+p+mid]*w;
                    a[j+p]=x+y;
                    a[j+p+mid]=x-y;
                }
            }
        }
        if(tag==-1){
            m98 inv=m98(lim).inv();
            for(int i=0;i<lim;++i)a[i]*=inv;
        }
    }

}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
    static Poly b[2],c;
    int lim=2,len=1,bas=1;
    int pos=0;
    b[pos][0]=x[0].inv();
    while(bas<(n<<1)){
        init(lim,len);
        pos^=1;
        for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
        for(int i=0;i<=lim;++i)c[i]=x[i];
        for(int i=lim>>1;i<=lim;++i)c[i]=0;
        b[pos^1].NTT(1,lim);c.NTT(1,lim);
        for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
        b[pos^1].NTT(-1,lim);
        for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
        bas<<=1,lim<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=b[pos][i];
    for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetSqrt(Poly &y,Poly &x,int n){

    static Poly t1,t2,s;
    s[0]=1;
    int lim=2,len=1,bas=1;
    while(bas<(n<<1)){
        for(int i=0;i<bas;++i)t1[i]=x[i];
        for(int i=0;i<bas;++i)t2[i]=s[i]<<1;
        GetInv(t2,t2,bas);init(lim,len);
        t2.NTT(1,lim);t1.NTT(1,lim);
        for(int i=0;i<=lim;++i)t1[i]=t1[i]*t2[i];
        t1.NTT(-1,lim);
        for(int i=0;i<bas;++i)s[i]=s[i]/2+t1[i];
        bas<<=1,lim<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=s[i];
    for(int i=0;i<bas;++i)t1[i]=t2[i]=s[i]=0;

}

void GetLn(Poly &y,Poly&x,int n){
    static Poly t1,t2;
    for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
    t1.der(n);GetInv(t2,t2,n);
    int lim=1,len=0;
    while(lim<n<<1)lim<<=1,++len;
    init(lim,len);
    t1.NTT(1,lim),t2.NTT(1,lim) ;
    for(int i=0;i<lim;++i)t1[i]*=t2[i];
    t1.NTT(-1,lim);t1.inte(n) ;
    for(int i=0;i<n;++i)y[i]=t1[i];
    for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}

void GetExp(Poly &y,Poly &x,int n){
    static Poly s,tmp;s[0]=1;
    int bas=1,lim=2,len=1;
    while(bas<n<<1){
        for(int i=0;i<bas;++i)tmp[i]=s[i];
        GetLn(tmp,tmp,bas);
        for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
        tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
        for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
        s.NTT(-1,lim) ;
        for(int i=bas;i<lim;++i)s[i]=0;
        lim<<=1,bas<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=s[i];
    for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}

int Si=0;
void expow(Poly &y,Poly &x,int n,m98 k,int k1){
    static Poly s;
    int z=n+1;
    for(int i=0;i<n;++i)if(x[i]!=0){z=i;break;}
    if(z==n+1||(z!=0&&Si)){
        for(int i=0;i<n;++i)y[i]=0;
        return;
    }
    for(itn i=z;i<n;++i)s[i-z]=x[i]/x[z];GetLn(s,s,n-z);
    for(int i=0;i<n-z;++i)s[i]*=k;GetExp(s,s,n-z);
    for(int i=0;i<n-z;++i)s[i]=s[i]*ksm(x[z],k1);
    for(int i=0;i<n;++i)y[i]=0;
    for(int i=0;1ll*z*k.val+i<n;++i)y[i+z*k.val]=s[i];
    for(int i=0;i<n;++i)s[i]=0; 
}

Poly b;

m98 k;
ll k1;

void In(){
    k=k1=0;
    char c=getchar();
    while(c<'0'||c>'9')c=getchar();
    while(c>='0'&&c<='9'){Si=Si|(k*10+(c^48)<k),k=(k<<3)+(k<<1)+(c^48),k1=(k1*10%(mod-1)+(c^48))%(mod-1),c=getchar();}
}

bool Med;
signed main(){
    init(300000);

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}

第二类斯特林数

定义

第二类斯特林数 \begin{Bmatrix}n\\m\end{Bmatrix} 也记作 S(n,m)，表示将 n 个区分的小球放入 m 个相同的盒子且不允许有空盒子的方案数。

递推公式

\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix}

通过组合意义来解释：

当加入一个球时可以将其单独放入一个空盒子中，方案数为 \begin{Bmatrix}n-1\\m-1\end{Bmatrix} ，

也可以将其放入一个非空的盒子中，方案数为 m\begin{Bmatrix}n-1\\m\end{Bmatrix}。

通过加法原理得到递推公式。

通向公式

\begin{Bmatrix}n\\m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!}

可以使用容斥证明该式子：

  1. 将 n 个区分的球放入 m 个区分的盒子中且允许有空盒子的方案为 A_m，
  2. 将 n 个区分的球放入 m 个区分的盒子但不能有空盒子的方案是 B_m。

那么有 :

A_m=m^n=\sum_{i=0}^m\binom{m}{i}B_i

通过二项式反演可得：

\begin{aligned} B_m&=\sum_{i=0}^m(-1)^{m-i}\binom{m}{i}A_i\\&=\sum_{i=0}^m (-1)^{m-i}\frac{m!i^n}{i!(m-i)!} \end{aligned}

因为 B_m 正好 \begin{Bmatrix}n\\m\end{Bmatrix} 的 i! 倍，所以得到通向公式：

\begin{Bmatrix}n\\m\end{Bmatrix}=\sum_{i=0}^m\frac{(-1)^{m-i}i^n}{i!(m-i)!}

第二类斯特林数·行

发现其通向公式

\begin{Bmatrix}n\\i\end{Bmatrix}=\sum\limits_{j=0}^i\frac{(-1)^{i-j}j^n}{j!(i-j)!}\\ =\sum_{j=0}^i\frac{j^n}{j!}.\frac{(-1)^{i-j}}{(i-j)!}

这是一个卷积的形式，所以可以 NTT 求。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=1e6+10;
const int mod=167772161;
const int inf=1e9+10;

template <int mod>
struct modint{

    int val;

    static int norm(const int &x){return x<0?x+mod:x;}

    modint inv()const{
        int a=val,b=mod,u=1,v=0,t;
        while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
        return modint(u);
    }

    modint():val(0){}
    modint(const int &m):val(norm(m)){}
    modint(const ll &m):val(norm(m%mod)){}
    modint operator -()const{return modint(norm(-val));}
    bool operator ==(const modint &x){return val==x.val;}
    bool operator !=(const modint &x){return val!=x.val;}
    bool operator <=(const modint &x){return val<=x.val;}
    bool operator >=(const modint &x){return val>=x.val;}
    bool operator >(const modint &x){return val>x.val;}
    bool operator <(const modint &x){return val<x.val;}
    modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
    modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
    modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
    modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
    modint& operator >>=(const modint &x){return val>>=x.val,*this;}
    modint& operator ^=(const modint &x){return val^=x.val,*this;}
    modint operator >>(const modint &x){return modint(*this)>>=x;}
    modint operator <<(const modint &x){return modint(*this)<<=x;}
    modint& operator /=(const modint &x){return *this*=x.inv();}
    modint operator +(const modint &x){return modint(*this)+=x;}
    modint operator -(const modint &x){return modint(*this)-=x;}
    modint operator *(const modint &x){return modint(*this)*=x;}
    modint operator /(const modint &x){return modint(*this)/=x;}
    modint operator ^(const modint &x){return modint(*this)^=x;}
    friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
    friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<167772161>m16;

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

int rev[Max];

void init(int lim,int len){
    for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
}

m16 ksm(m16 a,int b){
    m16 ans=1;
    for(;b;b>>=1){
        if(b&1)ans=ans*a;
        a=a*a;
    }
    return ans;
}

struct Poly{

    const m16 g=3;
    const m16 invg=(mod+1)/3;
    m16 a[Max];

    m16& operator [](const int x){return a[x];}

    void NTT(int tag,int lim){
        for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
            for(int R=mid<<1,j=0;j<lim;j+=R){
                m16 w=1;
                for(int p=0;p<mid;++p,w=w*W){
                    m16 x=a[p+j],y=a[p+j+mid]*w;
                    a[p+j]=x+y;
                    a[p+j+mid]=x-y;
                }
            }
        }
        if(tag==-1){
            m16 inv=m16(lim).inv();
            for(int i=0;i<lim;++i)a[i]*=inv;
        }
    }

};

Poly a,b;

m16 frac[Max],inv[Max];

bool Med;
signed main(){
    int n=read();
    inv[0]=frac[0]=1;
    for(int i=1;i<=n;++i)frac[i]=frac[i-1]*i;
    inv[n]=frac[n].inv();
    for(int i=n-1;i>=1;--i)inv[i]=inv[i+1]*(i+1);
    for(int i=0;i<=n;++i){
        a[i]=ksm(i,n)*inv[i];
        b[i]=inv[i]*(i%2==0?1:-1);
    }
    int lim=1,len=0;
    while(lim<=n+n)lim<<=1,++len;
    init(lim,len);
    a.NTT(1,lim),b.NTT(1,lim) ;
    for(int i=0;i<lim;++i)a[i]*=b[i];
    a.NTT(-1,lim) ;
    for(int i=0;i<=n;++i)cotu << a[i] << ' ';
    cotu << "\n";

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}

第二类斯特林数·列

先假设每个盒子也要区分，考虑对于一个盒子计算，其生成函数是 G(x)=\sum\limits_{i\geq 1}\frac{x^i}{i!}=(e^x-1)，其第 i 项表示在这个盒子放入了 i 个球。 这样对于 m 个盒子，F(x)=G(x)^m，那么如果有 n 个球，答案就是 [x^n]G(x)^m 由于目前的盒子是区分的，所以要乘以一个 m!，故有 n 个球的答案是 m![x^n]G(x)^m。这样可以多项式快速幂求，都是发现多项式 G(x) 的首项是 0 所以要多项式平移。这里使用的是 EGF 其卷积形式如下 \sum\limits_{i\geq 0}x^i\sum\limits_{j=0}^i\frac{1}{i!(i-j)!}a_ib_j=\sum\limits_{i\geq 0} \frac{x^i}{i!}\sum\limits_{j=0}^i\binom{i}{j}a_ib_j 就是考虑了 a,b 合并的顺序，而最后要乘一个 i! 去消掉 EGF 的 i! 这个系数。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=2e6+10;
const int mod=167772161;
const int inf=1e9+10;

template <int mod>
struct modint{

    int val;

    static int norm(const int &x){return x<0?x+mod:x;}

    modint inv()const{
        int a=val,b=mod,u=1,v=0,t;
        while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
        return modint(u);
    }

    modint():val(0){}
    modint(const int &m):val(norm(m)){}
    modint(const ll &m):val(norm(m%mod)){}
    modint operator -()const{return modint(norm(-val));}
    bool operator ==(const modint &x){return val==x.val;}
    bool operator !=(const modint &x){return val!=x.val;}
    bool operator <=(const modint &x){return val<=x.val;}
    bool operator >=(const modint &x){return val>=x.val;}
    bool operator >(const modint &x){return val>x.val;}
    bool operator <(const modint &x){return val<x.val;}
    modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
    modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
    modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
    modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
    modint& operator >>=(const modint &x){return val>>=x.val,*this;}
    modint& operator ^=(const modint &x){return val^=x.val,*this;}
    modint operator >>(const modint &x){return modint(*this)>>=x;}
    modint operator <<(const modint &x){return modint(*this)<<=x;}
    modint& operator /=(const modint &x){return *this*=x.inv();}
    modint operator +(const modint &x){return modint(*this)+=x;}
    modint operator -(const modint &x){return modint(*this)-=x;}
    modint operator *(const modint &x){return modint(*this)*=x;}
    modint operator /(const modint &x){return modint(*this)/=x;}
    modint operator ^(const modint &x){return modint(*this)^=x;}
    friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
    friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<167772161>m16;

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

m16 ksm(m16 a,int b){
    m16 ans=1;
    for(;b;b>>=1){
        if(b&1)ans=ans*a;
        a=a*a;
    }
    return ans;
}

m16 inv[Max],frac[Max],Inv[Max];

void init(int len) {
    frac[0]=inv[0]=1; 
    for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
    Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
    for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
    const m16 g=3;
    const m16 invg=(mod+1)/3;

    m16 a[Max];
    void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
    void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
    void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
    void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);}
    m16& operator [](const int x){return a[x];}

    void NTT(int tag,int lim){
        for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
            for(int R=mid<<1,j=0;j<lim;j+=R){
                m16 w=1;
                for(int p=0;p<mid;++p,w=w*W){
                    m16 x=a[j+p],y=a[j+p+mid]*w;
                    a[j+p]=x+y;
                    a[j+p+mid]=x-y;
                }
            }
        }
        if(tag==-1){
            m16 inv=m16(lim).inv();
            for(int i=0;i<lim;++i)a[i]*=inv;
        }
    }

}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
    static Poly b[2],c;
    int lim=2,len=1,bas=1;
    int pos=0;
    b[pos][0]=x[0].inv();
    while(bas<=(n<<1)){
        init(lim,len);
        pos^=1;
        for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
        for(int i=0;i<=lim;++i)c[i]=x[i];
        for(int i=lim>>1;i<=lim;++i)c[i]=0;
        b[pos^1].NTT(1,lim);c.NTT(1,lim);
        for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
        b[pos^1].NTT(-1,lim);
        for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
        bas<<=1,lim<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=b[pos][i];
    for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetLn(Poly &y,Poly&x,int n){
    static Poly t1,t2;
    for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
    t1.der(n);GetInv(t2,t2,n);
    int lim=1,len=0;
    while(lim<=n<<1)lim<<=1,++len;
    init(lim,len);
    t1.NTT(1,lim),t2.NTT(1,lim) ;
    for(int i=0;i<lim;++i)t1[i]*=t2[i];
    t1.NTT(-1,lim);t1.inte(n) ;
    for(int i=0;i<n;++i)y[i]=t1[i];
    for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}

void GetExp(Poly &y,Poly &x,int n){
    static Poly s,tmp;s[0]=1;
    int bas=1,lim=2,len=1;
    while(bas<=n<<1){
        for(int i=0;i<bas;++i)tmp[i]=s[i];
        GetLn(tmp,tmp,bas);
        for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
        tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
        for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
        s.NTT(-1,lim) ;
        for(int i=bas;i<lim;++i)s[i]=0;
        lim<<=1,bas<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=s[i];
    for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}

void expow(Poly &y,Poly &x,int n,int k){
    static Poly s;
    for(itn i=0;i<=n;++i)s[i]=x[i];GetLn(s,s,n);
    for(int i=0;i<=n;++i)s[i]*=k;GetExp(s,s,n);
    for(int i=0;i<=n;++i)y[i]=s[i];
    for(int i=0;i<=n;++i)s[i]=0;
}

Poly b;

bool Med;
signed main(){
    int n=read()+1;
    int k=read();
    init(n<<3);
    for(int i=1;i<n;++i)b[i-1]=Inv[i];
    expow(b,b,n,k);
    for(int i=0;i<n;++i){
        if(i<k)cout << 0 << " ";
        else cout << (b[i-k]*Inv[k]*frac[i])<< ' ';
    }

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}
/*

第一类斯特林数

定义

第一类斯特林数 \begin{bmatrix}n\\k\end{bmatrix} 表示将 n 个区分的元素划分为 k 个互不区分的非空轮换中的方案数。

递推公式

\begin{bmatrix}n\\m\end{bmatrix}=\begin{bmatrix}n-1\\m-1\end{bmatrix}+(n-1)\begin{bmatrix}n-1\\m\end{bmatrix}

因为加入一个元素时有以下两种情况：

  1. 加在之前一个轮换中，此时它可以加在 n-1 个元素的右边，故有 (n-1)\begin{bmatrix}n-1\\m\end{bmatrix} 种方案
  2. 新开一个轮换，故有 \begin{bmatrix}n-1\\m-1\end{bmatrix} 种方案

第一类斯特林数·行

考虑像第二类斯特林数

第一类斯特林数·列

可考虑对于每个轮换生成函数为 G(x)=\sum\limits_{i\geq 0}\frac{x^i}i 所有可以直接多项式快速幂求 F(x)=G(x)^m。

#include <bits/stdc++.h>
#define pii pair<int,int>
#define ll long long
#define mk make_pair
#define reaD read
#define raed read
#define haed head
#define cotu cout
#define se second
#define fi first
#define itn int
using namespace std;
bool Mst;
const int Max=2e6+10;
const int mod=167772161;
const int inf=1e9+10;

template <int mod>
struct modint{

    int val;

    static int norm(const int &x){return x<0?x+mod:x;}

    modint inv()const{
        int a=val,b=mod,u=1,v=0,t;
        while(b>0)t=a/b,swap(a-=t*b,b),swap(u-=t*v,v);
        return modint(u);
    }

    modint():val(0){}
    modint(const int &m):val(norm(m)){}
    modint(const ll &m):val(norm(m%mod)){}
    modint operator -()const{return modint(norm(-val));}
    bool operator ==(const modint &x){return val==x.val;}
    bool operator !=(const modint &x){return val!=x.val;}
    bool operator <=(const modint &x){return val<=x.val;}
    bool operator >=(const modint &x){return val>=x.val;}
    bool operator >(const modint &x){return val>x.val;}
    bool operator <(const modint &x){return val<x.val;}
    modint& operator *=(const modint &x){return val=static_cast<int>(1ll*val*x.val%mod),*this;}
    modint& operator <<=(const modint &x){return val=(1ll*val<<x.val)%mod,*this;}
    modint& operator +=(const modint &x){return val=(1ll*val+x.val)%mod,*this;}
    modint& operator -=(const modint &x){return val=norm(1ll*val-x.val),*this;}
    modint& operator >>=(const modint &x){return val>>=x.val,*this;}
    modint& operator ^=(const modint &x){return val^=x.val,*this;}
    modint operator >>(const modint &x){return modint(*this)>>=x;}
    modint operator <<(const modint &x){return modint(*this)<<=x;}
    modint& operator /=(const modint &x){return *this*=x.inv();}
    modint operator +(const modint &x){return modint(*this)+=x;}
    modint operator -(const modint &x){return modint(*this)-=x;}
    modint operator *(const modint &x){return modint(*this)*=x;}
    modint operator /(const modint &x){return modint(*this)/=x;}
    modint operator ^(const modint &x){return modint(*this)^=x;}
    friend std::ostream& operator<<(std::ostream& os,const modint &a){return os<<a.val;}
    friend std::istream& operator>>(std::istream& is,modint &a){return is>>a.val;}
};

typedef modint<167772161>m16;

inline int read(){
    int res=0,v=1;
    char c=getchar();
    while(c<'0'||c>'9'){v=(c=='-'?-1:1);c=getchar();}
    while(c>='0'&&c<='9'){res=(res<<3)+(res<<1)+(c^48);c=getchar();}
    return res*v;
}

m16 ksm(m16 a,int b){
    m16 ans=1;
    for(;b;b>>=1){
        if(b&1)ans=ans*a;
        a=a*a;
    }
    return ans;
}

m16 inv[Max],frac[Max],Inv[Max];

void init(int len) {
    frac[0]=inv[0]=1; 
    for(int i=1;i<=len;++i)frac[i]=frac[i-1]*i;
    Inv[len]=frac[len].inv();for(int i=len-1;i;--i)Inv[i]=Inv[i+1]*(i+1);
    for(int i=1;i<=len;++i)inv[i]=Inv[i]*frac[i-1];
}

int rev[Max];

void init(int lim,int len){for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));}

struct Poly{  //导和积都到n，否则到n-1; 
    const m16 g=3;
    const m16 invg=(mod+1)/3;

    m16 a[Max];
    void out(int n,int x=0){for(int i=0;i<n;++i)cout << a[i].val << " ";if(!x)cout << "\n";}
    void inte(int n){for(int i=n;i>=0;--i)a[i]=a[i-1]*inv[i];a[0]=0;}
    void ReadIn(int n){for(int i=0;i<n;++i)a[i].val=read();}
    void der(int n){for(int i=0;i<=n;++i)a[i]=a[i+1]*(i+1);a[n]=0;}
    m16& operator [](const int x){return a[x];}

    void NTT(int tag,int lim){
        for(int i=0;i<lim;++i)if(i<rev[i])swap(a[i],a[rev[i]]);
        for(int mid=1;mid<lim;mid<<=1){
            m16 W=ksm(tag==1?g:invg,(mod-1)/(mid<<1));
            for(int R=mid<<1,j=0;j<lim;j+=R){
                m16 w=1;
                for(int p=0;p<mid;++p,w=w*W){
                    m16 x=a[j+p],y=a[j+p+mid]*w;
                    a[j+p]=x+y;
                    a[j+p+mid]=x-y;
                }
            }
        }
        if(tag==-1){
            m16 inv=m16(lim).inv();
            for(int i=0;i<lim;++i)a[i]*=inv;
        }
    }

}a;

void GetInv(Poly &y,Poly &x,int n){    //第二个为参数，第一个为返回值
    static Poly b[2],c;
    int lim=2,len=1,bas=1;
    int pos=0;
    b[pos][0]=x[0].inv();
    while(bas<(n<<1)){
        init(lim,len);
        pos^=1;
        for(int i=0;i<bas;++i)b[pos][i]=b[pos^1][i]<<1;
        for(int i=0;i<=lim;++i)c[i]=x[i];
        for(int i=lim>>1;i<=lim;++i)c[i]=0;
        b[pos^1].NTT(1,lim);c.NTT(1,lim);
        for(int i=0;i<=lim;++i)b[pos^1][i]=b[pos^1][i]*b[pos^1][i]*c[i];
        b[pos^1].NTT(-1,lim);
        for(int i=0;i<bas;++i)b[pos][i]-=b[pos^1][i],b[pos^1][i]=0;
        bas<<=1,lim<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=b[pos][i];
    for(int i=0;i<bas;++i)b[0][i]=b[1][i]=c[i]=0;
}

void GetLn(Poly &y,Poly&x,int n){
    static Poly t1,t2;
    for(int i=0;i<n;++i)t1[i]=t2[i]=x[i];
    t1.der(n);GetInv(t2,t2,n);
    int lim=1,len=0;
    while(lim<n<<1)lim<<=1,++len;
    init(lim,len);
    t1.NTT(1,lim),t2.NTT(1,lim) ;
    for(int i=0;i<lim;++i)t1[i]*=t2[i];
    t1.NTT(-1,lim);t1.inte(n) ;
    for(int i=0;i<n;++i)y[i]=t1[i];
    for(int i=0;i<lim;++i)t1[i]=t2[i]=0;
}

void GetExp(Poly &y,Poly &x,int n){
    static Poly s,tmp;s[0]=1;
    int bas=1,lim=2,len=1;
    while(bas<n<<1){
        for(int i=0;i<bas;++i)tmp[i]=s[i];
        GetLn(tmp,tmp,bas);
        for(int i=0;i<bas;++i)tmp[i]=x[i]-tmp[i];
        tmp[0]+=1;init(lim,len);tmp.NTT(1,lim),s.NTT(1,lim) ;
        for(int i=0;i<lim;++i)s[i]=s[i]*tmp[i];
        s.NTT(-1,lim) ;
        for(int i=bas;i<lim;++i)s[i]=0;
        lim<<=1,bas<<=1,++len;
    }
    for(int i=0;i<n;++i)y[i]=s[i];
    for(itn i=0;i<bas;++i)s[i]=tmp[i]=0;
}

void expow(Poly &y,Poly &x,int n,int k){
    static Poly s;
    for(itn i=0;i<n;++i)s[i]=x[i];GetLn(s,s,n);
    for(int i=0;i<n;++i)s[i]*=k;GetExp(s,s,n);
    for(int i=0;i<n;++i)y[i]=s[i];
    for(int i=0;i<n;++i)s[i]=0;
}

Poly b;

bool Med;
signed main(){
    int n=read()+1;
    int k=read();
    init(n);
    for(int i=0;i<n-1;++i)b[i]=inv[i+1];
    expow(b,b,n,k);
    m16 z=frac[k].inv();
    for(int i=0;i<n;++i){
        if(i<k)cout << 0 << " ";
        else cout << (b[i-k]*z*frac[i]).val<< ' ';
    }

    cerr<< "Time: "<<clock()/1000.0 << "s\n";
    cerr<< "Memory: " << (&Mst-&Med)/1000000.0 << "MB\n";
    return 0;
}
/*

*/

生成函数

定义

生成函数又称母函数，是一种形式幂级数，即只关心未知数幂次和常数，不关心未知数取值。

一般的生成函数形式为

F(x)=\sum_{n\geq0}a_nk_n(x)

其中 a_n 为常数，k_n(x)称作核函数，不同核函数会产生不同生成函数。

对于两个生成函数 A(x) 和B(x)，如果满足 A(x) 中a_0\neq 0 并且 A(x)\times B(x)=1 ，那么称 B(x) 是生成函数 A(x) 的逆元。

性质

求导

我们定义 F(x)=\sum\limits_{i\geq0}a_ix^i 的导数是：

F'(x)=\sum_{i\geq 0}(i+1)a_{i+1}x^i=a_1+2a_2x^2+\dots

对于一般的求导式子依然成立：

Maclaurin 公式

麦克劳林展开式泰勒展开在函数 0 点的特殊情况，对于泰勒展开，我们构造一个多项式使其在若干次导下都与需要逼近的函数这么多次导相同（当然式该函数一个点的取值）。

我们设 P(0) 是函数 p(x) 在 0 处的取值，那么有：

P(x)=\sum_{i\geq 0}\frac{P^{(k)}(0)}{i!}x^i=P(0)+P'(0)x+\frac{P''(0)}{2!}x^2\dots

分类

常见生成函数

形如 F(x)=\sum_{n\geq0}a_nx^n 的生成函数称作常见生成函数，又称普通生成函数（OGF）。

当序列a 全部等于1时：

F(x)=1+x+x^2+x^3+\cdots +x^n

我们可以找到其逆元为B(x)=1-x ,所以:

F(x)=\frac{1}{B(x)}\\ \ \ \ \ \ \ \ \ \ \ =\frac{1}{1-x}

常见的普通生成函数有 Fibonacci 和 Catalan:

Fibonacci

其递归式定义为

f_x=\begin{cases} 1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ,x\in \{1,2\}\\ f_{x-1}+f_{x-2},x>2 \end{cases}

我们令F(x)=\sum_{n\geq0}f_nx^n 为f_x的OGF，那么我们可以发现

F(x)=f_0x^0+f_1x^1+f_2x^2+f_3x^3+f_4x^4+\cdots\\ xF(x)=\ \ \ \ \ \ \ \ \ \ f_0x^1+f_1x^2+f_2x^3+f_3x^4+\cdots\\ x^2F(x)= \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ f_0x^2+f_1x^3+f_2x^4+\cdots\\

将其作差我们可以发现

F(x)=x+xF(x)+x^2F(x)\\ F(x)=\frac{x}{1-x-x^2}

设：

F(x)=\frac{c}{1-ax}+\frac{d}{1-bx}\\ =\frac{c+d-(ad+cb)x}{1-(a+b)x+abx^2}=\frac{x}{1-x-x^2}\\ \therefore \begin{cases} ab=-1\\ a+b=1\\ c+d=0\\ ad+cb=-1\\ \end{cases}\\ \therefore \begin{cases} a=\frac{1+\sqrt{5}}{2}\\ b=\frac{1-\sqrt{5}}{2}\\ c=\frac{\sqrt{5}}{5}\\ d=-\frac{\sqrt{5}}{5}\\ \end{cases}\\ \therefore F(x)=\sum_{n\geq0}(c\cdot a^{n-1}+d\cdot b^{n-1})x^n\\ =\sum_{n\geq0}(\frac{\sqrt{5}}{5}\cdot ((\frac{1+\sqrt{5}}{2})^{n-1}-(\frac{1-\sqrt{5}}{2})^{n-1}))x^n\\
Catalan

定义卡特兰数序列 c_n 为：

c_n=\begin{cases} 1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ n\in \{0,1 \} \\ \sum_{i=0}^{n-1} c_ic_{n-i-1} \ \ \ n>1 \end{cases}

同理，我们令C(x)=\sum_{n\geq0}c_nx^n 为 c_x 的生成函数，那么我们可以发现：

C(x)=\sum_{n\geq0}c_nx^n\\ =\sum_{n\geq0}(\sum_{i=0}^{n-1}c_ic_{n-i-1})x^n\\ =1+\sum_{n\geq1}(\sum_{i=0}^{n-1}c_ic_{n-i-1})x^n\\

令 n'=n-1 ,

=1+x\sum_{n'\geq0}(\sum_{i=0}^{n'}c_ic_{n'-i})x^{n'}

我们可以发现，式子括号中的部分类似卷积形式，那么我们可以将其转换为：

=1+xC(x)^2\\ \therefore C(x)=1+xC(x)^2

解出这个一元二次方程（x 为待定系数）：

C(x)=\frac{1\pm \sqrt{1-4x}}{2x}\\

对其分子有理化:

\frac{1\pm \sqrt{1-4x}}{2x}=\frac{1-(1-4x)}{2x(1\mp \sqrt{1-4x})}=\frac{4x}{2x \mp 2x(\sqrt{1-4x})}=\frac{2}{1\mp\sqrt{1-4x}}

我们可以发现，当 x 趋近于 0 时取负号是收敛的，但是取正号是发散的，所以：

\therefore C(x)=\frac{1-\sqrt{1-4x}}{2x}

后面就是广义二项式定理（不会），只知道最终结果为：

c_n=\frac{(_{\ n}^{2n})}{n+1}

一些生成函数的封闭形式