学习笔记 - NTT(快速数论变换)

和上一篇 FFT 一样,其实还是复习笔记。


Part0. 参考资料

Tiw_Air_OAO 的 多项式乘法总结

Hint. 因为我上一篇博客写了 FFT,关于 FFT 的东西就不再赘述……


Part1. From FFT to NTT

FFT虽然跑得快,但是因为是浮点数运算,终究还是有 精度、常数 等问题。 - Tiw_Air_OAO

而且处理模数还比较麻烦……

于是就引入了 NTT —— 在模的意义下(多项式的系数对某数取模)将多项式和点值表达之间等价转化。分析 FFT 必须用到 double,会卡精度的原因,发现主要是 单位复根 在大多数情况下是一个涉及浮点数的复数。

NTT 即考虑在模的意义下,用别的东西来替换掉单位复根……


Part2. 原根

对于素数 $p$,定义它的原根为 $G^0,G^1,G^2,…,G^{p-2}$ 使得这 $p-1$ 个数互不相同(其中 $G$ 是整数)。

再根据原根定义出:对于一个 $n-1$ 次多项式(也就是 $n$ 项的多项式), $g_n^k=(G^{\frac {p-1}n})^k\bmod p$($k=0\text ~(n-1)$)。其中需要保证 $\frac{p-1}n$ 是一个整数 ——
只要学过 FFT,应该会知道一定要先把多项式的项数改为一个 $2$ 的幂才能 FFT ——
NTT 也一样,所以 $n$ 很有可能是一个比较大的 $2$ 的幂,要使得 $\frac {p-1}n$ 是整数,那么 $p-1$ 就需要包含很多个因子 $2$,但是这样的 $p$ 并不多。这就造成了 NTT 的实际应用非常受限……

那么我们能否用 $g_n^k$ 来替换掉单位复根?这需要检验它在模意义下是否具有单位复根在 FFT 中的全部性质……

Part2/1. 与单位复根的相似性

在 FFT 中,我们总共用到了单位复根的这些性质:

  1. $n$ 个单位复根互不相同;
  2. $\omega_n^k=\omega_{2n}^{2k}$;
  3. $\omega_n^{k}=-\omega_n^{k+n/2}$;
  4. $\omega_n^a\times\omega_n^b=\omega_n^{a+b}$。

于是再来检验一下 $g_n^k$ 的性质:

  1. 根据原根的定义,$g_n^k$ 也互不相同;

  2. $g_{2n}^{2k}=(G^{\frac{p-1}{2n}})^{2k}=(G^{\frac{p-1}{n}})^k=g_n^k$;

  3. $g_n^{k+n/2}=(G^{\frac{p-1}n})^{k+n/2}=(G^{\frac{p-1}n})^k\cdot(G^{\frac{p-1}2})=g_n^k\cdot G^{\frac{p-1}2}$;

    因为 $(G^{\frac{p-1}2})^2=G^{p-1}=1$ (费马小定理)

    然后又根据原根的定义 $G^{\frac{p-1}2}\not=G^{p-1}$;

    所以得到 $G^{\frac{p-1}2}=-1$;

    所以 $g_n^{k+n/2}=-g_n^k$;

  4. $g_n^a\times g_n^b=(G^{\frac{p-1}n})^a\times(G^{\frac{p-1}n})^b=(G^{\frac{p-1}n})^{a+b}=g_n^{a+b}$;

开心的发现它具有单位复根的所有性质~

Part2/2. 常见的NTT模数及其原根

最为常见的 $998244353$ (建议背住),原根是 $3$;但是并不一定出现这个质数就是 NTT……

Copy 自 “not_exist” 的博客 CSDN

【点击展开/收起表格
质数 原根
3 2
5 2
17 3
97 5
193 5
257 3
7681 17
12289 11
40961 3
65537 3
786433 10
5767169 3
7340033 3
23068673 3
104857601 3
167772161 3
469762049 3
998244353 3
1004535809 3
2013265921 31
2281701377 3
3221225473 5
75161927681 3
77309411329 7
206158430209 22
2061584302081 7
2748779069441 3
6597069766657 5
39582418599937 5
79164837199873 5
263882790666241 7
1231453023109120 3
1337006139375610 3
3799912185593850 5
4222124650659840 19
7881299347898360 6
31525197391593400 3
180143985094819000 6
1945555039024050000 5
4179340454199820000 3

Part3. 源代码

打包成了一个 namespace……

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
namespace NTT{
typedef long long ll;
const ll MOD=998244353;
ll QPow(ll Va,ll Vb){ //一般的快速幂
ll Rr=1;
while(Vb){
if(Vb&1) Rr=Rr*Va%MOD;
Va=Va*Va%MOD;
Vb>>=1;
}
return Rr;
}
void Change(int *Cc,int len,int tmp){ //NTT
int logl=(int)round(log2(len));
for(int i=0;i<len;i++){
int Revi=0;
for(int j=0;j<logl;j++) Revi|=(((i>>j)&1)<<(logl-1-j));
if(Revi<i) swap(Cc[i],Cc[Revi]);
}
for(int s=2;s<=len;s<<=1){
int t=s/2;
ll Vper= tmp==1? QPow(3,(MOD-1)/s) : QPow(3,(MOD-1)-(MOD-1)/s);
for(int i=0;i<len;i+=s){
ll Vr=1;
for(int j=0;j<t;j++,Vr=Vr*Vper%MOD){
ll Fe=Cc[i+j],Fo=Cc[i+j+t];
Cc[i+j]=(Fe+Fo*Vr)%MOD;
Cc[i+j+t]=((Fe-Fo*Vr)%MOD+MOD)%MOD;
}
}
}
ll Einv=QPow(len,MOD-2);
if(tmp==-1)
for(int i=0;i<len;i++)
Cc[i]=Cc[i]*Einv%MOD;
}
void QPow(int *Ca,ll Vb,int len,int *Cr){ //多项式的快速幂
for(int i=0;i<len;i++) Cr[i]=0;
Cr[0]=1;
Change(Ca,len,1);Change(Cr,len,1);
while(Vb){
if(Vb&1)
for(int i=0;i<len;i++)
Cr[i]=1ll*Cr[i]*Ca[i]%MOD;
for(int i=0;i<len;i++)
Ca[i]=1ll*Ca[i]*Ca[i]%MOD;
Vb>>=1;
}
Change(Cr,len,-1);
}
}

The End

Thanks for reading!

Email: lucky_glass@foxmail.com,欢迎提问~

0%