本文已参与「新人创作礼」活动， 一起开启掘金创作之路。
快速傅里叶变换 (Fast Fourier Transform) ，即利用计算机计算离散傅里叶变换（DFT）的快速计算方法，简称 FFT，于1965年由 J.W.库利和 T.W.图基提出。我们在这里只讨论其在算法竞赛中的运用。
问题定义
对于两个多项式 f ( x ) = ∑ i = 0 n − 1 a i × x i f(x)=\sum_{i=0}^{n-1}a^i\times x^i f ( x ) = ∑ i = 0 n − 1 a i × x i g ( x ) = ∑ i = 0 m − 1 b i × x i g(x)=\sum_{i=0}^{m-1}b^i\times x^i g ( x ) = ∑ i = 0 m − 1 b i × x i f ( x ) × g ( x ) f(x)\times g(x) f ( x ) × g ( x )
问题分析
f ( x ) × g ( x ) = ( f × g ) ( x ) = ( ∑ i = 0 n − 1 a i × x i ) ( ∑ i = 0 m − 1 b i × x i ) = ∑ i = 0 n + m − 2 c i × x i \begin{aligned}
f(x)\times g(x)&=(f\times g)(x)\\
&=(\sum_{i=0}^{n-1}a^i\times x^i)(\sum_{i=0}^{m-1}b^i\times x^i)\\
&=\sum_{i=0}^{n+m-2}c_i\times x^i
\end{aligned} f ( x ) × g ( x ) = ( f × g ) ( x ) = ( i = 0 ∑ n − 1 a i × x i ) ( i = 0 ∑ m − 1 b i × x i ) = i = 0 ∑ n + m − 2 c i × x i 显然我们可以以 O ( n 2 ) O(n^2) O ( n 2 )
for (int i=0 ;i<=n;++i)
for (int j=0 ;j<=m;++j) c[i+j]=a[i]*b[j];
考虑优化。
刚才我们定义了 f ( x ) = ∑ i = 0 n − 1 a i × x i f(x)=\sum_{i=0}^{n-1}a^i\times x^i f ( x ) = ∑ i = 0 n − 1 a i × x i n n n 系数表示 形式，如 f ( x ) f(x) f ( x )
( a 0 , a 1 , a 2 , . . . , a n − 1 ) (a_0,a_1,a_2,...,a_{n-1}) ( a 0 , a 1 , a 2 , ... , a n − 1 ) 通过线性代数的相关知识我们知道，如果我们有该多项式所描述的线上的 n + 1 n+1 n + 1 点值表示 形式。所以 f ( x ) f(x) f ( x )
{ ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , . . . , ( x n − 1 , f ( x n − 1 ) ) } \{(x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1}))\} {( x 0 , f ( x 0 )) , ( x 1 , f ( x 1 )) , ... , ( x n − 1 , f ( x n − 1 ))} 如果我们分别在 f ( x ) f(x) f ( x ) g ( x ) g(x) g ( x ) n + m − 1 n+m-1 n + m − 1
{ ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , . . . , ( x n + m − 2 , f ( x n + m − 2 ) ) } { ( x 0 , g ( x 0 ) ) , ( x 1 , g ( x 1 ) ) , . . . , ( x n + m − 2 , g ( x n + m − 2 ) ) } \{(x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n+m-2},f(x_{n+m-2}))\}\\
\{(x_0,g(x_0)),(x_1,g(x_1)),...,(x_{n+m-2},g(x_{n+m-2}))\} {( x 0 , f ( x 0 )) , ( x 1 , f ( x 1 )) , ... , ( x n + m − 2 , f ( x n + m − 2 ))} {( x 0 , g ( x 0 )) , ( x 1 , g ( x 1 )) , ... , ( x n + m − 2 , g ( x n + m − 2 ))} 则 f ( x ) × g ( x ) f(x)\times g(x) f ( x ) × g ( x )
{ ( x 0 , f ( x 0 ) × g ( x 0 ) ) , ( x 1 , f ( x 1 ) × g ( x 1 ) ) , . . . , ( x n + m − 2 , f ( x n + m − 2 ) × g ( x n + m − 2 ) ) } \{(x_0,f(x_0)\times g(x_0)),(x_1,f(x_1)\times g(x_1)),...,(x_{n+m-2},f(x_{n+m-2})\times g(x_{n+m-2}))\} {( x 0 , f ( x 0 ) × g ( x 0 )) , ( x 1 , f ( x 1 ) × g ( x 1 )) , ... , ( x n + m − 2 , f ( x n + m − 2 ) × g ( x n + m − 2 ))} 如果我们可以做到快速转换多项式的系数表示和点值表示，我们就可以用 O ( n ) O(n) O ( n ) f ( x ) × g ( x ) f(x)\times g(x) f ( x ) × g ( x ) n + m − 1 n+m-1 n + m − 1 O ( n ) O(n) O ( n ) O ( n 2 ) O(n^2) O ( n 2 )
而 FFT 可以实现在 O ( n l o g n ) O(nlogn) O ( n l o g n )
DFT
DFT 的任务是把一个多项式从系数表示转化为点值表示。
我们考虑将多项式从系数表示转化为点值表示的时间瓶颈是什么。显然我们必须要最终得到 n n n O ( n ) O(n) O ( n ) 1 1 1 n n n
单位根
相信看到这里的大家都知道什么是复平面吧_(:3
在复平面中的单位圆上的点模长均为 1 1 1
则单位圆上所有的点可以表示为 ( c o s ( θ ) , i × s i n ( θ ) ) (cos(\theta),i\times sin(\theta)) ( cos ( θ ) , i × s in ( θ ))
n n n
由欧拉公式
e i × x = c o s ( x ) + i × s i n ( x ) e^{i\times x}=cos(x)+i\times sin(x) e i × x = cos ( x ) + i × s in ( x ) 当 x = π x=\pi x = π n , k ∈ N + n,k\in N+ n , k ∈ N +
e i π = − 1 e 2 i π = 1 e 2 k i π = 1 e 2 k i π n = 1 n e 2 k i π n = 1 n \begin{aligned}
e^{i\pi}&=-1\\
e^{2i\pi}&=1\\
e^{2ki\pi}&=1\\
\sqrt[n]{e^{2ki\pi}}&=\sqrt[n]{1}\\
e^\frac{2ki\pi}{n}&=\sqrt[n]{1}\\
\end{aligned} e iπ e 2 iπ e 2 kiπ n e 2 kiπ e n 2 kiπ = − 1 = 1 = 1 = n 1 = n 1 所以 n n n n n n k k k
e 2 k i π n = c o s ( 2 k π n ) + i × s i n ( 2 k π n ) e^\frac{2ki\pi}{n}=cos(\frac{2k\pi}{n})+i\times sin(\frac{2k\pi}{n}) e n 2 kiπ = cos ( n 2 kπ ) + i × s in ( n 2 kπ ) 由上式我们可以发现，单位根实际上就是把 ( 1 , 0 ) (1,0) ( 1 , 0 ) n n n n n n e 2 i π n = c o s ( 2 π n ) + i × s i n ( 2 π n ) e^\frac{2i\pi}{n}=cos(\frac{2\pi}{n})+i\times sin(\frac{2\pi}{n}) e n 2 iπ = cos ( n 2 π ) + i × s in ( n 2 π ) ω n \omega_n ω n
由三角函数的性质和基本数学运算法则可以得到单位根具有以下性质：
共有 n n n n n n k k k ω n k \omega_n^k ω n k
对所有的整数 p p p ω n k = ω p n p k \omega_n^k=\omega_{pn}^{pk} ω n k = ω p n p k
ω n k + n 2 = − ω n k \omega_n^{k+\frac{n}{2}}=-\omega_n^{k} ω n k + 2 n = − ω n k ω n k + n = ω n k \omega_n^{k+n}=\omega_n^{k} ω n k + n = ω n k ∑ i = 0 n − 1 ω n i = 0 \sum_{i=0}^{n-1}\omega_n^i=0 ∑ i = 0 n − 1 ω n i = 0
前四条都比较显然，我们证明一下第五条性质：
首先令
S = 1 + x + x 2 + x 3 + . . . + x n − 1 S=1+x+x^2+x^3+...+x^{n-1} S = 1 + x + x 2 + x 3 + ... + x n − 1 则
x S = x + x 2 + x 3 + . . . + x n − 1 + x n xS=x+x^2+x^3+...+x^{n-1}+x^n x S = x + x 2 + x 3 + ... + x n − 1 + x n 两式作差可以得到
x S − S = ( x + x 2 + x 3 + . . . + x n − 1 + x n ) − ( 1 + x + x 2 + x 3 + . . . + x n − 1 ) = x n − 1 = ( x − 1 ) × ( 1 + x + x 2 + x 3 + . . . + x n − 1 ) \begin{aligned}
xS-S&=(x+x^2+x^3+...+x^{n-1}+x^n)-(1+x+x^2+x^3+...+x^{n-1}
)\\
&=x^n-1\\
&=(x-1)\times(1+x+x^2+x^3+...+x^{n-1})
\end{aligned} x S − S = ( x + x 2 + x 3 + ... + x n − 1 + x n ) − ( 1 + x + x 2 + x 3 + ... + x n − 1 ) = x n − 1 = ( x − 1 ) × ( 1 + x + x 2 + x 3 + ... + x n − 1 ) 由 x n − 1 = ( x − 1 ) × ( 1 + x + x 2 + x 3 + . . . + x n − 1 ) x^n-1=(x-1)\times(1+x+x^2+x^3+...+x^{n-1}) x n − 1 = ( x − 1 ) × ( 1 + x + x 2 + x 3 + ... + x n − 1 )
当 x = 1 x=1 x = 1
否则，若 x n = 1 x^n=1 x n = 1
得证 ∑ i = 0 n − 1 ω n i = 0 \sum_{i=0}^{n-1}\omega_n^i=0 ∑ i = 0 n − 1 ω n i = 0
那么我们把单位根作为自变量有什么用处呢？
FFT
对于一个多项式 f ( x ) = ∑ i = 0 n − 1 a i × x i f(x)=\sum_{i=0}^{n-1}a^i\times x^i f ( x ) = ∑ i = 0 n − 1 a i × x i n n n 2 k 2^k 2 k 2 k 2^k 2 k
f ( x ) = ∑ i = 0 n − 1 a i × x i = a 0 × x 0 + a 1 × x 1 + . . . + a n × x n = ( a 0 × x 0 + a 2 × x 2 + . . . + a n − 2 × x n − 2 ) + ( a 1 × x 1 + a 3 × x 3 + . . . + a n − 1 × x n − 1 ) = ( 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 ) \begin{aligned}
f(x)&=\sum_{i=0}^{n-1}a^i\times x^i\\
&=a_0\times x^0+a_1\times x^1+...+a_n\times x^n\\
&=(a_0\times x^0+a_2\times x^2+...+a_{n-2}\times x^{n-2})+(a_1\times x^1+a_3\times x^3+...+a_{n-1}\times x^{n-1})\\
&=(a_0\times x^0+a_2\times x^2+...+a_{n-2}\times x^{n-2})+x \times(a_1\times x^0+a_3\times x^2+...+a_{n-1}\times x^{n-2})
\end{aligned} f ( x ) = i = 0 ∑ n − 1 a i × x i = a 0 × x 0 + a 1 × x 1 + ... + a n × x n = ( a 0 × x 0 + a 2 × x 2 + ... + a n − 2 × x n − 2 ) + ( a 1 × x 1 + a 3 × x 3 + ... + a n − 1 × x n − 1 ) = ( 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 ) 令
f 1 ( x ) = a 0 × x 0 + a 2 × x 1 + . . . + a n − 2 × x n 2 − 1 f 2 ( x ) = a 1 × x 0 + a 3 × x 1 + . . . + a n − 1 × x n 2 − 1 f_1(x)=a_0\times x^0+a_2\times x^1+...+a_{n-2}\times x^{\frac{n}{2}-1}\\
f_2(x)=a_1\times x^0+a_3\times x^1+...+a_{n-1}\times x^{\frac{n}{2}-1} f 1 ( x ) = a 0 × x 0 + a 2 × x 1 + ... + a n − 2 × x 2 n − 1 f 2 ( x ) = a 1 × x 0 + a 3 × x 1 + ... + a n − 1 × x 2 n − 1 则
f ( x ) = f 1 ( x 2 ) + x × f 2 ( x 2 ) f(x)=f_1(x^2)+x\times f_2(x^2) f ( x ) = f 1 ( x 2 ) + x × f 2 ( x 2 ) 将 n n n
f ( ω n k ) = f 1 ( ω n 2 k ) + ω n k × f 2 ( ω n 2 k ) = f 1 ( ω n 2 k ) + ω n k × f 2 ( ω n 2 k ) \begin{aligned}
f(\omega_n^k)&=f_1(\omega_n^{2k})+\omega_n^{k}\times f_2(\omega_n^{2k})\\
&=f_1(\omega_{\frac{n}{2}}^{k})+\omega_n^{k}\times f_2(\omega_{\frac{n}{2}}^{k})\\
\end{aligned} f ( ω n k ) = f 1 ( ω n 2 k ) + ω n k × f 2 ( ω n 2 k ) = f 1 ( ω 2 n k ) + ω n k × f 2 ( ω 2 n k ) 当 k ≥ n 2 k\ge\frac{n}{2} k ≥ 2 n k ′ = k − n 2 k'=k-\frac{n}{2} k ′ = k − 2 n
f ( ω n k ) = f 1 ( ω n 2 2 ( k ′ + n 2 ) ) + ω n k ′ + n 2 × f 2 ( ω n 2 2 ( k ′ + n 2 ) ) = f 1 ( ω n 2 k ′ ) − ω n k ′ × f 2 ( ω n 2 k ′ )
\begin{aligned}
f(\omega_n^k)&=f_1(\omega_{\frac{n}{2}}^{2(k'+\frac{n}{2})})+\omega_n^{k'+\frac{n}{2}}\times f_2(\omega_{\frac{n}{2}}^{2(k'+\frac{n}{2})})\\
&=f_1(\omega_{\frac{n}{2}}^{k'})-\omega_n^{k'}\times f_2(\omega_{\frac{n}{2}}^{k'})
\end{aligned} f ( ω n k ) = f 1 ( ω 2 n 2 ( k ′ + 2 n ) ) + ω n k ′ + 2 n × f 2 ( ω 2 n 2 ( k ′ + 2 n ) ) = f 1 ( ω 2 n k ′ ) − ω n k ′ × f 2 ( ω 2 n k ′ ) 如果我们递归地做下去，我们只要知道两个多项式 f 1 ( x ) f_1(x) f 1 ( x ) f 2 ( x ) f_2(x) f 2 ( x ) ( ω n 2 0 , ω n 2 1 , . . . , ω n 2 n 2 − 1 ) (\omega_\frac{n}{2}^0,\omega_\frac{n}{2}^1,...,\omega_\frac{n}{2}^{\frac{n}{2}-1}) ( ω 2 n 0 , ω 2 n 1 , ... , ω 2 n 2 n − 1 ) f ( x ) f(x) f ( x ) ( ω n 0 , ω n 1 , . . . , ω n n − 1 ) (\omega_n^0,\omega_n^1,...,\omega_n^{n-1}) ( ω n 0 , ω n 1 , ... , ω n n − 1 ) n = 1 n=1 n = 1 ω 1 = 1 \omega_1=1 ω 1 = 1
总体复杂度是 O ( n l o g n ) O(nlogn) O ( n l o g n )
IDFT
接下来我们考虑怎么把多项式从点值表示转化回系数表示。
如果我们把 DFT 表示为矩阵形式，就可以得到
[ ( ω n 0 ) 0 ( ω n 0 ) 1 ( ω n 0 ) 2 . . . ( ω n 0 ) n − 1 ( ω n 1 ) 0 ( ω n 1 ) 1 ( ω n 1 ) 2 . . . ( ω n 1 ) n − 1 ( ω n 2 ) 0 ( ω n 2 ) 1 ( ω n 2 ) 2 . . . ( ω n 2 ) n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ ( ω n n − 1 ) 0 ( ω n n − 1 ) 1 ( ω n n − 1 ) 2 . . . ( ω n n − 1 ) n − 1 ] [ a 0 a 1 a 2 ⋮ a n − 1 ] = [ y 0 y 1 y 2 ⋮ y n − 1 ] \begin{bmatrix}
(\omega_n^0)^0&(\omega_n^0)^1&(\omega_n^0)^2&...&(\omega_n^0)^{n-1}\\
(\omega_n^1)^0&(\omega_n^1)^1&(\omega_n^1)^2&...&(\omega_n^1)^{n-1}\\
(\omega_n^2)^0&(\omega_n^2)^1&(\omega_n^2)^2&...&(\omega_n^2)^{n-1}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&(\omega_n^{n-1})^2&...&(\omega_n^{n-1})^{n-1}\\
\end{bmatrix}
\begin{bmatrix}
a_0\\a_1\\a_2\\\vdots\\a_{n-1}
\end{bmatrix}
=
\begin{bmatrix}
y_0\\y_1\\y_2\\\vdots\\y_{n-1}
\end{bmatrix} ⎣ ⎡ ( ω n 0 ) 0 ( ω n 1 ) 0 ( ω n 2 ) 0 ⋮ ( ω n n − 1 ) 0 ( ω n 0 ) 1 ( ω n 1 ) 1 ( ω n 2 ) 1 ⋮ ( ω n n − 1 ) 1 ( ω n 0 ) 2 ( ω n 1 ) 2 ( ω n 2 ) 2 ⋮ ( ω n n − 1 ) 2 ... ... ... ⋱ ... ( ω n 0 ) n − 1 ( ω n 1 ) n − 1 ( ω n 2 ) n − 1 ⋮ ( ω n n − 1 ) n − 1 ⎦ ⎤ ⎣ ⎡ a 0 a 1 a 2 ⋮ a n − 1 ⎦ ⎤ = ⎣ ⎡ y 0 y 1 y 2 ⋮ y n − 1 ⎦ ⎤ 其中 a i a_i a i i i i y j y_j y j ω n j \omega_n^j ω n j
Ω = [ ( ω n 0 ) 0 ( ω n 0 ) 1 ( ω n 0 ) 2 . . . ( ω n 0 ) n − 1 ( ω n 1 ) 0 ( ω n 1 ) 1 ( ω n 1 ) 2 . . . ( ω n 1 ) n − 1 ( ω n 2 ) 0 ( ω n 2 ) 1 ( ω n 2 ) 2 . . . ( ω n 2 ) n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ ( ω n n − 1 ) 0 ( ω n n − 1 ) 1 ( ω n n − 1 ) 2 . . . ( ω n n − 1 ) n − 1 ] \Omega=
\begin{bmatrix}
(\omega_n^0)^0&(\omega_n^0)^1&(\omega_n^0)^2&...&(\omega_n^0)^{n-1}\\
(\omega_n^1)^0&(\omega_n^1)^1&(\omega_n^1)^2&...&(\omega_n^1)^{n-1}\\
(\omega_n^2)^0&(\omega_n^2)^1&(\omega_n^2)^2&...&(\omega_n^2)^{n-1}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&(\omega_n^{n-1})^2&...&(\omega_n^{n-1})^{n-1}\\
\end{bmatrix} Ω = ⎣ ⎡ ( ω n 0 ) 0 ( ω n 1 ) 0 ( ω n 2 ) 0 ⋮ ( ω n n − 1 ) 0 ( ω n 0 ) 1 ( ω n 1 ) 1 ( ω n 2 ) 1 ⋮ ( ω n n − 1 ) 1 ( ω n 0 ) 2 ( ω n 1 ) 2 ( ω n 2 ) 2 ⋮ ( ω n n − 1 ) 2 ... ... ... ⋱ ... ( ω n 0 ) n − 1 ( ω n 1 ) n − 1 ( ω n 2 ) n − 1 ⋮ ( ω n n − 1 ) n − 1 ⎦ ⎤ 则我们只需要找到 Ω \Omega Ω Ω − 1 \Omega^{-1} Ω − 1 Ω \Omega Ω
我们刚才证明了 n n n 0 0 0
[ ( ω n n − i ) 0 ( ω n n − i ) 1 ( ω n n − i ) 2 . . . ( ω n n − i ) n − 1 ] [ ( ω n j ) 0 ( ω n j ) 1 ( ω n j ) 2 ⋮ ( ω n j ) n − 1 ] = ∑ k = 0 n − 1 ( ω n n − i ) k ( ω n j ) k = ∑ k = 0 n − 1 ( ω n n − i × ω n j ) k = ∑ k = 0 n − 1 ( ω n n − i + j ) k = { n ( i = j ) 0 ( i ≠ j ) \begin{aligned}
\begin{bmatrix}
(\omega_n^{n-i})^0&(\omega_n^{n-i})^1&(\omega_n^{n-i})^2&...&(\omega_n^{n-i})^{n-1}
\end{bmatrix}
\begin{bmatrix}
(\omega_n^j)^0\\(\omega_n^j)^1\\(\omega_n^j)^2\\\vdots\\(\omega_n^j)^{n-1}
\end{bmatrix}
&=
\sum_{k=0}^{n-1}(\omega_n^{n-i})^k(\omega_n^j)^k\\
&=\sum_{k=0}^{n-1}(\omega_n^{n-i}\times \omega_n^j)^k\\
&=\sum_{k=0}^{n-1}(\omega_n^{n-i+j})^k\\
&=
\begin{cases}
n\ \ \ (i=j)\\
0\ \ \ (i\neq j)
\end{cases}
\end{aligned} [ ( ω n n − i ) 0 ( ω n n − i ) 1 ( ω n n − i ) 2 ... ( ω n n − i ) n − 1 ] ⎣ ⎡ ( ω n j ) 0 ( ω n j ) 1 ( ω n j ) 2 ⋮ ( ω n j ) n − 1 ⎦ ⎤ = k = 0 ∑ n − 1 ( ω n n − i ) k ( ω n j ) k = k = 0 ∑ n − 1 ( ω n n − i × ω n j ) k = k = 0 ∑ n − 1 ( ω n n − i + j ) k = { n ( i = j ) 0 ( i = j ) 所以我们可以得到 Ω \Omega Ω
Ω − 1 = 1 n [ ( ω n n ) 0 ( ω n n ) 1 ( ω n n ) 2 . . . ( ω n n ) n − 1 ( ω n n − 1 ) 0 ( ω n n − 1 ) 1 ( ω n n − 1 ) 2 . . . ( ω n n − 1 ) n − 1 ( ω n n − 2 ) 0 ( ω n n − 2 ) 1 ( ω n n − 2 ) 2 . . . ( ω n n − 2 ) n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ ( ω n 1 ) 0 ( ω n 1 ) 1 ( ω n 1 ) 2 . . . ( ω n 1 ) n − 1 ] = 1 n [ ( ω n 0 ‾ ) 0 ( ω n 0 ‾ ) 1 ( ω n 0 ‾ ) 2 . . . ( ω n 0 ‾ ) n − 1 ( ω n 1 ‾ ) 0 ( ω n 1 ‾ ) 1 ( ω n 1 ‾ ) 2 . . . ( ω n 1 ‾ ) n − 1 ( ω n 2 ‾ ) 0 ( ω n 2 ‾ ) 1 ( ω n 2 ‾ ) 2 . . . ( ω n 2 ‾ ) n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ ( ω n n − 1 ‾ ) 0 ( ω n n − 1 ‾ ) 1 ( ω n n − 1 ‾ ) 2 . . . ( ω n n − 1 ‾ ) n − 1 ] \begin{aligned}
\Omega^{-1}&=
\frac{1}{n}
\begin{bmatrix}
(\omega_n^n)^0&(\omega_n^n)^1&(\omega_n^n)^2&...&(\omega_n^n)^{n-1}\\
(\omega_n^{n-1})^0&(\omega_n^{n-1})^1&(\omega_n^{n-1})^2&...&(\omega_n^{n-1})^{n-1}\\
(\omega_n^{n-2})^0&(\omega_n^{n-2})^1&(\omega_n^{n-2})^2&...&(\omega_n^{n-2})^{n-1}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
(\omega_n^1)^0&(\omega_n^1)^1&(\omega_n^1)^2&...&(\omega_n^1)^{n-1}\\
\end{bmatrix}\\
&=
\frac{1}{n}
\begin{bmatrix}
(\overline{\omega_n^0})^0&(\overline{\omega_n^0})^1&(\overline{\omega_n^0})^2&...&(\overline{\omega_n^0})^{n-1}\\
(\overline{\omega_n^1})^0&(\overline{\omega_n^1})^1&(\overline{\omega_n^1})^2&...&(\overline{\omega_n^1})^{n-1}\\
(\overline{\omega_n^2})^0&(\overline{\omega_n^2})^1&(\overline{\omega_n^2})^2&...&(\overline{\omega_n^2})^{n-1}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
(\overline{\omega_n^{n-1}})^0&(\overline{\omega_n^{n-1}})^1&(\overline{\omega_n^{n-1}})^2&...&(\overline{\omega_n^{n-1}})^{n-1}\\
\end{bmatrix}
\end{aligned} Ω − 1 = n 1 ⎣ ⎡ ( ω n n ) 0 ( ω n n − 1 ) 0 ( ω n n − 2 ) 0 ⋮ ( ω n 1 ) 0 ( ω n n ) 1 ( ω n n − 1 ) 1 ( ω n n − 2 ) 1 ⋮ ( ω n 1 ) 1 ( ω n n ) 2 ( ω n n − 1 ) 2 ( ω n n − 2 ) 2 ⋮ ( ω n 1 ) 2 ... ... ... ⋱ ... ( ω n n ) n − 1 ( ω n n − 1 ) n − 1 ( ω n n − 2 ) n − 1 ⋮ ( ω n 1 ) n − 1 ⎦ ⎤ = n 1 ⎣ ⎡ ( ω n 0 ) 0 ( ω n 1 ) 0 ( ω n 2 ) 0 ⋮ ( ω n n − 1 ) 0 ( ω n 0 ) 1 ( ω n 1 ) 1 ( ω n 2 ) 1 ⋮ ( ω n n − 1 ) 1 ( ω n 0 ) 2 ( ω n 1 ) 2 ( ω n 2 ) 2 ⋮ ( ω n n − 1 ) 2 ... ... ... ⋱ ... ( ω n 0 ) n − 1 ( ω n 1 ) n − 1 ( ω n 2 ) n − 1 ⋮ ( ω n n − 1 ) n − 1 ⎦ ⎤ 所以将一个多项式在分治的过程中乘上的单位根变为其共轭复数，分治完的每一项除以 n n n
递归实现
#include <iostream>
#include <cstdio>
#include <cmath>
#include <vector>
using namespace std;
const double Pi=acos (-1.0 );
struct asdf {
long double x,y;
asdf operator + (const asdf b) const { return (asdf){x+b.x,y+b.y}; }
asdf operator - (const asdf b) const { return (asdf){x-b.x,y-b.y}; }
asdf operator * (const asdf b) const { return (asdf){x*b.x-y*b.y,x*b.y+y*b.x}; }
asdf operator * (const int b) const {return (asdf){x,y*b}; }
};
int N,M,l,tot=1 ;
void FFT (int len,vector<asdf> &a,int qwq)
if (len==1 ) return ;
vector<asdf> a1 (len>>1 ) ,a2 (len>>1 ) ;
for (int i=0 ;i<len;++i)
if (i%2 ) a2[i>>1 ]=a[i];
else a1[i>>1 ]=a[i];
FFT (len>>1 ,a1,qwq);
FFT (len>>1 ,a2,qwq);
asdf Wn={cos (2.0 *Pi/len),qwq*sin (2.0 *Pi/len)};
asdf w={1 ,0 };
for (int i=0 ;i<(len>>1 );i++)
{
a[i]=a1[i]+w*a2[i];
a[i+(len>>1 )]=a1[i]-w*a2[i];
w=w*Wn;
}
}
int main ()
int N,M;
scanf ("%d%d" ,&N,&M);
for (;tot<=N+M;l++) tot<<=1 ;
vector<asdf> a (tot+1 ) ,b (tot+1 ) ;
for (int i=0 ;i<=N;i++) scanf ("%Lf" ,&a[i].x);
for (int i=0 ;i<=M;i++) scanf ("%Lf" ,&b[i].x);
FFT (tot,a,1 ),FFT (tot,b,1 );
for (int i=0 ;i<=tot;i++) a[i]=a[i]*b[i];
FFT (tot,a,-1 );
for (int i=0 ;i<=N+M;i++)
{
printf ("%lld " ,(long long )(a[i].x/tot+0.5 ));
}
return 0 ;
}
非递归实现
但是这样计算常数很大，我们试图用非递归方法实现。
考虑一下奇偶分开时发生了什么。比如一个长度为 8 8 8
( a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 ) ( a 0 , a 2 , a 4 , a 6 ) ( a 1 , a 3 , a 5 , a 7 ) ( a 0 , a 4 ) ( a 2 , a 6 ) ( a 1 , a 5 ) ( a 3 , a 7 ) ( a 0 ) ( a 4 ) ( a 2 ) ( a 6 ) ( a 1 ) ( a 5 ) ( a 3 ) ( a 7 ) \
(a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7})\\
(a_{0}, a_{2}, a_{4}, a_{6})(a_{1}, a_{3}, a_{5}, a_{7})\\
(a_{0}, a_{4})(a_{2}, a_{6})(a_{1}, a_{5})(a_{3}, a_{7})\\
(a_{0})(a_{4})(a_{2})(a_{6})(a_{1})(a_{5})(a_{3})(a_{7}) ( a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 ) ( a 0 , a 2 , a 4 , a 6 ) ( a 1 , a 3 , a 5 , a 7 ) ( a 0 , a 4 ) ( a 2 , a 6 ) ( a 1 , a 5 ) ( a 3 , a 7 ) ( a 0 ) ( a 4 ) ( a 2 ) ( a 6 ) ( a 1 ) ( a 5 ) ( a 3 ) ( a 7 ) 如果我们可以直接将数组变成最后一步的样子，就可以自下而上递推计算了。
观察上述过程可以发现，第一次奇偶分列时，下标最低位是 0 0 0 2 2 2 0 0 0
所以系数 a i a_i a i i i i r [ i ] r[i] r [ i ] i i i t o t tot t o t [ 0 , t o t − 1 ] [0,tot-1] [ 0 , t o t − 1 ] t o t tot t o t 2 2 2 n = l o g 2 ( t o t ) n=log_2(tot) n = l o g 2 ( t o t )
假设我们当前枚举到了 i = ∑ j = 0 n − 1 s j × 2 j i=\sum_{j=0}^{n-1}s_j\times 2^j i = ∑ j = 0 n − 1 s j × 2 j
i = ( s n − 1 , s n − 2 , . . . , s 1 , s 0 ) 2 i=(s_{n-1},s_{n-2},...,s_1,s_0)_2 i = ( s n − 1 , s n − 2 , ... , s 1 , s 0 ) 2 将其右移一位得到
i > > 1 = ( 0 , s n − 1 , s n − 2 , . . . , s 1 ) 2 i>>1=(0,s_{n-1},s_{n-2},...,s_1)_2 i >> 1 = ( 0 , s n − 1 , s n − 2 , ... , s 1 ) 2 那么我们已经求得了将 i > > 1 i>>1 i >> 1
r [ i > > 1 ] = ( s 1 , s 2 , . . . , s n − 2 , s n − 1 , 0 ) 2 r[i>>1]=(s_1,s_2,...,s_{n-2},s_{n-1},0)_2 r [ i >> 1 ] = ( s 1 , s 2 , ... , s n − 2 , s n − 1 , 0 ) 2 所以我们只需要将 r [ i > > 1 ] r[i>>1] r [ i >> 1 ] s 0 < < ( n − 1 ) s_{0}<<(n-1) s 0 << ( n − 1 )
r [ i ] = ( s 0 , s 1 , s 2 , . . . , s n − 2 , s n − 1 ) 2 r[i]=(s_0,s_1,s_2,...,s_{n-2},s_{n-1})_2 r [ i ] = ( s 0 , s 1 , s 2 , ... , s n − 2 , s n − 1 ) 2
综上，转移方程为：
r [ i ] = ( r [ i > > 1 ] > > 1 ) ∣ ( ( ( i & 1 ) < < ( log 2 ( t o t ) − 1 ) ) r[i]=(r[i>>1]>>1)|(((i\&1)<<(\log_2(tot)-1)) r [ i ] = ( r [ i >> 1 ] >> 1 ) ∣ ((( i &1 ) << ( log 2 ( t o t ) − 1 )) 下面给出非递归实现的代码。
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const double Pi=acos (-1.0 );
struct asdf
{
double x,y;
}a[10000010 ],b[10000010 ];
asdf operator + (asdf a,asdf b)
{
return (asdf){a.x+b.x,a.y+b.y};
}
asdf operator - (asdf a,asdf b)
{
return (asdf){a.x-b.x,a.y-b.y};
}
asdf operator * (asdf a,asdf b)
{
return (asdf){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
int N,M,l,r[10000010 ],tot=1 ;
void FFT (asdf *A,int qwq)
for (int i=0 ;i<tot;i++)
if (i<r[i]) swap (A[i],A[r[i]]);
for (int mid=1 ;mid<tot;mid<<=1 )
{
asdf Wn=(asdf){cos (Pi/mid),qwq*sin (Pi/mid)};
for (int R=mid<<1 ,j=0 ;j<tot;j+=R)
{
asdf w=(asdf){1 ,0 };
for (int k=0 ;k<mid;k++,w=w*Wn)
{
asdf x=A[j+k],y=w*A[j+mid+k];
A[j+k]=x+y,A[j+mid+k]=x-y;
}
}
}
}
int main ()
int N,M;
scanf ("%d%d" ,&N,&M);
for (int i=0 ;i<=N;i++) scanf ("%lf" ,&a[i].x);
for (int i=0 ;i<=M;i++) scanf ("%lf" ,&b[i].x);
for (;tot<=N+M;l++) tot<<=1 ;
for (int i=0 ;i<tot;i++) r[i]=(r[i>>1 ]>>1 )|((i&1 )<<(l-1 ));
FFT (a,1 ),FFT (b,1 );
for (int i=0 ;i<=tot;i++) a[i]=a[i]*b[i];
FFT (a,-1 );
for (int i=0 ;i<=N+M;i++) printf ("%d " ,(int )(a[i].x/tot+0.5 ));
return 0 ;
}
预处理优化精度
队友用我之前写的上面那个板子在代码源提交大 WA 特 WA，发现精度损失主要出现在我们求第 k k k n n n
如果我们将其改成每次使用单位根时用用公式直接计算，需要多次计算三角函数值，很容易 TLE。
我们可以直接预估程序需要处理的项数最多的多项式有多少项，并将其补成 2 2 2 F F T N FFTN FFTN F F T N FFTN FFTN
对所有的整数 p p p ω n k = ω p n p k \omega_n^k=\omega_{pn}^{pk} ω n k = ω p n p k
快速且精度较高地求出我们需要的单位根。
代码如下：
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const double Pi=acos (-1.0 );
const int FFTN=262144 ;
struct asdf {
long double x,y;
asdf operator + (const asdf b) const { return (asdf){x+b.x,y+b.y}; }
asdf operator - (const asdf b) const { return (asdf){x-b.x,y-b.y}; }
asdf operator * (const asdf b) const { return (asdf){x*b.x-y*b.y,x*b.y+y*b.x}; }
asdf operator * (const int b) const {return (asdf){x,y*b}; }
}a[10000010 ],b[10000010 ],nw[10000010 ];
int N,M,l,r[10000010 ],tot=1 ;
void FFT (asdf *A,int qwq)
for (int i=0 ;i<tot;i++)
if (i<r[i]) swap (A[i],A[r[i]]);
for (int mid=1 ;mid<tot;mid<<=1 )
{
for (int R=mid<<1 ,j=0 ;j<tot;j+=R)
{
int t=FFTN/2 /mid;
for (int k=0 ;k<mid;k++)
{
asdf x=A[j+k],y=nw[t*k]*qwq*A[j+mid+k];
A[j+k]=x+y,A[j+mid+k]=x-y;
}
}
}
}
int main ()
for (int i=0 ;i<FFTN;++i)
nw[i]=(asdf){cosl (2 *Pi*i/FFTN),sinl (2 *Pi*i/FFTN)};
int N,M;
scanf ("%d%d" ,&N,&M);
for (int i=0 ;i<=N;i++) scanf ("%Lf" ,&a[i].x);
for (int i=0 ;i<=M;i++) scanf ("%Lf" ,&b[i].x);
for (;tot<=N+M;l++) tot<<=1 ;
for (int i=0 ;i<tot;i++) r[i]=(r[i>>1 ]>>1 )|((i&1 )<<(l-1 ));
FFT (a,1 ),FFT (b,1 );
for (int i=0 ;i<=tot;i++) a[i]=a[i]*b[i];
FFT (a,-1 );
for (int i=0 ;i<=N+M;i++) printf ("%lld " ,(long long )(a[i].x/tot+0.5 ));
return 0 ;
}
其他
虽然上面的代码已经比较优秀了，但是我们还可以对常数进行进一步的优化，日后遇到卡常的题目我再来进行补充。
