[學習筆記][省選算法]多項式乘法之 FFT & NTT 及卷積一、開頭二、預備知識三、多項式乘法四、多項式的點值表示及 DFT五、求 DFT六、求傅裡葉逆變換 IDFT （ DFT−1 DFT − 1 \text{DFT}^{-1} ）

一、開頭

（Place ：山東省神犇協會第 998244353 會議廳）

SD 神犇 998244353 号（會長）：所有神犇，洛谷黑題都切完了嗎？

全體神犇：切完了！

神犇 1004535809 号： 998244353 ，我昨天看到一個弱雞，洛谷名叫 xyz32768 ，他一道黑題都沒做，您認為應該怎樣呢？

神犇 998244353 号：這個問題您們考慮一下，應該用怎樣的方式把 xyz32768 D 一遍。時刻記住，我們協會的口号：見一弱雞， D 一弱雞！

神犇 167772161 号：要不，考他 FFT 吧。這是我們協會裡每個人都會的算法！

神犇 998244353 号：Good idea ！

（Place ： FFT 星球 DFT 國 IDFT 省 NTT 市）

神犇 167772161 号：哈哈，好久沒見，我們搞 OI 多年，能回憶起的比賽，至少也從三年前的 SDOI 2015 開始吧，那次比賽有一道十分模闆的題。題目名稱：序列統計。題面你自己去 BZOJ 或洛谷上去查，如果你 1h 内寫不出來我會對整個 FFT 星球說你菜！

xyz32768 ：BZOJ 3992 / 洛谷 P3321 ？

神犇 1004535809 号：對。這是一道紫題，但此大佬已經将之視為宇宙之水題了。

（蒟蒻 xyz32768 當然不會，看了看算法标簽）

xyz32768 ： 什麼？？？？？？ 快速傅裡葉變換,DFT,FFT ？？？？？？

神犇 998244353 号：哈哈哈哈哈，你真是菜啊， FFT 這麼水的算法都不會，再見蒟蒻！

神犇 998244353 号 & 神犇 1004535809 号 & 神犇 167772161 号（大聲喊）： xyz32768 菜， xyz32768 好菜， xyz32768 最菜！

二、預備知識

1、複數

下面的 i i ，除非作為 ∑∑ 求和的變量，其餘都表示虛數機關 −1−−−√ − 1 。

一個複數可以表示成 a+bi a + b i 的形式，也可以表示成 r(cosθ+isinθ) r ( cos ⁡ θ + i sin ⁡ θ ) 的形式。

上面 a a 為實部， bb 為虛部， r=a2+b2−−−−−−√ r = a 2 + b 2 為模長， θ=arctanba θ = arctan ⁡ b a 為輻角。

C++ 中，

#include <bits/stdc++.h>
using namespace std;
           

複數類為（ T 為實部和虛部的資料類型）：

complex<T>
           

定義一個複數（ a a 為實部， bb 為虛部）

complex <T> x(a, b)
           

實部

x.real()
           

一個複數 a+bi a + b i 對應複平面上一個點 (a,b) ( a , b ) ，模長 r r 為點到原點的距離，輻角 θθ 為終邊為 (0,0)−(a,b) ( 0 , 0 ) − ( a , b ) 的角（即該邊與 x 軸正半軸形成的角）。

2、機關根

ωn=cos2πn+isin2πn ω n = cos ⁡ 2 π n + i sin ⁡ 2 π n

那麼

ω0n,ω1n,ω2n,...,ωn−1n ω n 0 , ω n 1 , ω n 2 , . . . , ω n n − 1

稱為 n n 次機關根。

如：

ω04=1,ω14=i,ω24=−1,ω34=−iω40=1,ω41=i,ω42=−1,ω43=−i

n n 次機關根的性質 ：

（1） nn 次機關根所表示的點分布在複平面的機關圓上，并且把圓周等分成 n n 段。

（2） nn 次機關根的 n n 次方等于 11 。

（3） ω2i2n=ωin,ωn2n=−1,ωn+i2n=−ωi2n ω 2 n 2 i = ω n i , ω 2 n n = − 1 , ω 2 n n + i = − ω 2 n i

3、原根

在模 p p （為質數）意義下，如果存在數 gg ，滿足

g0,g1,g2,...,gp−2 g 0 , g 1 , g 2 , . . . , g p − 2

互不相同，那麼 g g 為模 pp 的原根。

原根都比較小，可以暴力枚舉。将 p−1 p − 1 分解質因數，對于一個 g g ，如果對于 p−1p−1 的任何一個質因子 x x ，都有

gp−1x≠1gp−1x≠1

那麼 g g 是模 pp 的一個原根。

三、多項式乘法

步入正題。一個 n n 次多項式 AA 和一個 n n 次多項式 BB 相乘，設 A A 的系數為 a0...na0...n （ xi x i 的項系數為 ai a i ）， B B 的系數為 b0...nb0...n ，那麼

abk=∑i+j=kaibj a b k = ∑ i + j = k a i b j

複雜度顯然是 O(n2) O ( n 2 ) 的。快速傅裡葉變換（ FFT ）則能做到 O(nlogn) O ( n log ⁡ n ) 的複雜度。

四、多項式的點值表示及 DFT

多項式的點值表示：選取 n+1 n + 1 個互不相同的值 x0...n x 0... n 代入 A A 對多項式進行求值，得到 n+1n+1 個多項式值，那麼這 n+1 n + 1 個自變量和這 n+1 n + 1 個多項式值就是多項式的點值表示。

舉例：一個四項式 A=3+2x−x2+x3 A = 3 + 2 x − x 2 + x 3 ，

選取自變量 x0=0,x1=−1,x2=1,x3=2 x 0 = 0 , x 1 = − 1 , x 2 = 1 , x 3 = 2 ，

求得 A(0)=3,A(−1)=−1,A(1)=5,A(2)=11 A ( 0 ) = 3 , A ( − 1 ) = − 1 , A ( 1 ) = 5 , A ( 2 ) = 11 。

這就是 A A 的一種點值表示。

一個多項式的點值表示有無窮多種，但是一個包含 n+1n+1 個自變量和多項式值的點值表示能唯一确定一個 n n 次多項式。

回到多項式乘法的内容，考慮到如果把 A,BA,B 變成點值表示，那麼就可以在 O(n) O ( n ) 時間裡得到 AB A B 的點值表示，然後再把 AB A B 轉化成系數表示。

于是，這個算法的瓶頸在于系數與點值表示下的轉換。

能否給這 n+1 n + 1 個變量取合适的值，使這個轉換能夠得到優化呢？

答案是肯定的。自變量取值就是 n+1 n + 1 次機關根。

将 ω0n,ω1n,...,ωn−1n ω n 0 , ω n 1 , . . . , ω n n − 1 代入一個 n n 項（注意，不是 nn 次）式 A A ，得到的點值表示，稱為 AA 的離散傅氏變換（ DFT DFT ）。

五、求 DFT

FFT 采用的是分治的思想。（這裡首先要将 n n 變成 22 的整數次幂）先将 A A 按照系數奇偶性分類：

A=(a0x0+a2x2+...+an−2xn−2)+(a1x1+a3x3+...+an−1xn−1)A=(a0x0+a2x2+...+an−2xn−2)+(a1x1+a3x3+...+an−1xn−1)

=(a0x0+a2x2+...+an−2xn−2)+x(a1x0+a3x2+...+an−1xn−2) = ( a 0 x 0 + a 2 x 2 + . . . + a n − 2 x n − 2 ) + x ( a 1 x 0 + a 3 x 2 + . . . + a n − 1 x n − 2 )

設：

B=a0x0+a2x1+...+an−2xn2−1 B = a 0 x 0 + a 2 x 1 + . . . + a n − 2 x n 2 − 1

C=a1x0+a3x1+...+an−1xn2−1 C = a 1 x 0 + a 3 x 1 + . . . + a n − 1 x n 2 − 1

假設已經求得了 B B 和 CC 的 DFT DFT ，那麼有：

A(ωin)=B((ωin)2)+ωinC((ωin)2) A ( ω n i ) = B ( ( ω n i ) 2 ) + ω n i C ( ( ω n i ) 2 )

=B(ω2in)+ωinC(ω2in) = B ( ω n 2 i ) + ω n i C ( ω n 2 i )

=B(ωin2)+ωinC(ωin2) = B ( ω n 2 i ) + ω n i C ( ω n 2 i )

是以，對于任何一個 i∈[0,n2) i ∈ [ 0 , n 2 ) ，

A(ωin)=B(ωin2)+ωinC(ωin2) A ( ω n i ) = B ( ω n 2 i ) + ω n i C ( ω n 2 i )

而由于 ωi+n2n2=ωin2 ω n 2 i + n 2 = ω n 2 i 且 ωi+n2n=−ωin ω n i + n 2 = − ω n i ，是以

A(ωi+n2n)=B(ωin2)−ωinC(ωin2) A ( ω n i + n 2 ) = B ( ω n 2 i ) − ω n i C ( ω n 2 i )

于是，我們實作了在 O(n) O ( n ) 的時間内，從 DFT(B) DFT ( B ) 和 DFT(C) DFT ( C ) 得到 DFT(A) DFT ( A ) 。

複雜度：設求一個 n n 項式的 DFTDFT 的複雜度為 T(n) T ( n ) ，那麼

T(n)=O(n)+2T(n2) T ( n ) = O ( n ) + 2 T ( n 2 )

根據主定理得到

T(n)=nlogn T ( n ) = n log ⁡ n

六、求傅裡葉逆變換 IDFT （ DFT−1 DFT − 1 ）

已知 A A 的 DFTDFT ：

A(ω0n),A(ω1n),...,A(ωn−1n) A ( ω n 0 ) , A ( ω n 1 ) , . . . , A ( ω n n − 1 )

求 A A 的系數表示。

解關于 a0,a1,...,an−1a0,a1,...,an−1 的線性方程組

⎧⎩⎨⎪⎪⎪⎪⎪⎪a0(ω0n)0+a1(ω0n)1+...+an−1(ω0n)n−1=A(ω0n)a0(ω1n)0+a1(ω1n)1+...+an−1(ω1n)n−1=A(ω1n)...a0(ωn−1n)0+a1(ωn−1n)1+...+an−1(ωn−1n)n−1=A(ωn−1n) { a 0 ( ω n 0 ) 0 + a 1 ( ω n 0 ) 1 + . . . + a n − 1 ( ω n 0 ) n − 1 = A ( ω n 0 ) a 0 ( ω n 1 ) 0 + a 1 ( ω n 1 ) 1 + . . . + a n − 1 ( ω n 1 ) n − 1 = A ( ω n 1 ) . . . a 0 ( ω n n − 1 ) 0 + a 1 ( ω n n − 1 ) 1 + . . . + a n − 1 ( ω n n − 1 ) n − 1 = A ( ω n n − 1 )

把方程組寫成矩陣形式：

⎡⎣⎢⎢⎢⎢(ω0n)0(ω1n)0...(ωn−1n)0(ω0n)1(ω1n)1...(ωn−1n)1............(ω0n)n−1(ω1n)n−1...(ωn−1n)n−1⎤⎦⎥⎥⎥⎥×⎡⎣⎢⎢⎢a0a1...an−1⎤⎦⎥⎥⎥=⎡⎣⎢⎢⎢⎢A(ω0n)A(ω1n)...A(ωn−1n)⎤⎦⎥⎥⎥⎥ [ ( ω n 0 ) 0 ( ω n 0 ) 1 . . . ( ω n 0 ) n − 1 ( ω n 1 ) 0 ( ω n 1 ) 1 . . . ( ω n 1 ) n − 1 . . . . . . . . . . . . ( ω n n − 1 ) 0 ( ω n n − 1 ) 1 . . . ( ω n n − 1 ) n − 1 ] × [ a 0 a 1 . . . a n − 1 ] = [ A ( ω n 0 ) A ( ω n 1 ) . . . A ( ω n n − 1 ) ]

設上面三個矩陣分别為 X,Y,Z X , Y , Z ，即

X×Y=Z X × Y = Z

考慮求 X−1 X − 1 。

設矩陣

R=⎡⎣⎢⎢⎢⎢(ω−0n)0(ω−1n)0...(ω−(n−1)n)0(ω−0n)1(ω−1n)1...(ω−(n−1)n)1............(ω−0n)n−1(ω−1n)n−1...(ω−(n−1)n)n−1⎤⎦⎥⎥⎥⎥ R = [ ( ω n − 0 ) 0 ( ω n − 0 ) 1 . . . ( ω n − 0 ) n − 1 ( ω n − 1 ) 0 ( ω n − 1 ) 1 . . . ( ω n − 1 ) n − 1 . . . . . . . . . . . . ( ω n − ( n − 1 ) ) 0 ( ω n − ( n − 1 ) ) 1 . . . ( ω n − ( n − 1 ) ) n − 1 ]

即 R R 第 ii 行第 j j 列（從 00 開始計數）的值為 ω−ijn ω n − i j 。

那麼 X×R X × R 第 i i 行第 ii 列的值為

∑j=0n−1ωijnω−ijn ∑ j = 0 n − 1 ω n i j ω n − i j

=∑j=0n−11=n = ∑ j = 0 n − 1 1 = n

第 i i 行第 jj 列（ i≠j i ≠ j ）的值為：

∑k=0n−1ωiknω−kjn=∑k=0n−1ω(i−j)kn=∑k=0n−1(ωi−j)k ∑ k = 0 n − 1 ω n i k ω n − k j = ∑ k = 0 n − 1 ω n ( i − j ) k = ∑ k = 0 n − 1 ( ω i − j ) k

=(ωi−jn)n−1ωi−jn−1=0 = ( ω n i − j ) n − 1 ω n i − j − 1 = 0

是以 X−1=1nR X − 1 = 1 n R 。

是以 Y=1nR×Z Y = 1 n R × Z

求 A A 的 DFT−1DFT−1 ，隻需要用 ω−in ω n − i 代替 ωin ω n i 做一遍 DFT DFT ，然後把求得的結果除以 n n 就能得到系數。

注：

ω−1n=ωn−1n=cos2πn−isin2πnωn−1=ωnn−1=cos⁡2πn−isin⁡2πn

七、FFT 的遞歸實作

cyx 為複數類， n 為項數， d 為機關根（如果是逆變換則為機關根的倒數），下為求 a0=y[le],a1=y[le+st],a2=[le+2st],... a 0 = y [ l e ] , a 1 = y [ l e + s t ] , a 2 = [ l e + 2 s t ] , . . . 的 DFT 。

void FFT(int n, cyx *y, int le, int st, cyx *d) {
    if (n == ) return; int m = n >> ;
    FFT(m, y, le, st << , d);
    FFT(m, y, le + st, st << , d); int i;
    for (i = ; i < m; i++) {
        int pos =  * st * i;
        tmp[i] = y[le + pos] + d[i * st] * y[le + pos + st];
        tmp[i + m] = y[le + pos] - d[i * st] * y[le + pos + st];
    }
    for (i = ; i < n; i++)
        y[le + i * st] = tmp[i];
}
           

八、FFT 的疊代實作

遞歸實作的 FFT 會遇到效率問題。我們發現， FFT 每次不是把區間直接分割成兩個子區間進行分治，是按照奇偶性進行分治。于是考慮如果 n=2k n = 2 k ，那麼設 rev[i] r e v [ i ] 表示 i i 在 kk 位二進制數下的各位數倒轉，如：

rev[(11001)2]=(10011)2 r e v [ ( 11001 ) 2 ] = ( 10011 ) 2

遞推公式（可以自行了解）：

rev[0]=0 r e v [ 0 ] = 0

rev[i]=(rev[i>>1]>>1) or ((i and 1)<<(k−1)) r e v [ i ] = ( r e v [ i >> 1 ] >> 1 )  or  ( ( i  and  1 ) << ( k − 1 ) )

如果在這個過程中，對 m=n2w m = n 2 w 個數求 DFT ，那麼參加求 DFT 的數的 下标有什麼特點？

————二進制表示的後 w w 位相等！

如果令 brev[i]=aibrev[i]=ai ，那麼這個特點就變成了二進制表示的前 w w 位相等，而二進制前 ww 位相等的數構成了一個連續區間！

這樣，從小區間合并到大區間，就是 FFT 的疊代實作了。

下面，如果 op = 1 則為 DFT ，否則 op = -1 ，為 IDFT 。

void FFT(int n, cyx *y, int op) {
    int i, j, k;
    for (i = ; i < n; i++) if (i < rev[i]) swap(y[i], y[rev[i]]);
    for (k = ; k < n; k <<= ) {
        cyx x(cos(pi / k), op * sin(pi / k));
        for (i = ; i < n; i += (k << )) {
            cyx w(, );
            for (j = ; j < k; j++) {
                cyx u = y[i + j], v = w * y[i + j + k];
                y[i + j] = u + v; y[i + j + k] = u - v;
                w = w * x;
            }
        }
    }
}
           

九、快速數論變換 NTT

FFT 涉及到複數運算，難免有精度問題。考慮如果題目要求在模意義下怎麼做。

思考原根是否有機關根的性質：

（1）根據費馬小定理， (gi)p−1≡1(modp) ( g i ) p − 1 ≡ 1 ( mod p ) 。

（2） gp−12≡−1(modp) g p − 1 2 ≡ − 1 ( mod p )

答案是肯定的。

是以當 p=a×2k+1 p = a × 2 k + 1 ， n=2x n = 2 x 且 x≤k x ≤ k 時，隻要将 ga×2k−x g a × 2 k − x 替換掉機關根，就能實作 NTT 。

下面是 ZZQ=p=1004535809 Z Z Q = p = 1004535809 時：（注： 3×334845270=1004535810 3 × 334845270 = 1004535810 ）

inline void NTT(const int &n, int *a, const int &op) {
    int i, j, k, r = op ==  ?  : ;
    for (i = ; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (k = ; k < n; k <<= ) {
        int x = qpow(r, (ZZQ - ) / (k << ), ZZQ);
        for (i = ; i < n; i += (k << )) {
            int w = ;
            for (j = ; j < k; j++) {
                int u = a[i + j], v = l * w * a[i + j + k] % ZZQ;
                a[i + j] = (u + v) % ZZQ; a[i + j + k] = (u - v + ZZQ) % ZZQ;
                w = l * w * x % ZZQ;
            }
        }
    }
}
           

十、NTT 的任意模數問題

如果 NTT 的模數不能分解成 a×2k+1 a × 2 k + 1 的形式，上面的做法就不可用。

解決辦法是選取幾個質數，使他們的積大于 n(p−1)2 n ( p − 1 ) 2 。

在這幾個質數的模意義下做 NTT ，然後用中國剩餘定理 CRT 合并後對 p p 取模。

十一、快速卷積

兩個定義域在自然數上的函數 f(n),g(n)f(n),g(n) ，他們的卷積為：

(f⨂g)(n)=∑i=0nf(i)g(n−i) ( f ⨂ g ) ( n ) = ∑ i = 0 n f ( i ) g ( n − i )

可以發現這是 f(n) f ( n ) 和 g(n) g ( n ) 對應兩個多項式（ xn x n 的系數分别為 f(n) f ( n ) 及 g(n) g ( n ) ）相乘之後的 n n 次項。用 FFT 可以求得 f⨂gf⨂g 的所有項值。

如果是求

∑i=0xf(i)g(i+k) ∑ i = 0 x f ( i ) g ( i + k )

那麼将 g g 翻轉，令 g′(i)=g(n−i)g′(i)=g(n−i) ，那麼就是 f⨂g′ f ⨂ g ′ ，仍然是一個卷積的形式。

十二、題目

1、[BZOJ3527][Zjoi2014]力：

https://www.lydsy.com/JudgeOnline/problem.php?id=3527

2、[BZOJ3992][SDOI2015]序列統計：

https://www.lydsy.com/JudgeOnline/problem.php?id=3992

3、[BZOJ4827][Hnoi2017]禮物：

https://www.lydsy.com/JudgeOnline/problem.php?id=4827

4、[BZJ3160]萬徑人蹤滅：

https://www.lydsy.com/JudgeOnline/problem.php?id=3160

5、[BZOJ4555][Tjoi2016&Heoi2016]求和：

https://www.lydsy.com/JudgeOnline/problem.php?id=4555

6、[BZOJ3456]城市規劃：

https://www.lydsy.com/JudgeOnline/problem.php?id=3456