NTT算法小结

从理论上说，经过人们优化的FFT已经十分优秀，能够处理大部分的多项式乘法，但是有的时候仍然会出现下面的情况：

1）常数仍然比较大

2）在进行与整数有关的FFT时，发现得到的结果是一堆诡异的数，你需要不停的和精度搏斗

那么在这时，你就需要学会快速数论变换（NTT）

前置芝士

快速傅里叶变换

你可以上网百度，或者看我的博客

阶与原根

我们由欧拉定理可以知道，对于任意的正整数(a、m)，如果满足(gcd(a,m)=1)，就有(a^{varphi(m)}equiv 1(mod m))

但我们发现，还有一些数满足(a^pequiv 1(mod m))且(p< varphi(m))，因此人们定义了阶

设(m>1)，且(gcd(a,m)=1)，则满足(a^pequiv 1(mod m))的最小正整数(p)成为(a)对模(m)的阶，记作(delta_m(a))

于是就会有(delta_m(a)|p)，充分性很显然，我们证一下必要性

我们设(p=delta_m(a)·q+r)（其中(0leq r < delta_m(a))）

那么(a^p=a^{delta_m(a)·q}·a^requiv a^requiv 1(mod m))

由(p)是最小正整数知(r=0)，所以(delta_m(a)|p)

然后你就有了(delta_m(a)|varphi(p))

好吧这个性质并没有什么用

由上面的欧拉定理，我们不难理解数学家们为什么搞出这么一个蛋疼的定义——原根

如果(delta_m(a)=varphi(m))，那么称(a)是模(m)的一个原根，为了下面表述的方便我们将它记作(g)

我们发现原根有一个这样的性质：(g^0,g^1,g^2,cdots,g^{varphi(m)-1})构成了一个模(m)的完全剩余系

证明：考虑反证法，即假设c存在(i,j)((i>j))满足(g^iequiv g^j(mod m))

​ 两边同时除以(g^j)，有(g^{i-j}equiv 1)，而很明显(i-j<varphi(m))，这与原根的定义相矛盾

​ 于是性质得证

性质验证

在FFT中，我们使用单位根的原因就是单位根满足的一些性质可以加速计算，如果原根也满足的话，那么我们在计算时可以直接替换

性质1

(w_n^0,w_n^1,w_n^2,cdots,w_n^{n-1})两两不同

这在上面已经得到了证明

性质2

(w_{2n}^{2p}=w_n^p)

如果(w_{n}=g^p),那么就应该有(w_{2n}=g^{frac{p}{2}})，他们在乘上两倍的次幂之后值相等

性质3

(w_{n}^{frac{n}{2}+p}=-w_n^p)

因为有((g^{frac{n}{2}})^2equiv 1)，为了保持原根的定义，就会有(g^{frac{n}{2}}equiv -1(mod n))

带回去运算即可

综上所述，原根满足原来的单位根所具有的性质，因此我们可以考虑用原根来代替单位根

实际运用

一点注意事项

1、原根的话在模数不确定的情况下需要自己求，不过如果模数是(998244353)或者(1004535809)的话，它们的原根是3

2、注意在IDFT的时候，原来直接除以的地方要换做求逆元

代码

#include<iostream>
#include<string>
#include<string.h>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
using namespace std;
#define rep(i,a,b) for (i=a;i<=b;i++)
typedef long long ll;
#define maxd 998244353
const double pi=acos(-1.0);
#define int long long 
ll n,m,a[5005000],b[5005000];
int lim=1,r[5005000];

int qpow(int x,int y)
{
    int ans=1,sum=x;
    while (y)
    {
        int tmp=y%2;y/=2;
        if (tmp) ans=(1ll*ans*sum)%maxd;
        sum=(1ll*sum*sum)%maxd;
    }
    return ans;
}

void ntt(int lim,ll *a,int typ)
{
    int i;
    for (i=0;i<lim;i++)
        if (i<r[i]) swap(a[i],a[r[i]]);
    int mid;
    for (mid=1;mid<lim;mid<<=1)
    {
        int gn=qpow(3,(maxd-1)/(mid<<1));
        int sta,len=mid<<1,j;
        for (sta=0;sta<lim;sta+=len)
        {
            int g=1;
            for (j=0;j<mid;j++,g=(g*gn)%maxd)
            {
                int x1=a[j+sta],y1=(g*a[j+sta+mid])%maxd;
                a[j+sta]=(x1+y1)%maxd;
                a[j+sta+mid]=(x1-y1+maxd)%maxd;
            }
        }
    }
    if (typ==-1) reverse(&a[1],&a[lim]);
}

int read()
{
    int x=0,f=1;char ch=getchar();
    while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
    while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
    return x*f;
}

signed main()
{
    n=read();m=read();int i,cnt=0;
    for (i=0;i<=n;i++) a[i]=read();
    for (i=0;i<=m;i++) b[i]=read();
    while (lim<=n+m) {lim<<=1;cnt++;}
    for (i=0;i<=lim;i++)
        r[i]=((r[i>>1]>>1)|((i&1)<<(cnt-1)));
    ntt(lim,a,1);
    ntt(lim,b,1);
    for (i=0;i<=lim;i++) a[i]=(a[i]*b[i])%maxd;
    ntt(lim,a,-1);
    int tmp=qpow(lim,maxd-2);
    for (i=0;i<=n+m;i++) 
    {
        a[i]=(a[i]*tmp)%maxd;
        printf("%lld ",a[i]);
    }
    return 0;
}