NTT初步

NTT可以理解为带模意义下的FFT。

原根 代替了 单位根

原根 的定义:

$(a,m)=1$,定义 $a$ 对模 $m$ 的质数 $\delta_ m(a)$ 为使得 $a^d\equiv 1(\text{mod} ~m)$ 的最小正整数 $d$,

若 $\delta_ m(a)=\varphi(m)$,则 $a$ 为模 $m$ 的原根。

当模数为 $m=k2^t+1$ 且 $m\in \text{prime}$ 时,$\varphi(m)=m-1$

找出任意一个原根 $g$ ,

$g^{0\sim m-1}$ 两两不相同。

之后将 $\omega_{m}^{i}=g^i$ 套入之后,发现符合所有单位根用到的性质。

于是就有了快速数论变换(NTT)。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<ctime>
#include<map>
#include<cstdlib>
#include<vector> 
#define gc getchar()
#define ll long long
#define ull unsigned long long
#define file(s) freopen(s".in","r",stdin);freopen(s".out","w",stdout)
#define I inline 
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n)) 
using namespace std;
const int N=2e6+5,mod=998244353,W=mod-1;
const ull p0=31,p1=37,G=3;
template<class o>I void qr(o &x)
{
    char c=gc;int f=1;x=0;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=gc;}
    while(c>='0'&&c<='9'){x=x*10+(c^48);c=gc;}
    x*=f;
}
template<class o>I void qw(o x)
{
    if(x<0)x=-x,putchar('-');
    if(x/10)qw(x/10);
    putchar(x%10+48);
}
I int ksm(int a,int b=mod-2){int ans=1;for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;return ans;}
int r[N<<1],cnt;
I void rev(int n)
{
    if(cnt==n)return ;cnt=n;
    for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1|((i&1)?n>>1:0));
}
void dft(int *g,bool op,int n)
{
    rev(n);
    static ull f[N<<1],w[N<<1]={1};
    for(int i=0;i<n;i++)f[i]=((1ll*mod<<5)+g[r[i]])%mod;
    for(int p=2,l=1;p<=n;l=p,p<<=1)
    {
        w[1]=ksm(G,op?W-W/p:W/p);
        for(int i=2;i<l;i++)w[i]=w[i-1]*w[1]%mod;
        for(int i=0;i<n;i+=p)for(int j=0,t;j<l;j++)
            t=w[j]*f[i|j|l]%mod,f[i|j|l]=f[i|j]+mod-t,f[i|j]=f[i|j]+t;
        if(l==(1<<10))for(int i=0;i<n;i++)f[i]%=mod;
    }
    if(op){int t=ksm(n);for(int i=0;i<n;i++)g[i]=f[i]%mod*t%mod;}
    else for(int i=0;i<n;i++)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n){for(int i=0;i<n;i++)f[i]=1ll*f[i]*g[i]%mod;}
int _g1[N<<1];
#define sav _g1
void times(int *f,int *g,int m,int limit)
{
    int n=1;for(;n<(m<<1);n<<=1);
    cpy(sav,g,n);for(int i=m;i<n;i++)sav[i]=0;
    dft(f,0,n);dft(sav,0,n);
    px(f,sav,n);dft(f,1,n);
    for(int i=limit;i<n;i++)f[i]=0;
}
#undef sav
int F[N<<1]; 
int main()
{
    return 0;
}

多项式全家桶

多项式求逆

首先得明白一个概念:

多项式乘法:( $f[i]$ 表示 $x^i$ 的系数)

$(f*g)[i]=\sum_{j=0}^i f[j]g[i-j]$

多项式求逆:

$(f^{-1}*f)[i]=\sum_{j=0}^if[j]f^{-1}[i-j]=[i=0]$

因此多项式求逆不仅仅只是将$f$ 的的每一项系数求逆那么简单。

对于多项式求逆,

有一个很粗暴的方法,就是

$\text{mod }x^n\rightarrow \text{mod}~ x^{n+1}$

优化就采用倍增,

$\text{mod}~x^{n/2}\rightarrow \text{mod}~x^n$

怎么实现呢?

假设现在已经求好了 $R_*=F^{-1}(\text{mod} ~x^{n/2})$

意义为:多项式 $R_*$ 与多项式 $F$ 在不考虑卷积之后的 $x^{n/2}$ 及 $x^{n/2}$ 以上的位置的系数时,$R_*$ 为 $F^{-1}$ ,即 $\forall i\in[0,n/2),(R_**F)[i]=[i=0]$ 。

现在要求: $R=F^{-1}(\text{mod}~x^n)$

显然,$R=R_*(\text{mod } x^{n/2})$

即$\forall i\in[0,n/2),(R-R_*)[i]=0$

然后把,$(R-R_*)$ 这个多项式平方。

就有一个很美妙的性质,

$\forall i\in[n/2,n),(R-R_*)^2[i]=\sum_{j=0}^{i}(R-R_*)[j](R-R_*)[i-j]=0$

即 $(R-R_*)^2=0(\text{mod}~x^n)$ ,这里的 $0$ 表示系数均为 $0$ 的多项式。

$(R-R_*)^2=0(\text{mod}~x^n)$

$R^2-2RR_*+{R_*}^2=0(\text{mod}~x^n)$

同乘 $F$ ,

$R-2R_*+{R_*}^2*F=0(\text{mod}~x^n)$

$R=2R_*-{R_*}^2*F(\text{mod}~x^n)$

这个时候仅需要快速算出 $({R_*}^2*F)$ 这个多项式的每一个系数就可以得出结果了。

由于 ${R_*}^2*F$ 的最高次幂为 $x^{2n-1}$,所以要开两倍。

int _w2[N<<1],_r2[N<<1];
void invp(int *f,int m)
{
    int n=1;for(;n<m;n<<=1);
    #define w _w2
    #define r _r2
    w[0]=ksm(f[0]);
    for(int p=2;p<=n;p<<=1)
    {
        for(int i=0;i<(p>>1);i++)
            r[i]=(w[i]<<1)%mod;
        cpy(sav,f,p);
        dft(w,0,p<<1);px(w,w,p<<1);
        dft(sav,0,p<<1);px(w,sav,p<<1);
        dft(w,1,p<<1);clr(w+p,p);
        for(int i=0;i<p;i++)    
            w[i]=(r[i]+mod-w[i])%mod;
    }cpy(f,w,m);clr(w,n<<1),clr(r,n<<1),clr(sav,n<<1);
    #undef w
    #undef r
}

多项式牛顿迭代

$G(x)$ 是一个比较好算的函数,可以转化为常量,$F$ 是多项式。

若 $G(F)=0(\text{mod} ~x^n)$,且 $G(F_*)=0(\text{mod}~x^{n/2})$

$F=F_*-\frac{G(F_*)}{G'(F_*)}(\text{mod}~x^n)$

证明:

将 $G(F)$ 在关于 $F_*=F(\text{mod}~x^n)$ 处进行泰勒展开,

$G(F)=G(F_*)+\frac{G^{(1)}(F_*)}{1!}(F-F_*)+\frac{G^{(2)}(F_*)}{2!}(F-F_*)^2+\ldots+R_n(F)$

由于 $(F-F_*)$ 在$i\in[0,n/2)$ 项上的系数都为 $0$,故 $(F-F_*)^2$ 在$i\in[0,n)$ 项上的系数都为 $0$。

那么,$(F-F_*)^m,m\ge 3$ 同理。

则 $G(F)=G(F_*)+G^{(1)}(F_*)(F-F_*)(\text{mod}~x^n)$。

由于 $G(F)=0(\text{mod} ~x^n)$,则 $G(F_*)+G^{(1)}(F_*)(F-F_*)=0(\text{mod}~x^n)$ 。

$G^{(1)}(F_*)F=-G(F_*)+G^{(1)}(F_*)F_*(\text{mod}~x^n)$

$F=F_*-\frac{G(F_*)}{G^{(1)}(F_*)}(\text{mod}~x^n)$

有什么用呢?

方便推式子:

比如上文的多项式求逆:

已知 $A*B=1(\text{mod} ~x^n), A$ 是已知多项式。

设 $G(B)=A*B-1=0(\text{mod}~x^n)$ ,即 $G(x)=A x-1$

那么 $G'=A(\text{mod }x^n)$,

由牛顿迭代可知:

$B=B_*-\frac{G(B_*)}{G^{(1)}(B_*)}=B_* -\frac{A*B_*-1}{A}=B_*-B(A*B_*-1)$

因为 $A*B_*-1$ 的低 $n/2$ 位均为 $0$ ,而 $B$ 中的高 $n/2$ 位仅与 $A*B_*-1$ 仅有的高 $n/2$ 位 相乘,

但最终还是被 $x^n$ 模掉。

因此 $B_*$ 在这里是与 $B$ 是等价的。

故 $B=B_*(2-A*B_*)$ 。

多项式开根

$B^2-A=0(\text{mod}~x^n)$

$G(B)=B^2-A=0$,$G(x)=x^2-A$

$G’(x)=2x,G’(B)=2B$

则 $B=B_*-\frac{G(B_*)}{G^{(1)}(B_*)}=B_*-\frac{{B_*}^2-A}{2B_*}=\frac{A+{B_*}^2}{2B_*}(\text{mod }x^n)$

相当于两个$n$次多项式相乘。

int _b1[N<<1],_b2[N<<1];
void sqrp(int *f,int m)
{
    int n=1;for(;n<m;n<<=1);
    #define b1 _b1
    #define b2 _b2
    b1[0]=1;
    for(int p=2;p<=n;p<<=1)
    {
        for(int i=0;i<(p>>1);i++)
            b2[i]=(b1[i]<<1)%mod;
        invp(b2,p);
        dft(b1,0,p);px(b1,b1,p);dft(b1,1,p);
        for(int i=0;i<p;i++)b1[i]=(b1[i]+f[i])%mod;
        times(b1,b2,p,p);
    }cpy(f,b1,m);clr(b1,n),clr(b2,n);
    #undef b1
    #undef b2
}

多项式带余除法

$F(x)=Q(x)*G(x)+R(x)$

其中,$F(x)$ 为 $n$ 次多项式 , $G(x)$ 为 $m$ 次多项式,$Q(x)$ 为$n-m$次多项式, $R(x)$ 为不足 $m-1$ 次多项式。

想办法把余数削掉。

但显然,如果不进行任何变换,是完全没有办法削掉的。

考虑翻转多项式的系数。

自己动手一下,就可以发现 $x^nF(x^{-1})$ 的系数恰好就是原多项式系数的翻转,记为 $F_{r}(x)$。

于是,$x^nF(x^{-1})=x^nQ(x^{-1})*G(x^{-1})+x^nR(x^{-1})$

$F_{r}(x)=(x^{n-m}Q(x^{-1})*x^{m}G(x^{-1}))+x^{n-m+1}x^{m-1}R(x^{-1})$

$~~~~~~~~~=Q_r(x)*G_r(x)+x^{n-m+1}R_r(x)$

之后就发现只要模上 $x^{n-m+1}$ 就可以削去余项,即

$F_r(x)=Q_r(x)*G_R(x)(\text{mod }x^{n-m+1})$

$Q_r(x)=\frac{F_r(x)}{G_r(x)}(\text{mod} ~x^{n-m+1})$

于是发现只要将 $G_r(x)$ 求逆,再进行一次多项式乘法就可以得到 $Q_r(x)$ 这个$n-m$ 次多项式,进而翻转求得 $Q(x)$ 。

再乘$G(x)$,之后就是原多项式之间的加减问题了。

void rever(int *f,int n)
{
    for(int i=0;i<n;i++)sav[i]=f[i];
    for(int i=0;i<n;i++)f[i]=sav[n-i-1];
    clr(sav,n);
}
int _q[N<<1],_t[N<<1]; 
void print(int *f,int n)
{
    for(int i=0;i<n;i++)    
        qw(f[i]),putchar(' ');
    puts("");
}
void div(int *f,int *g,int n,int m)
{
    #define q _q
    #define t _t
    rever(g,m);cpy(q,g,n-m+1);rever(g,m);
    rever(f,n);cpy(t,f,n-m+1);rever(f,n);
    invp(q,n-m+1);
    times(q,t,n-m+1,n-m+1,n-m+1);
    rever(q,n-m+1);print(q,n-m+1);
    times(q,g,n-m+1,m,n);
    for(int i=0;i<m-1;i++)
        q[i]=(f[i]-q[i]+mod)%mod;
    print(q,m-1);
    #undef q
    #undef t
}

多项式ln

$\ln(A(x))=B(x)$

$[\ln(A(x))]'=B'(x)$

$A'(x)\ln'(A(x))=B'(x)$

由 $ln'(x)=\frac{1}{x}$ 可得: $A'(x)A^{-1}(x)=B'(x)$

同时积分,

$\int A'(x)A^{-1}(x)dx=B(x)$

由于是一个多项式,运用一下 $(x^n)'=nx^{n-1}$ 就可以得到积分的解了。

int inv[N<<1];
void deri(int *f,int n){for(int i=1;i<n;i++)f[i-1]=1ll*f[i]*i%mod;f[n-1]=0;}
void init(int lim){inv[1]=1;for(int i=2;i<=lim;i++)inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;}
void integ(int *f,int n){for(int i=n;i;i--)f[i]=1ll*f[i-1]*inv[i]%mod;f[0]=0;} 
int _s3[N<<2]; 
void lnp(int *f,int n)
{
    #define g _s3
    cpy(g,f,n);
    invp(g,n);deri(f,n);
    times(f,g,n,n,n);
    integ(f,n-1);
    clr(g,n);
    #undef g
}

多项式exp

$e^{A(x)}=B(x),A(x)$ 为已知多项式

$G(B(x))=ln(B(x))-A(x)=0,G'(B(x))=B^{-1}(x)$

$B_*(x)=B(x)(\text{mod }x^{n/2})$

$B(x)=B_*(x)-\frac{G(B_*(x))}{G^{(1)}(B_*(x))}(\text{mod }x^n)$

$B(x)=B_*(x)-B_*(x)(ln(B_*(x))-A(x))=B_*(x)(1-ln(B_*(x))+A(x))$

int _s1[N<<1],_s2[N<<1];
void Exp(int *f,int m)
{
    int n=1;for(;n<m;n<<=1);
    #define s _s1
    #define s2 _s2
    s[0]=1;
    for(int p=2;p<=n;p<<=1)
    {
        cpy(s2,s,p>>1);lnp(s2,p);
        for(int i=0;i<p;i++)
            s2[i]=(f[i]-s2[i]+mod)%mod;
        s2[0]=(s2[0]+1)%mod;
        times(s,s2,p,p,p);
    }cpy(f,s,m);clr(s,n),clr(s2,n);
    #undef s
    #undef s2
}

多项式快速幂

$A(x)^k=e^{ln(A(x))*k}$,$A(x)$ 已知。

仅需要将 $F(x)=k*\ln(A(x))$ 求出来

再套入 $B(x)=e^{F(x)}$ 的多项式exp即可

int main()
{
    int n,k=0;qr(n);scanf("%s",str+1);int m=strlen(str+1);init(n);
    for(int i=1;i<=m;i++)
        k=(1ll*k*10+(str[i]^48))%mod; 
    for(int i=0;i<n;i++)qr(F[i]);
    lnp(F,n);
    for(int i=0;i<n;i++)F[i]=1ll*F[i]*k%mod;
    Exp(F,n);
    for(int i=0;i<n;i++)
        qw(F[i]),putchar(' ');
    puts("");
    return 0;
}