NTT 入门_随笔
JueFan 一只绝帆

NTT 入门

牢记目标:对于一个给定的 m1m-1 次的多项式(m=2nm=2^n),我们要求出 mm 个点 ωmk\omega_{m}^{k} 代进去的点值,有了点值我们可以进行多项式乘法等内容。

把多项式奇偶分类,F(x)=F0(x2)+xF1(x2)F(x)=F_0(x^2)+xF_1(x^2)

由于我们递归下去要规模减半,考虑对点值分成两半,显然 ωmk+(m/2)=ωmk\omega_{m}^{k+(m/2)}=-\omega_{m}^k,所以偶次多项式的后一半点值等价于前一半点值,奇次多项式的后一半点值是前一半点值取反。

把点值还原成多项式(IDFT\rm IDFT)时,结论是指数上多一个负号,得到的结果应除以长度。

考虑构造一个多项式 G(x)=i=0m1F(ωmi)xiG(x)=\sum_{i=0}^{m-1}F(\omega_{m}^i)x^i,然后对其推式子。

G(ωk)=i=0m1ωmkij=0m1ajωmij=(i,j)ajωmi(jk)\begin{aligned}G(\omega^{-k})&=\sum_{i=0}^{m-1}\omega^{-ki}_m\sum_{j=0}^{m-1}a_j\omega_{m}^{ij} \\&=\sum_{(i,j)}a_j\omega_m^{i(j-k)}\end{aligned}

这里实质是一个单位根反演,我们发现 j=kj=k 时每个 ii 的贡献均为 11jkj\ne k 时所有 ii 的贡献和为 00,所以 G(ωk)=makG(\omega^{-k})=ma_k,所以求解 IDFT\rm IDFT 时只需把 ω\omega 换成 ω1\omega^{-1} 并在最后除掉总长度即可。

值得注意的是,如果我们并不把指数取反,而是代入原来的指数,那根据单位根的性质,其实我们相当于倒着代入了负指数,所以最后加一个 reverse 即可。

注意是从 11 开始的 reverse\rm reverse,因为 ω0=ω0\omega^0=\omega^{-0}

单位根换成原根后我们可以得到一个递归写法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
const int g=3,p=998244353;
ll ksm(ll a,ll n=p-2,ll res=1) {
for(;n;n>>=1,a=a*a%p) if(n&1) res=res*a%p;
return res;
}
void NTT(ll w,vec<ll> &a) {
if(sz(a)==1) return;
int n=sz(a);
vec<ll> b[2]{vec<ll>(n/2),vec<ll>(n/2)};
F(i,0,n-1) b[i&1][i>>1]=a[i];
NTT(w*w%p,b[0]);NTT(w*w%p,b[1]);
ll r=1;
F(i,0,n/2-1) {
a[i]=(b[0][i]+r*b[1][i])%p;
a[i+(n/2)]=(b[0][i]+(p-r)*b[1][i])%p;
(r*=w)%=p;
}
}
void DFT(vec<ll> &a) {NTT(ksm(g,(p-1)/sz(a)),a);}
void IDFT(vec<ll> &a) {
NTT(ksm(ksm(g,(p-1)/sz(a))),a);
ll iv=ksm(sz(a));
for(ll &x:a) (x*=iv)%=p;
}

为了减小常数,我们通常要将分治的过程循环化,我们发现我们先从上往下进行了奇偶分类,再自底向上合并。

如果我们要用循环代替递归,此时我们如何处理“奇偶分类”?

我们发现,奇偶分类恰恰是把最低位提到了最高位,所有奇偶分类都做完之后相当于二进制 reverse\rm reverse

这个过程叫蝴蝶变换,我们可以用 i2\frac i 2reverse\rm reverse 推出 iireverse\rm reverse,由此我们可以求出蝴蝶变换的位置:

1
F(i,0,(1<<n)-1) r[i]=r[i>>1]>>1|((i&1)<<n-1);

处理完奇偶分类后,我们就从底向上进行合并过程,这部分就直接抄过来即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void NTT(vec<ll> &a,bool type) {
int n=sz(a),l=__lg(n);
vec<int> r(n);
F(i,0,n-1) r[i]=r[i>>1]>>1|(i&1)<<l-1;
F(i,0,n-1) if(i<r[i]) swap(a[i],a[r[i]]);
for(int m=1;m<n;m*=2) {
ll w=ksm(type?g:ksm(g),(p-1)/(m*2));
for(int j=0;j<n;j+=m*2) {
ll r=1;
F(i,j,j+m-1) {
ll x=a[i],y=r*a[i+m]%p;
a[i]=(x+y)%p;a[i+m]=(x+p-y)%p;
(r*=w)%=p;
}
}
}
if(!type) {
ll iv=ksm(n);
F(i,0,n-1) (a[i]*=iv)%=p;
}
}
void DFT(vec<ll> &a) {NTT(a,1);}
void IDFT(vec<ll> &a) {NTT(a,0);}

这是一份符合大多数人刻板印象的 NTT 模板。

DIT 与 DIF

上文所述的 DFT\rm DFT 方式属于 DIT\rm DIT,按时域抽取,这里的按时域抽取指的是把多项式系数进行奇偶分类。

事实上,我们存在另一种 DFT\rm DFT 方式:DIF\rm DIF,按频域抽取,指的是把代入的点进行奇偶分类。

考虑 ω\omega 的奇偶性,也就是把多项式按照前后劈开:

DFT(a)k=i=0n1ωnkiai=i=0n/21(ai+(1)kai+n/2)ωnki\begin{aligned}\text{DFT}(a)_k&=\sum_{i=0}^{n-1}\omega_{n}^{ki}a_i \\&=\sum_{i=0}^{n/2-1}(a_i+(-1)^ka_{i+n/2})\omega_{n}^{ki}\end{aligned}

2k2\mid k 时:

DFT(a)k=i=0n/21(ai+ai+n/2)ωnki=i=0n/21(ai+ai+n/2)ωn/2(k/2)i=DFT(a 前后对应位置相加)k/2\begin{aligned}\text{DFT}(a)_k&=\sum_{i=0}^{n/2-1}(a_i+a_{i+n/2})\omega_{n}^{ki} \\&=\sum_{i=0}^{n/2-1}(a_i+a_{i+n/2})\omega_{n/2}^{(k/2)i} \\&=\text{DFT}(a\text{ 前后对应位置相加})_{k/2}\end{aligned}

2k2\nmid k 时略有不同:

DFT(a)k=i=0n/21(aiai+n/2)ωnki=i=0n/21(aiai+n/2)ωn/2iωn/2[k/2]i=DFT({(aiai+n/2)ωn/2i})k/2\begin{aligned}\text{DFT}(a)_k&=\sum_{i=0}^{n/2-1}(a_i-a_{i+n/2})\omega_{n}^{ki} \\&=\sum_{i=0}^{n/2-1}(a_i-a_{i+n/2})\omega_{n/2}^i\omega_{n/2}^{[k/2]i} \\&=\text{DFT}(\{(a_i-a_{i+n/2})\omega_{n/2}^i\})_{k/2}\end{aligned}

想象它的递归写法,我们在向下递归时进行值的变换(实质是乘上范德蒙德矩阵),在回溯时进行蝴蝶变换,刚好与 DIT\rm DIT 的过程是相反的。

那么,只要我们在 DFT\rm DFT 时使用 DIF\rm DIFIDFT\rm IDFT 时使用 DIT\rm DIT,就可以免去蝴蝶变换了。

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
void DIF(vec<ll> &a) {
int n=sz(a);
for(int m=n/2;m;m/=2) {
ll w=ksm(g,(p-1)/(m*2));
for(int j=0;j<n;j+=m*2) {
ll r=1;
F(i,j,j+m-1) {
ll x=a[i],y=a[i+m];
a[i]=(x+y)%p;a[i+m]=(x+p-y)*r%p;
(r*=w)%=p;
}
}
}
}
void DIT(vec<ll> &a) {
int n=sz(a);
for(int m=1;m<n;m*=2) {
ll w=ksm(ksm(g),(p-1)/(m*2));
for(int j=0;j<n;j+=m*2) {
ll r=1;
F(i,j,j+m-1) {
ll x=a[i],y=r*a[i+m]%p;
a[i]=(x+y)%p;a[i+m]=(x+p-y)%p;
(r*=w)%=p;
}
}
}
ll iv=ksm(n);
F(i,0,n-1) (a[i]*=iv)%=p;
}

看上去代码长,实际上二者逻辑相似,复制粘贴就很容易敲出。

比较简单的卡常方法:预处理单位根幂次,提交记录

虽然经过仔细卡常,10610^6 规模的数据还可以把速度提升至原来的 两倍,但具体内容已脱离有趣的理论范畴,且实战中并无用处,所以舍弃。

MTT

也就是任意模数 NTT\rm NTT

给一个使用 55FFT\rm FFT 的做法,较为舒服。

首先我们把系数变小,F(x)=A(x)+cB(x),G(x)=C(x)+cD(x)F(x)=A(x)+cB(x),G(x)=C(x)+cD(x),其中 c=Vc=\sqrt V,取 2152^{15} 即可。

然后多项式相乘展开,相当于我们要求 A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x)A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x),构造 T(x)=C(x)+iD(x)T(x)=C(x)+iD(x),则我们能在 A(x)T(x)A(x)T(x)B(x)T(x)B(x)T(x) 中找到所需的项。

只需三次 DFT\rm DFT,两次 IDFT\rm IDFT 即可。

由于我们要做很多次 FFT\rm FFT,所以预处理单位根就能很大程度上提高效率。

三次变两次优化

也就是只需要一遍 DFT\rm DFT 一遍 IDFT\rm IDFT 求出 fgf*g,它利用了 (f+gi)2=f2g2+2fgi(f+gi)^2=f^2-g^2+2fgi,所以我们只需对 (f+gi)(f+gi) 进行 DFT\rm DFT 即可。这玩意跑的比 NTT\rm NTT 快。

但是能 FFT\rm FFT 还无需 MTT\rm MTT 的场合太少了。

参考资料

再探 FFT – DIT 与 DIF,另种推导和优化

模板 P4245题解 - Prean

 评论
评论插件加载失败
正在加载评论插件
由 Hexo 驱动 & 主题 Keep
总字数 231.7k 访客数 访问量