快速傅里叶变换(FFT) 基本概念 我们在计网中曾学过信号在时域和频域之间的转化需要傅里叶变换(Fourier Transform)
离散傅里叶变换(Discrete Fourier Transform，DFT),表示傅里叶变换在时域和频域上呈离散形式
而在算法竞赛上，DFT主要是用来解决多项式乘法等问题
FFT(快速傅里叶变换)即为优化版的DFT
前置知识 以处理多项式乘法为例
对于任意多项式,我们可以使用系数表示法,即按顺序排列其系数 ( a 0 , a 1 , a 2 , . . . , a n ) (a_0, a_1, a_2,..., a_n) ( a 0 , a 1 , a 2 , . . . , a n ) 对于一个n阶多项式,我们只需要n+1个点便可以确定(证明: 范德蒙特行列式) 所以要计算两个n阶多项式乘法f ( x ) , g ( x ) f(x),g(x) f ( x ) , g ( x )
如果我们知道多项式f(x)的点值,也知道g(x)的点值,那么我们就可以通过O(1)计算出f(x)和g(x)的乘积多项式的点值
系数转点值 离散傅里叶变换通过在复平面上找特定的点以将系数转为点值,公式如下:
X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ e i 2 π N k n X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{i\frac{2\pi}{N}kn}
X [ k ] = n = 0 ∑ N − 1 x [ n ] ⋅ e i N 2 π k n
但是DFT没有很好使用单位根的性质,依然复杂度是O ( n 2 ) O(n^2) O ( n 2 ) O ( n l o g n ) O(nlogn) O ( n l o g n )
对于多项式f ( x ) = a 0 + a 1 x 2 + a 2 x 3 + . . . + a n x n f(x)=a_0+a_1x^2+a_2x^3+...+a_{n}x^n f ( x ) = a 0 + a 1 x 2 + a 2 x 3 + . . . + a n x n
我们可以拆分奇偶
f ( x ) = ( a 0 + a 2 x 2 + . . . + a n − 1 x n − 1 ) + ( a 1 + a 3 x 2 + . . . + a n − 1 x n − 1 ) x f(x) = (a_0+a_2x^2+...+a_{n-1}x^{n-1}) + (a_1+a_3x^2+...+a_{n-1}x^{n-1})x
f ( x ) = ( a 0 + a 2 x 2 + . . . + a n − 1 x n − 1 ) + ( a 1 + a 3 x 2 + . . . + a n − 1 x n − 1 ) x
现在奇偶的系数都间隔为2,即[ 0 , 2 , 4 , . . . , n − 1 ] [0,2,4,...,n-1] [ 0 , 2 , 4 , . . . , n − 1 ] x 2 x^2 x 2 [ 0 , 1 , 2 , . . . , n − 1 2 ] [0,1,2,...,\frac{n-1}{2}] [ 0 , 1 , 2 , . . . , 2 n − 1 ]
所以:
f ( x ) = F 0 ( x 2 ) + F 1 ( x 2 ) x f(x) = F_0(x^2) + F_1(x^2)x
f ( x ) = F 0 ( x 2 ) + F 1 ( x 2 ) x
现在我们发现F 0 F_0 F 0 F 1 F_1 F 1
如何利用单位根的性质呢?
首先我们可以回忆一下单位根的性质
为了方便书写,我们用ω n k \omega^k_n ω n k
e i 2 π N k n = ω N n k e^{i\frac{2\pi}{N}kn} = \omega^{nk}_N
e i N 2 π k n = ω N n k
此时有:
ω n k = ω c n c k \omega^{k}_n = \omega^{ck}_{cn}
ω n k = ω c n c k
ω n k = − ω n k + n 2 \omega^{k}_n = - \omega^{k+\frac{n}{2}}_{n}
ω n k = − ω n k + 2 n
ω n 0 = ω n n = 1 \omega^0_n=\omega^n_n=1
ω n 0 = ω n n = 1
因此:
f ( x ) = F 0 ( x 2 ) + F 1 ( x 2 ) x f(x) = F_0(x^2) + F_1(x^2)x
f ( x ) = F 0 ( x 2 ) + F 1 ( x 2 ) x
f ( ω n k ) = F 0 ( ω n 2 k ) + F 1 ( ω n 2 k ) ω n k f(\omega^k_n) = F_0(\omega^{2k}_n) + F_1(\omega^{2k}_n)\omega^k_n
f ( ω n k ) = F 0 ( ω n 2 k ) + F 1 ( ω n 2 k ) ω n k
f ( ω n k ) = F 0 ( ω n 2 k ) + F 1 ( ω n 2 k ) ω n k f(\omega^k_n) = F_0(\omega^{k}_{\frac{n}{2}}) + F_1(\omega^{k}_{\frac{n}{2}})\omega^k_n
f ( ω n k ) = F 0 ( ω 2 n k ) + F 1 ( ω 2 n k ) ω n k
你可能会疑问但是这样递归下去依然没有减少计算量啊?
别忘了我们这里F 1 , F 0 F_1,F_0 F 1 , F 0 g ( x ) = g ( − x ) g(x)=g(-x) g ( x ) = g ( − x )
同时这里还有个最重要的性质没用到:
ω n k = − ω n k + n 2 \omega^{k}_n = - \omega^{k+\frac{n}{2}}_{n}
ω n k = − ω n k + 2 n
所以
f ( ω n k + n 2 ) = F 0 ( ω n 2 k ) − F 1 ( ω n 2 k ) ω n k f(\omega^{k+\frac{n}{2}}_n) = F_0(\omega^{k}_{\frac{n}{2}}) - F_1(\omega^{k}_{\frac{n}{2}})\omega^k_n
f ( ω n k + 2 n ) = F 0 ( ω 2 n k ) − F 1 ( ω 2 n k ) ω n k
所以现在我们只需要计算一半的点即可得到F 0 ( ω n 2 k ) , F 1 ( ω n 2 k ) F_0(\omega^{k}_{\frac{n}{2}}),F_1(\omega^{k}_{\frac{n}{2}}) F 0 ( ω 2 n k ) , F 1 ( ω 2 n k )
所以如果我们知道对于点[ ω n 2 1 , ω n 2 2 , ω n 2 3 , . . . , ω n 2 n ] [\omega^1_\frac{n}{2},\omega^2_\frac{n}{2},\omega^3_\frac{n}{2},...,\omega^n_\frac{n}{2}] [ ω 2 n 1 , ω 2 n 2 , ω 2 n 3 , . . . , ω 2 n n ] F 0 , F 1 F_0,F_1 F 0 , F 1 f ( ω n k ) f(\omega^k_n) f ( ω n k )
如何计算?还记得我们上面说的吗,递归即可
由此每次递归我们都会将子问题的计算规模减半,显然此时的复杂度即为O ( n l o g n ) O(nlogn) O ( n l o g n )
点值转系数 在得到点值表达式后,我们可以通过FFT的逆变换得到系数表达式
我们先从DFT的逆变换IDFT开始推导
还记得前面的DFT公式吗?
X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ e i 2 π N k n X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{i\frac{2\pi}{N}kn}
X [ k ] = n = 0 ∑ N − 1 x [ n ] ⋅ e i N 2 π k n
我们可以把DFT公式视为矩阵乘法:
[ 1 1 1 ⋯ 1 1 ω N 1 ω N 2 ⋯ ω N N − 1 1 ω N 2 ω N 4 ⋯ ω N 2 ( N − 1 ) ⋮ ⋮ ⋮ ⋱ ⋮ 1 ω N N − 1 ω N 2 ( N − 1 ) ⋯ ω N ( N − 1 ) 2 ] [ a 0 a 1 a 2 ⋮ a N − 1 ] = [ p 0 p 1 p 2 ⋮ p N − 1 ] \begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega_N^1 & \omega_N^2 & \cdots & \omega_N^{N-1} \\
1 & \omega_N^2 & \omega_N^4 & \cdots & \omega_N^{2(N-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega_N^{N-1} & \omega_N^{2(N-1)} & \cdots & \omega_N^{{(N-1)}^2}
\end{bmatrix}
\space \space
\begin{bmatrix}
a_0 \\
a_1 \\
a_2 \\
\vdots \\
a_{N-1}
\end{bmatrix}
{=}
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_{N-1}
\end{bmatrix}
⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1 1 1 ⋮ 1 1 ω N 1 ω N 2 ⋮ ω N N − 1 1 ω N 2 ω N 4 ⋮ ω N 2 ( N − 1 ) ⋯ ⋯ ⋯ ⋱ ⋯ 1 ω N N − 1 ω N 2 ( N − 1 ) ⋮ ω N ( N − 1 ) 2 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ ⎣ ⎢ ⎢ ⎢ ⎢ ⎡ a 0 a 1 a 2 ⋮ a N − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎡ p 0 p 1 p 2 ⋮ p N − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎤
其中[ a 0 a 1 a 2 ⋮ a N − 1 ] \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{N-1} \end{bmatrix} ⎣ ⎢ ⎢ ⎢ ⎢ ⎡ a 0 a 1 a 2 ⋮ a N − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎤ [ p 0 p 1 p 2 ⋮ p N − 1 ] \begin{bmatrix}p_0 \\ p_1 \\ p_2 \\ \vdots \\ p_{N-1} \end{bmatrix} ⎣ ⎢ ⎢ ⎢ ⎢ ⎡ p 0 p 1 p 2 ⋮ p N − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎤
所以如果我们要求系数表达式,只需要求原来的逆矩阵即可
同时我们有:
[ 1 1 1 ⋯ 1 1 ω N 1 ω N 2 ⋯ ω N N − 1 1 ω N 2 ω N 4 ⋯ ω N 2 ( N − 1 ) ⋮ ⋮ ⋮ ⋱ ⋮ 1 ω N N − 1 ω N 2 ( N − 1 ) ⋯ ω N ( N − 1 ) 2 ] T = 1 N [ 1 1 1 ⋯ 1 1 ( ω N 1 ) − 1 ( ω N 2 ) − 1 ⋯ ( ω N N − 1 ) − 1 1 ( ω N 2 ) − 1 ( ω N 4 ) − 1 ⋯ ( ω N 2 ( N − 1 ) ) − 1 ⋮ ⋮ ⋮ ⋱ ⋮ 1 ( ω N N − 1 ) − 1 ( ω N 2 ( N − 1 ) ) − 1 ⋯ ( ω N ( N − 1 ) 2 ) − 1 ] \begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega_N^1 & \omega_N^2 & \cdots & \omega_N^{N-1} \\
1 & \omega_N^2 & \omega_N^4 & \cdots & \omega_N^{2(N-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega_N^{N-1} & \omega_N^{2(N-1)} & \cdots & \omega_N^{{(N-1)}^2}
\end{bmatrix}^T{=}\frac{1}{N}
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & (\omega_N^1)^{-1} & (\omega_N^2)^{-1} & \cdots & (\omega_N^{N-1})^{-1} \\
1 & (\omega_N^2)^{-1} & (\omega_N^4)^{-1} & \cdots & (\omega_N^{2(N-1)})^{-1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & (\omega_N^{N-1})^{-1} & (\omega_N^{2(N-1)})^{-1} & \cdots & (\omega_N^{{(N-1)}^2})^{-1}
\end{bmatrix}
⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1 1 1 ⋮ 1 1 ω N 1 ω N 2 ⋮ ω N N − 1 1 ω N 2 ω N 4 ⋮ ω N 2 ( N − 1 ) ⋯ ⋯ ⋯ ⋱ ⋯ 1 ω N N − 1 ω N 2 ( N − 1 ) ⋮ ω N ( N − 1 ) 2 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ T = N 1 ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1 1 1 ⋮ 1 1 ( ω N 1 ) − 1 ( ω N 2 ) − 1 ⋮ ( ω N N − 1 ) − 1 1 ( ω N 2 ) − 1 ( ω N 4 ) − 1 ⋮ ( ω N 2 ( N − 1 ) ) − 1 ⋯ ⋯ ⋯ ⋱ ⋯ 1 ( ω N N − 1 ) − 1 ( ω N 2 ( N − 1 ) ) − 1 ⋮ ( ω N ( N − 1 ) 2 ) − 1 ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤
综上则有:
X [ k ] = ∑ n = 0 N − 1 x [ n ] ⋅ e − i 2 π N k n X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-i\frac{2\pi}{N}kn}
X [ k ] = n = 0 ∑ N − 1 x [ n ] ⋅ e − i N 2 π k n
所以将一个多项式在分治的过程中乘上的单位根变为其共轭复数，分治完的每一项除以n即可得到原多项式的每一项系数
代码实现详细步骤 有时间补
优化:蝴蝶变换 朴素FFT需要复制数组,并且是递归,常数较大内容到时候补,先蒯板子:
int rev[MAXN * 3 ];
const comp I (0 , 1 ) const double PI = acos (-1 );
void fft (comp F[], int N, int sgn = 1 )
int bit = __lg(N);
for (int i = 0 ; i < N; ++i)
{
rev[i] = (rev[i >> 1 ] >> 1 ) | ((i & 1 ) << (bit - 1 ));
if (i < rev[i])
swap (F[i], F[rev[i]]);
}
for (int l = 1 ; l < N; l <<= 1 )
{
comp step = exp (sgn * PI / l * I);
for (int i = 0 ; i < N; i += l * 2 )
{
comp cur (1 , 0 ) ;
for (int k = i; k < i + l; ++k)
{
comp g = F[k], h = F[k + l] * cur;
F[k] = g + h, F[k + l] = g - h;
cur *= step;
}
}
}
}
优化:三次变两次 还不会
更多 FFT在计算时会使用三角函数,因此无法避免运算时存在的精度误差
而NTT使用有限域上的单位根(即原根)代替复平面的单位根
因此NTT在计算时不会存在精度误差,且支持取模,并且常数似乎更小一点
或许有空我会更新NTT的笔记
参考文章 参考视频 
OPEN17的个人小站
