NTT学习笔记

NTT 的中文名称为快速数论变换，英文全称为 Number-theoretic transform，是 FFT 在数论基础上的实现，可以避免 FFT 的复数运算中的精度损失。

NTT 实际上是将 FFT 中的复数运算巧妙的替换为了整数运算。

原根

若一个正整数 和一个质数 满足对于所有 ， 的值互不相同。

NTT原理

实际上，NTT 就是选取了适当的模数以及任何一个这个质数的原根，并用一些整数代替了复数（被代替的就是上面的 ）。

直接这样说估计没人懂，举个例子。

就以最常见的 NTT 模数 和它的原根（之一） 来举例。

首先一个比较重要的地方，就是 ，这是 NTT 模数的关键之一，这个模数减一必须是 的较大次幂的倍数。

假设现在要进行 NTT 的多项式共有 项（和 FFT 一样， 是 的整数次幂），那么用哪个数字代替 比较合适呢？

不难想到，用 代替 ，同理对于其他 ，用 代替 ，这样就大功告成。

此外将 FFT 的其他操作（加法，乘法，减法等）改成模意义下的运算，所有的复数都改成整数，就可以了。

同理，假如使用一个质数 与它的原根 进行 NTT 时，使用 代替 即可。

正确性证明

不难发现，FFT 能够快速运算的运算在于复数的那三个性质。

同时，这也是 NTT 可以快速计算的关键，因为我们刚刚对于 的取值也满足这些性质！

下面我们来一条一条看。 证明： 注意，对于所有质数 都有 。

---

证明：

---

证明：

注意到一定会有

所以这个式子的值为 。

额， 似乎不算很显然，我们来证一下。

首先，根据原根的性质，我们知道对于 ，肯定恰好存在一个 使得 。

那么我们就知道，一定会有 。

同时，我们知道，一定恰好存在一个 使得 （还是原根的性质），而我们知道 。

所以就会有 。

这个方程一共有两个解， 或 。

显然， 不满足条件，因为 。

所以一定会有 。

NTT 就说完了。

代码

#include<algorithm>
#include<cstdio>
using namespace std;
using ll=long long;
const int mod=998244353,G=3;
int wh[2097152];
int f[2097152];
int g[2097152];
int n,m;
ll qow(ll a,int b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
inline void NTT(int *f,int n,bool dft=true)
{
    for(int i=1;i<n;i++) if(i<wh[i]) swap(f[i],f[wh[i]]);
    for(int len=1;(1<<len)<=n;len++)
    {
        int step,num=1<<len;
        if(!dft) step=qow((mod+1)/G,(mod-1)/(1<<len));
        else step=qow(G,(mod-1)/(1<<len));
        for(int l=0,mid=num/2;l<n;l+=num,mid+=num) for(int i=0,s=1;i<num/2;i++,s=(ll)s*step%mod)
        {
            ll now=(ll)f[i+mid]*s%mod;
            f[i+mid]=(f[i+l]-now+mod)%mod;
            f[i+l]=(f[i+l]+now)%mod;
        }
    }
}
inline void NTT(int *f,int *g,int n)
{
    for(int i=1;i<n;i++) wh[i]=(wh[i>>1]>>1)+(i&1?n/2:0);
    NTT(f,n),NTT(g,n);
    for(int i=0;i<n;i++) f[i]=((ll)f[i]*g[i])%mod;
    NTT(f,n,false);ll num=qow(n,mod-2);
    for(int i=0;i<n;i++) f[i]=num*f[i]%mod;
}
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++) scanf("%d",f+i);
    for(int i=0;i<=m;i++) scanf("%d",g+i);
    //两个多项式乘起来一共有n+m+1项
    m+=n,n=1;for(;n<=m;n<<=1);
    //找到第一个大于(原来的)n+m的 2的次幂。
    NTT(f,g,n);
    for(int i=0;i<=m;i++) printf("%d ",f[i]);
    return 0;
}

提交记录

加上快读和快写还能快一点。

---

下面是几个模板呀。

多项式求逆

P4238 【模板】多项式乘法逆

设原多项式为 ，它的逆为 。

设 ，它的逆为 。

则一定会有以下几个性质： 则由 可得： 两边同时乘上 。（其实只有一边） 也就是说，只要我们知道了 ，就可以在 的时间内把 算出来（NTT 把上面的式子卷一下）。

时间复杂度 。

代码

#include<algorithm>
#include<cstdio>
using namespace std;
const int mod=998244353,G=3,RG=332748118;
using ll=long long;
inline int read()
{
    int s=0,w=1;char ch;
    while((ch=getchar())>'9'||ch<'0') if(ch=='-') w=-1;
    while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
    return s*w;
}
inline void write(int num)
{
    if(num>9) write(num/10);
    putchar(num-num/10*10+'0');
}
inline ll qow(ll a,int b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
int wh[262144];
int a[262144];
int b[262144];
int f[262144];
int g[262144];
int m;
inline void NTT(int *f,int n,bool dft=true)
{
    for(int i=1;i<n;i++) if(i<wh[i]) swap(f[i],f[wh[i]]);
    for(int len=1;(1<<len)<=n;len++)
    {
        int step,num=1<<len;
        step=qow(dft?G:RG,(mod-1)/(1<<len));
        for(int l=0,mid=num/2;l<n;l+=num,mid+=num) for(int i=0,s=1;i<num/2;i++,s=(ll)s*step%mod)
        {
            ll now=(ll)f[i+mid]*s%mod;
            f[i+mid]=(f[i+l]-now+mod)%mod;
            f[i+l]=(f[i+l]+now)%mod;
        }
    }
}
inline void NTT(int *f,int *g,int n)
{
    for(int i=1;i<n;i++) wh[i]=(wh[i>>1]>>1)+(i&1?n/2:0);
    NTT(f,n),NTT(g,n);
    for(int i=0;i<n;i++)
    {
        f[i]=(2-(ll)f[i]*g[i])%mod;
        f[i]=(f[i]+mod)%mod;
        f[i]=(ll)f[i]*g[i]%mod;
    }
    NTT(f,n,false);ll num=qow(n,mod-2);
    for(int i=0;i<n;i++) f[i]=num*f[i]%mod;
}
void getrev(int m)
{
    if(m==1)
    {
        b[0]=qow(a[0],mod-2);
        return ;
    }
    getrev((m+1)/2);
    int n=1;for(;n<=m*2;n<<=1);
    for(int i=0;i<m;i++) f[i]=a[i],g[i]=b[i];
    for(int i=m;i<n;i++) f[i]=g[i]=0;
    NTT(f,g,n);
    for(int i=0;i<m;i++) b[i]=f[i];
}
int main()
{
    m=read();
    for(int i=0;i<m;i++) a[i]=read();
    getrev(m);
    for(int i=0;i<m;i++) printf("%d ",b[i]);
    return 0;
}

提交记录

任意模数多项式乘法

P4245 【模板】任意模数多项式乘法

我们可以发现，这道题里面的模数可能不能满足之前我们对模数的要求，甚至不是一个质数，所以我们不能直接进行 NTT。

这里我使用的是三模数 NTT。

具体做法就是我们找三个可以进行的 NTT 的模数，然后做三遍 NTT，最后使用 CRT 求出答案。

这里我使用的三个模数分别为 ，，。

直接做就可以了。

然而最后的 CRT 还是有一些细节的()，因为模数个数比较少，所以我们可以手颓柿子：

首先我们知道这些：

所以一定存在三个自然数 使得 。 这样我们就可以求出 的值了，不妨记作 。 令 ，则存在 使得 。 和 都知道了，我们就可以求出 了。

当然，你也可以使用 __int128来做。

中间的取模还是有很多细节的，具体可以看代码。

代码

#include<algorithm>
#include<cstdio>
using namespace std;
using ll=long long;
const int mod1=998244353;
const int mod2=985661441;
const int mod3=1004535809;
const int RG1=332748118;
const int RG2=328553814;
const int RG3=334845270;
const int G=3;
inline int read()
{
    int s=0,w=1;char ch;
    while((ch=getchar())>'9'||ch<'0') if(ch=='-') w=-1;
    while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
    return s*w;
}
inline void write(int a)
{
    if(a>9) write(a/10);
    putchar(a-a/10*10+'0');
}
inline ll qow(ll a,int b,const int mod)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
int wh[262144];
ll f1[262144];
ll f2[262144];
ll f3[262144];
ll g1[262144];
ll g2[262144];
ll g3[262144];
int n,m,p;
inline void NTT(ll *f,const int mod,const int RG,bool dft=true)
{
    for(int i=1;i<n;i++) if(i<wh[i]) swap(f[i],f[wh[i]]);
    for(int len=1;(1<<len)<=n;len++)
    {
        ll step=qow(dft?G:RG,(mod-1)/(1<<len),mod);int num=1<<len;
        for(int l=0,mid=num>>1;l<n;l+=num,mid+=num)
        {
            ll s=1;
            for(int i=0;i<(num>>1);i++)
            {
                ll now=f[mid+i]*s%mod;
                f[mid+i]=(f[l+i]-now+mod)%mod;
                f[l+i]=(f[l+i]+now)%mod;
                s=s*step%mod;
            }
        }
    }
    if(!dft)
    {
        ll num=qow(n,mod-2,mod);
        for(int i=0;i<n;i++) f[i]=f[i]*num%mod;
    }
}
inline void NTT(ll *f,ll *g,const int mod,const int RG)
{
    NTT(f,mod,RG),NTT(g,mod,RG);
    for(int i=0;i<n;i++) f[i]=f[i]*g[i]%mod;
    NTT(f,mod,RG,false);
    // for(int i=0;i<n;i++) printf("%lld ",f[i]);
    // printf("\n");
}
inline ll get(ll a1,ll a2,ll a3)
{
    ll k1=(a2-a1+mod2)*qow(mod1,mod2-2,mod2)%mod2;
    ll a4=(k1*mod1+a1);
    ll k4=(a3-a4%mod3+mod3)%mod3*qow((ll)mod1*mod2%mod3,mod3-2,mod3)%mod3;
    return (k4%p*mod1%p*mod2%p+a4%p+p)%p;
}
int main()
{
    n=read(),m=read(),p=read();
    for(int i=0;i<=n;i++) f1[i]=f2[i]=f3[i]=read();
    for(int i=0;i<=m;i++) g1[i]=g2[i]=g3[i]=read();
    m+=n,n=1;for(;n<=m;n<<=1);
    for(int i=1;i<n;i++) wh[i]=(wh[i>>1]>>1)+(i&1?n/2:0);
    NTT(f1,g1,mod1,RG1),NTT(f2,g2,mod2,RG2),NTT(f3,g3,mod3,RG3);
    for(int i=0;i<=m;i++) printf("%lld ",get(f1[i],f2[i],f3[i]));
    return 0;
}

提交记录

把三个数字写成一个结构体，一起 NTT 还能更快亿点，真的是亿点。

代码

#include<iostream>
#include<cstdio>
using namespace std;
using ll=long long;
ll qow(ll a,int b,int mod)
{
    a%=mod;
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
const ll mod1=998244353;
const ll mod2=1004535809;
const ll mod3=469762049;
const ll inv1=qow(3,mod1-2,mod1);
const ll inv2=qow(3,mod2-2,mod2);
const ll inv3=qow(3,mod3-2,mod3);
inline int read()
{
    int s=0,w=1;char ch;
    while((ch=getchar())>'9'||ch<'0') if(ch=='-') w=-1;
    while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
    return s*w;
}
inline ll lread()
{
    ll s=0,w=1;char ch;
    while((ch=getchar())>'9'||ch<'0') if(ch=='-') w=-1;
    while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
    return s*w;
}
int p;
struct node
{
    ll A,B,C;
    node(ll a,ll b,ll c){ A=a,B=b,C=c;}
    node(ll a=0){ A=B=C=a;}
    void operator = (int a){ A=a%mod1,B=a%mod2,C=a%mod3;}
    node operator + (node b){ return node((A+b.A)%mod1,(B+b.B)%mod2,(C+b.C)%mod3);}
    node operator - (node b){ return node((A+mod1-b.A)%mod1,(B+mod2-b.B)%mod2,(C+mod3-b.C)%mod3);}
    node operator * (node b){ return node(A*b.A%mod1,B*b.B%mod2,C*b.C%mod3);}
    ll get()
    {
        ll k1=(B-A+mod2)%mod2*qow(mod1,mod2-2,mod2)%mod2;
        ll x4=A+k1*mod1;
        ll k4=(C-x4%mod3+mod3)%mod3*qow(mod1*mod2,mod3-2,mod3)%mod3;
        return (k4%p*mod1%p*mod2%p+x4)%p;
    }
}f[300001],g[300001];
ostream & operator << (ostream &out,node a)
{
    out<<'('<<a.A<<','<<a.B<<','<<a.C<<')';
    return out;
}
void swap(node &a,node &b){ node t=a;a=b,b=t;}
int wh[300001];
int a[300001];
int b[300001];
int n,m;
void ntt(node *f,int dft)
{
    for(int i=0;i<n;i++) if(i<wh[i]) swap(f[i],f[wh[i]]);
    for(int len=2;len<=n;len<<=1)
    {
        node step(qow(dft==1?3:inv1,(mod1-1)/len,mod1),qow(dft==1?3:inv2,(mod2-1)/len,mod2),qow(dft==1?3:inv3,(mod3-1)/len,mod3));
        int half=len>>1;
        for(int l=0;l<n;l+=len)
        {
            node now(1);
            for(int i=l;i<l+half;i++)
            {
                node num=now*f[i+half];
                f[i+half]=f[i]-num;
                f[i]=f[i]+num;
                now=now*step;
            }
        }	
    }
    if(!dft)
    {
        ll num1=qow(n,mod1-2,mod1),num2=qow(n,mod2-2,mod2),num3=qow(n,mod3-2,mod3);
        node num(num1,num2,num3);
        for(int i=0;i<n;i++) f[i]=f[i]*num;
    }
}
void do_ntt(node *f,node *g)
{
    for(int i=0;i<n;i++) wh[i]=(wh[i>>1]>>1)|((i&1)?n>>1:0);
    ntt(f,1),ntt(g,1);
    for(int i=0;i<n;i++) f[i]=f[i]*g[i];
    ntt(f,0);
}
int main()
{
    n=read(),m=read(),p=read();
    for(int i=0;i<=n;i++) a[i]=read();
    for(int i=0;i<=m;i++) b[i]=read();
    for(m+=n,n=1;n<=m;n<<=1);
    for(int i=0;i<n;i++) f[i]=a[i],g[i]=b[i];
    do_ntt(f,g);
    for(int i=0;i<=m;i++) printf("%lld ",f[i].get());
    return 0;
}

提交记录

分治NTT

P4721 【模板】分治 FFT

虽然写的是 FFT，但实际上是分治 NTT...

我们要求出一个满足条件的多项式 ，但是 的值和 有关。

考虑进行分治。

假设我们现在分治到了一个区间 ，下一步肯定是分治 这个区间。

然后呢？我们需要快速处理区间 向 的贡献。

实际上 NTT 卷一下就可以了。

时间复杂度 。

代码

#include<algorithm>
#include<cstdio>
using namespace std;
using ll=long long;
const int mod=998244353;
const int RG=332748118;
const int G=3;
inline int read()
{
    int s=0,w=1;char ch;
    while((ch=getchar())>'9'||ch<'0') if(ch=='-') w=-1;
    while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
    return s*w;
}
inline void write(int a)
{
    if(a>9) write(a/10);
    putchar(a-a/10*10+'0');
}
inline ll qow(ll a,int b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
int wh[262144];
int a[262144];
int b[262144];
ll f[262144];
ll g[262144];
int n,m;
inline void NTT(ll *f,int n,bool dft=true)
{
    for(int i=1;i<n;i++) if(i<wh[i]) swap(f[i],f[wh[i]]);
    for(int len=1;(1<<len)<=n;len++)
    {
        ll step=qow(dft?G:RG,(mod-1)/(1<<len));int num=1<<len;
        for(int l=0,mid=num>>1;l<n;l+=num,mid+=num)
        {
            ll s=1;
            for(int i=0;i<(num>>1);i++)
            {
                ll now=f[mid+i]*s%mod;
                f[mid+i]=(f[l+i]-now+mod)%mod;
                f[l+i]=(f[l+i]+now)%mod;
                s=s*step%mod;
            }
        }
    }
    if(!dft)
    {
        ll num=qow(n,mod-2);
        for(int i=0;i<n;i++) f[i]=f[i]*num%mod;
    }
}
inline void NTT(int n)
{
    for(int i=1;i<n;i++) wh[i]=(wh[i>>1]>>1)+(i&1?n>>1:0);
    NTT(f,n),NTT(g,n);
    for(int i=0;i<n;i++) f[i]=f[i]*g[i]%mod;
    NTT(f,n,false);
}
inline void run(int l,int r)
{
    if(l==r||l>=m) return ;
    int mid=(l+r)>>1,len=r-l+1;
    run(l,mid);

    for(int i=l;i<=mid;i++) f[i-l]=b[i];
    for(int i=1;i<=len;i++) g[i-1]=a[i];
    NTT(len<<1);
    for(int i=mid+1;i<=r;i++) b[i]=(b[i]+f[i-l-1])%mod;
    for(int i=0;i<(len<<1);i++) f[i]=g[i]=0;

    run(mid+1,r);
}
int main()
{
    m=read();
    for(int i=1;i<m;i++) a[i]=read();
    n=1;for(;n<m;n<<=1);
    b[0]=1,run(0,n-1);
    for(int i=0;i<m;i++) write(b[i]),putchar(' ');	
    return 0;
}

提交记录

例题

hdu7162 和 题解