浅谈$NTT$

Posted 2021-03-13 wxq1229

(NTT)，快速数论变换，可以理解为带模数的FFT。

---

原根 & 阶

先来补一点数论。（这里讲的应该很少，都是针对(ntt)胡的，具体的话可以去看《初等数论》那本小黄书）。

阶（指数）

如果(m > 1, (a,m) = 1)，那么必有整数(d)，使得下面这个柿子成立
[ a^d equiv 1 (mod m) ]
我们称令这个柿子成立的最小的(d)叫(a mod m)的 阶（或指数），记作(delta_m(a))。

一些性质

• ((1).) 因为(a,m)互质，所以根据欧拉定理(a^{varphi(m)} equiv 1 (mod m))，立即得到(delta_m(a) le varphi(m))

• $(2). $ 若(d_0)也满足(a^{d_0} equiv 1 (mod m))，那么必有(delta_m(a) mid d_0)，这个显然。
• $(3). $ 由((1))和((2))可以立即推出(delta_m(a) mid varphi(m))

• (cdots)

---

原根

若(m > 1, (a,m) = 1)，且(delta_m(a) = varphi(m))，那么则称(a)是(mod m)的一个原根。

特别(m)是质数(p)的时候若(delta_p(a) = p-1)，则(a)是原根。

性质

设(g)是(m)的原根，那么(g, g^2,g^3,cdots,g^{varphi(m)})一定是在(mod m)意义下与(m)互质的(varphi(m))个数的一个排列。这个根据定义也显然。

特别的当(m)是一个质数(p)的时候，有(g, g^2, g^3, cdots, g^{p-1})在(mod p)意义下是(1...p-1)这(p-1)个数的一个排列。

偷偷告诉你，原根一般都很小，默认<=100就好了

那对于一个素数，如何找原根呢？

直接暴力，枚举(g = 1...100)，然后(check)，根据阶性质((1))我们在枚举(d)的时候可以枚举(p-1)的因子，注意不要枚举到(p-1)。

#include <bits/stdc++.h>
using namespace std;
int p;
int fpow(int a,int b,int mod=p){
    int ret=1; for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) ret=1ll*ret*a%mod;
    return ret;
}
int chk(int g){
    for (int i=2;i*i<=p-1;i++) // !!从2开始
        if ((p-1)%i==0&&(fpow(g,i)==1||fpow(g,(p-1)/i)==1)) return 0;
    return 1;
}
int main(){
    scanf("%d",&p);
    for (int i=1;i<=100;i++)
        if (chk(i)){printf("G: %d
",i); return 0;}
    return 1;
    return 0;
}

安利一个巨佬原根表

---

NTT

P.S. 以下的原根都是在一个形如(k2^t + 1)的质数(p)的同余系下定义的。（(t)足够大，一般(t ge 20)）

设(g)是(mod p)的原根，那么(g, g^2, cdots, g^{p-1})在(mod p)意义下是(p-1)个不同的数。

---

回去看一下FFT，随便交一发就(T)了，你看这个FFT，他就是逊啦！。

发现(FFT)还有很多不足之处

• 复数运算，浮点数乘法，自带无穷大常数

• 多次调用三角函数，常数大，还掉精度

• 复数乘法异常难背，还要手算一遍

• (cdots)

总之我们得想办法把频繁的复数运算和单位根去掉。

在颓(FFT)的时候我们用到了单位根的几个性质

• 1. (omega_{n}^{k} = (omega_{n}^{1})^k)
• 1. (omega_{n}^{k} = omega_{n}^{k mod n})
• 1. (omega_{n}^{k + n/2} = -omega_{n}^{k})
• 1. (omega_{2n}^{2k} = omega_{n}^{k})

当然根据单位根的定义，(omega_{n}^{0}, omega_{n}^{1}, cdots, omega_{n}^{n-1})必须是(n-1)个不同的数。

结论就是我们在(mod p)意义下，我们用原根(g)的((p-1)/n)次方代替(omega_{n}^{1})就对了，即(omega_{n}^1 = g^{(p-1)/n})。

我们来证一下性质（一下运算均是在(mod p)意义下的）。（1. 就是定义，即(omega_{n}^{k} = (omega_{n}^1)^k = g^{k(p-1)/n}),我们不用证）。

• 1. (omega_{n}^{k} = omega_{n}^{k \\% n})

     因为(p)是质数，有小费马(a^{p-1} equiv 1 (mod p))，立刻得到(g^{(p-1)k/n} = g^{(p-1)k/n \\% (p-1)})

     这时候并不能直接取模，因为(k/n)不一定是整数。我们令({ x })表示在(mod p)意义下(x)的小数部分。

     有((p-1)k/n \\% (p-1) = (p-1){ k/n } = (p-1)(k \\% n) /n)，代上去即证。

• 1. (omega_{n}^{k+n/2} = -omega_{n}^{k})

     根据定义，有(omega_{n}^{k+n/2} = omega_{n}^{k} cdot omega_{n}^{n/2})

     由于((omega_{n}^{n/2})^2 = omega_{n}^{n} = 1)，所以(omega_{n}^{n/2})在(mod p)意义下只可能是(1)或(-1)，又因为(omega_{n}^{n/2} ot= omega_{n}^{0})，所以(omega_{n}^{n/2} equiv -1 ( mod p))。代上去就好了。

• 1. (omega_{2n}^{2k} = omega_{n}^{k})

     这个应该是最好证的。。。

     把(g^{(p-1)/n})代上去，(g^{(p-1)2k/2n} = g^{(p-1)k/n})。没了。

综上所述，在(mod p)意义下用(g^{(p-1)/n})代替(omega_{n}^{1})一点毛病都没有，在(IDFT)的时候用((omega_{n}^{1})^{-1})代替(omega_{n}^{-1})，算逆元就好了。

所以再来看一下板子题 才发现这里系数都<10

我们取一个信仰膜数(p = 998244353)，然后在(mod p)意义下做(NTT)，没了。

#include <bits/stdc++.h>
using namespace std;
inline int read(){
    char ch=getchar();int x=0,f=1;
    while (!isdigit(ch)){f=ch=='-'?-f:f;ch=getchar();}
    while (isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
const int N=4e6+10,P=998244353,G=3,IG=332748118;
int fpow(int a,int b,int mod=P){
    int ret=1; for (;b;b>>=1,a=1ll*a*a%P) if (b&1) ret=1ll*a*ret%P;
    return ret;
}
#define add(x,y) ((x)+(y)>=P?(x)+(y)-P:(x)+(y))
#define sub(x,y) ((x)-(y)<0?(x)-(y)+P:(x)-(y))
int rev[N];
void ntt(int *f,int n,int flg){
    for (int i=0;i<n;i++) if (i<rev[i]) swap(f[i],f[rev[i]]);
    for (int len=2;len<=n;len<<=1){
        int w1=fpow(flg==1?G:IG,(P-1)/len); // w^{-1} = (w^1)^{-1}
        for (int i=0;i<n;i+=len){
            int w=1;
            for (int j=i;j<(len>>1)+i;j++){
                int tmp=1ll*w*f[j+(len>>1)]%P;
                f[j+(len>>1)]=sub(f[j],tmp);
                f[j]=add(f[j],tmp);
                w=1ll*w*w1%P;
            }
        }
    }
}
int f[N],g[N];
int main(){
    int n=read(),m=read();
    for (int i=0;i<=n;i++) f[i]=read();
    for (int i=0;i<=m;i++) g[i]=read();
    int limit=1; while (limit<=n+m)limit<<=1;
    for (int i=0;i<limit;i++) rev[i]=rev[i>>1]>>1|((i&1)?limit>>1:0);
    ntt(f,limit,1),ntt(g,limit,1);
    for (int i=0;i<limit;i++) f[i]=1ll*f[i]*g[i]%P;
    ntt(f,limit,-1); int ivn=fpow(limit,P-2);
    for (int i=0;i<=m+n;i++) printf("%d ",1ll*f[i]*ivn%P);
    return 0;
}

提交记录，总共(1.7s)，nice!