初涉傅里叶变换

Fourier Transform

  • 傅里叶变换
  • 离散傅里叶变换
  • 快速傅里叶变换
    • 预备知识
    • FFT
  • 傅里叶逆变换
    • 预备知识
    • IDFT
  • 一些优化
    • 蝴蝶操作
    • 迭代实现
  • IDFT 另一种思路
  • FFT 的应用
    • A × B Problem 升级版


傅里叶变换

Fourier transform


  傅里叶变换(法语:Transformation de Fourier、英语:Fourier transform)是一种线性积分变换,用于信号在时域(或空域)和频域之间的变换,在物理学和工程学中有许多应用。

  摘自 wiki 百科

  简单来说, F T \mathrm{FT} FT 就是在某个域内将给定的信号映射到振幅,但此文意在简述信奥赛中, F F T \mathrm{FFT} FFT 的推导过程以及应用,故简单以此开头,不做更深入的解释。


离散傅里叶变换

Discrete Fourier transform


  离散傅里叶变换(Discrete Fourier transform,缩写为 DFT),是傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其 DTFT 的频域采样。

  摘自 wiki 百科

  定义序列 x = { x i ∣ i = 0 , 1 , 2 , ⋯   , N − 1 } x=\{x_i|i=0,1,2,\cdots,N-1\} x={xii=0,1,2,,N1},定义它的 D F T \mathrm{DFT} DFT x ^ = { x ^ i } \hat x=\{\hat x_i\} x^={x^i} x ^ k = ∑ n = 1 N − 1 e − i 2 π N n k x k k = 0 , 1 , 2 , ⋯   , N − 1 \hat x_k=\sum_{n=1}^{N-1}e^{-i\frac {2\pi}Nnk}x_k\quad k=0,1,2,\cdots,N-1 x^k=n=1N1eiN2πnkxkk=0,1,2,,N1  其中 i i i 为虚数单位,通常将 x ^ \hat x x^ 记为 x ^ = F x \hat x=\mathcal Fx x^=Fx y = y= y= D F T N \mathrm{DFT}_N DFTN ( x ) (x) (x)。离散傅里叶变换是可逆的线性变换,即有 F : C N → C N \mathcal F:\mathbb C^N\rightarrow\mathbb C^N F:CNCN

  定义序列 x ^ \hat x x^ 的离散傅里叶逆变换 ( I D F T \mathrm{IDFT} IDFT) 为 x k = 1 N ∑ n = 1 N − 1 e i 2 π N n k x ^ k k = 0 , 1 , 2 , ⋯   , N − 1 x_k=\frac 1N\sum_{n=1}^{N-1}e^{i\frac {2\pi}Nnk}\hat x_k\quad k=0,1,2,\cdots,N-1 xk=N1n=1N1eiN2πnkx^kk=0,1,2,,N1  通常将 x x x 记为 x = F − 1 x ^ x=\mathcal F^{-1}\hat x x=F1x^ x = x= x= I D F T N \mathrm{IDFT}_N IDFTN ( y ) (y) (y)

   D F T N \mathrm{DFT}_N DFTN I D F T N \mathrm{IDFT}_N IDFTN 前的归一化系数通常为 1 1 1 1 N \cfrac 1N N1,但有时我们也会统一使用 1 N \cfrac 1{\sqrt N} N 1


快速傅里叶变换

fast Fourier transform


  朴素的 D F T \mathrm{DFT} DFT,复杂度在 O ( n 2 ) O(n^2) O(n2),通过快速傅里叶变换 ( F F T \mathrm{FFT} FFT) 可以优化至 O ( n log ⁡ n ) O(n\log n) O(nlogn),不过嗯讲可能太干了,

  先介绍下作用,在信奥赛中, F F T \mathrm{FFT} FFT 通常用于加速多项式乘法。

  于是乎,我们先从一点基础内容开始。


预备知识

getting started


多项式基础


  在数学中,若干个关于 x x x 的单项式相加组成的代数式被称为多项式 f f f,其中最高次项的系数被称为多项式的度,记做 deg ⁡ f \deg f degf,于是 f f f 可以被系数表示为: f ( x ) = a 0 + a 1 x 1 + ⋯ + a deg ⁡ f x deg ⁡ f = ∑ i = 0 deg ⁡ f a i x i \begin{aligned}f(x) &= a_0 + a_1x^1+\cdots+a_{\deg f}x^{\deg f}\\&=\sum_{i=0}^{\deg f}a_ix^i\end{aligned} f(x)=a0+a1x1++adegfxdegf=i=0degfaixi


多项式系数表示法


  多项式系数表示法,就是用多项式的各个项的系数,组成一个系数序列来表示这个多项式,当然,没有要补 0 0 0,即 : : f ( x ) = ∑ i = 0 deg ⁡ f a i x i ⇔ { a 0 , a 1 , a 2 , ⋯   , a deg ⁡ f } f(x) =\sum_{i=0}^{\deg f}a_ix^i\quad\Leftrightarrow\quad\{a_0,a_1,a_2,\cdots,a_{\deg f}\} f(x)=i=0degfaixi{a0,a1,a2,,adegf}  以下,若无特殊说明,我们均将 f f f 的系数表示视为列向量 ( a 0 , a 1 , a 2 , ⋯   , a deg ⁡ f ) T (a_0,a_1,a_2,\cdots,a_{\deg f})^T (a0,a1,a2,,adegf)T


多项式点值表示法


  多项式 f f f 显然是一个函数,我们可以从图像上取出 deg ⁡ f + 1 \deg f + 1 degf+1 个不同的点值对来表示它,即 : : f ( x ) = ∑ i = 0 deg ⁡ f a i x i ⇔ { ( x 0 , y 0 ) , ( x 1 , y 0 ) , ( x 2 , y 2 ) , ⋯   , ( x deg ⁡ f , y deg ⁡ f ) } f(x) =\sum_{i=0}^{\deg f}a_ix^i\quad\Leftrightarrow\quad\{(x_0,y_0),(x_1,y_0),(x_2,y_2),\cdots,(x_{\deg f},y_{\deg f})\} f(x)=i=0degfaixi{(x0,y0),(x1,y0),(x2,y2),,(xdegf,ydegf)}   x k x_k xk 各不相同且都有 y k = f ( x k ) y_k = f(x_k) yk=f(xk),并且多项式点值表示是可逆的,即对于任意 n n n x k x_k xk 各不相同的点值对,都能唯一确定一个 n − 1 n-1 n1 次多项式,下面给出定理。


定理:多项式插值的唯一性


  对于任意 n + 1 n + 1 n+1 个点值对组成的集合 { ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n , y n ) } \{(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\} {(x0,y0),(x1,y1),,(xn,yn)},都可以唯一确定一个度数为 n n n 的多项式 f ( x ) = ∑ i = 0 n a i x i f(x) = \sum_{i=0}^na_ix^i f(x)=i=0naixi,满足 y k = f ( x k ) y_k=f(x_k) yk=f(xk)

  简单证明:

  显然有矩阵 [ 1 x 0 x 0 2 ⋯ x 0 n 1 x 1 x 1 2 ⋯ x 1 n ⋮ ⋮ ⋮ ⋱ ⋮ 1 x n x n 2 ⋯ x n n ] [ a 0 a 1 ⋮ a n ] = [ y 0 y 1 ⋮ y n ] \begin{bmatrix}1 & x_0 & x_0^2& \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix}y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix} 111x0x1xnx02x12xn2x0nx1nxnna0a1an=y0y1yn  左项为范德蒙矩阵,若对于任意 i , j i,j i,j 0 ≤ i , j ≤ n 0 \leq i,j \leq n 0i,jn 都有 x i ≠ x j x_i \neq x_j xi=xj,则左项可逆,左乘这个逆矩阵,即可得到 f ( x ) f(x) f(x) 系数表示,从而唯一的确定多项式 f ( x ) f(x) f(x)


多项式乘法


  给定多项式 A ( x ) = ∑ i = 0 deg ⁡ A a i x i A(x) = \sum_{i=0}^{\deg A}a_ix^i A(x)=i=0degAaixi B ( x ) = ∑ i = 0 deg ⁡ B b i x i B(x) = \sum_{i=0}^{\deg B}b_ix^i B(x)=i=0degBbixi

  定义 A × B A × B A×B 的积 C C C : : C ( x ) = A ( x ) × B ( x ) = ∑ i = 0 deg ⁡ A ∑ j = 0 deg ⁡ B a i b j x i + j \begin{aligned}C(x) &= A(x) × B(x)\\&=\sum_{i=0}^{\deg A}\sum_{j=0}^{\deg B}a_ib_jx^{i+j}\end{aligned} C(x)=A(x)×B(x)=i=0degAj=0degBaibjxi+j  对于度数分别为 n 、 m n、m nm 的多项式,它们的乘积度数为 n + m n + m n+m

  朴素的去计算多项式乘法,复杂度在 O ( n 2 ) O(n^2) O(n2),但现在我们考虑这样一个情况,若我们知道 : : A ( x ) = { ( x 0 , A ( x 0 ) ) , ( x 1 , A ( x 1 ) ) , ⋯   , ( x deg ⁡ A + deg ⁡ B , A ( x deg ⁡ A + deg ⁡ B ) ) } B ( x ) = { ( x 0 , B ( x 0 ) ) , ( x 1 , B ( x 1 ) ) , ⋯   , ( x deg ⁡ A + deg ⁡ B , B ( x deg ⁡ A + deg ⁡ B ) ) } A(x) =\{(x_0,A(x_0)),(x_1,A(x_1)),\cdots,(x_{\deg A + \deg B},A(x_{\deg A + \deg B}))\}\\B(x) =\{(x_0,B(x_0)),(x_1,B(x_1)),\cdots,(x_{\deg A + \deg B},B(x_{\deg A + \deg B}))\} A(x)={(x0,A(x0)),(x1,A(x1)),,(xdegA+degB,A(xdegA+degB))}B(x)={(x0,B(x0)),(x1,B(x1)),,(xdegA+degB,B(xdegA+degB))}  自然能在线性复杂度下计算出 : : C ( x ) = { ( x 0 , A ( x 0 ) B ( x 0 ) ) , ( x 1 , A ( x 1 ) B ( x 1 ) ) , ⋯   , ( x deg ⁡ A + deg ⁡ B , A ( x deg ⁡ A + deg ⁡ B ) ) B ( x deg ⁡ A + deg ⁡ B ) ) } C(x)=\{(x_0,A(x_0)B(x_0)),(x_1,A(x_1)B(x_1)),\cdots,(x_{\deg A + \deg B},A(x_{\deg A + \deg B}))B(x_{\deg A + \deg B}))\} C(x)={(x0,A(x0)B(x0)),(x1,A(x1)B(x1)),,(xdegA+degB,A(xdegA+degB))B(xdegA+degB))}   D F T \mathrm{DFT} DFT I D F T \mathrm{IDFT} IDFT 在多项式乘法中,扮演的就是多项式系数表示和点值表示之间,转换的这样一个角色,而 F F T \mathrm{FFT} FFT 可以将这个过程优化至 O ( n log ⁡ n ) O(n \log n) O(nlogn),故 F F T \mathrm{FFT} FFT 实现的多项式乘法,复杂度在 O ( n log ⁡ n ) O(n\log n) O(nlogn)


复数


  形如 z = a + b i z = a + bi z=a+bi a , b ∈ R a,b \in\mathbb R a,bR 的数被称为虚数, a a a 被称为实部, b b b 被称为虚部, i i i 被称为虚数单位。

  引入虚数相关的两个公式:

  虚数单位: i = − 1 i = \sqrt{-1} i=1
  欧拉函数: e i x = cos ⁡ ( x ) + i sin ⁡ ( x ) e^{ix} = \cos(x) + i\sin(x) eix=cos(x)+isin(x)


单位根


  数学上,n 次单位根是 n 次幂为 1 的复数。它们位于复平面的单位圆上,构成正多边形的顶点,但最多只可有两个顶点同时标在实数线上。

  摘自 wiki 百科

  定义 x n = 1 x^n= 1 xn=1 的复数根 x x x n n n 次单位根,根据欧拉函数,显然它的一个解为 ω n = e i 2 π n \omega_n = e^{i\frac{2\pi}{n}} ωn=ein2π,则它的解集可以表示为 ω n k \omega_n^k ωnk k = 0 , 1 , 2 , ⋯   , n − 1 k = 0,1,2,\cdots,n-1 k=0,1,2,,n1,故单位的 n n n 次根有 n n n 个。


单位根的性质


   ω n 0 = ω n n = 1 \omega_n^0 = \omega_n^n = 1 ωn0=ωnn=1
   ω 2 n n + k = − ω 2 n k \omega_{2n}^{n+k} = -\omega_{2n}^{k} ω2nn+k=ω2nk
   ω 2 n 2 k = ω n k \omega_{2n}^{2k}=\omega_n^k ω2n2k=ωnk

  简单推算一下式 2 2 2

   e i 2 ( n + k ) π 2 n = cos ⁡ ( π + k n π ) + i sin ⁡ ( π + k n π ) = − cos ⁡ ( k n π ) − i sin ⁡ ( k n π ) e^{i\frac{2(n+k)\pi}{2n}}=\cos(\pi +\frac kn \pi) + i\sin(\pi +\frac kn \pi) = -\cos(\frac kn \pi) - i\sin(\frac kn \pi) ei2n2(n+k)π=cos(π+nkπ)+isin(π+nkπ)=cos(nkπ)isin(nkπ)

   − e i 2 k π 2 n    = − cos ⁡ ( k n π ) − i sin ⁡ ( k n π ) -e^{i\frac{2k\pi}{2n}}\:\: =-\cos(\frac kn \pi) - i\sin(\frac kn \pi) ei2n2kπ=cos(nkπ)isin(nkπ)

  当然也可以从在坐标轴上的以原点为中心的单位单位圆上观察出上述三个公式,总之,无论推导方法,它们都是接下来 F F T \mathrm{FFT} FFT 中最为核心的一部分,所以熟悉它是掌握 F F T \mathrm{FFT} FFT 不可或缺的一部分。


FFT


  快速傅里叶变换的基本思想就是分治,我们取多项式 f ( x ) f(x) f(x) 上的 2 ⌈ log ⁡ 2 ( deg ⁡ f ) ⌉ 2^{\lceil\log_2(\deg f)\rceil} 2log2(degf) 个点,运用单位根的性质,分治的去计算 f ( x ) f(x) f(x) ω n k \omega_n^k ωnk k = 0 , 1 , 2 , ⋯   , n − 1 k = 0,1,2,\cdots,n-1 k=0,1,2,,n1 上的值,从而做到在 O ( n log ⁡ n ) O(n\log n) O(nlogn) 的复杂度下将 f f f 的系数表示转换为点值表示。

  具体地说,对于多项式 C ( x ) C(x) C(x),我们设 n = 2 ⌈ log ⁡ 2 ( deg ⁡ C ) ⌉ n = 2^{\lceil\log_2(\deg C)\rceil} n=2log2(degC) 则有 : : C ( x ) = a 0 + a 1 x 1 + ⋯ + a n − 1 x n − 1 C(x) = a_0 + a_1x^1+\cdots+a_{n-1}x^{n-1} C(x)=a0+a1x1++an1xn1  我们按奇偶次将 C C C 分为两组 : : C ( x ) = ( a 0 + a 1 x 2 + ⋯ + a n − 2 x n − 2 ) + ( a 1 x + a 3 x 3 + ⋯ + a n − 1 x n − 1 ) C(x) = (a_0 + a_1x^2+\cdots+a_{n-2}x^{n-2}) + (a_1x + a_3x^3+\cdots+a_{n-1}x^{n-1}) C(x)=(a0+a1x2++an2xn2)+(a1x+a3x3++an1xn1)  然后从右项中提出一个 x : x: x C ( x ) = ( a 0 + a 1 x 2 + ⋯ + a n − 2 x n − 2 ) + x ( a 1 + a 3 x 2 + ⋯ + a n − 1 x n − 2 ) C(x) = (a_0 + a_1x^2+\cdots+a_{n-2}x^{n-2}) + x(a_1 + a_3x^2+\cdots+a_{n-1}x^{n-2}) C(x)=(a0+a1x2++an2xn2)+x(a1+a3x2++an1xn2)  分别令 : : : A ( x ) = a 0 + a 1 x 2 + ⋯ + a n − 2 x n − 2 B ( x ) = a 1 + a 3 x 2 + ⋯ + a n − 1 x n − 2 \begin{aligned}A(x)&=a_0 + a_1x^2+\cdots+a_{n-2}x^{n-2}\\B(x)&=a_1 + a_3x^2+\cdots+a_{n-1}x^{n-2}\end{aligned} A(x)B(x)=a0+a1x2++an2xn2=a1+a3x2++an1xn2  此时 : : C ( x )    = A ( x 2 ) + x B ( x 2 ) C ( ω n x ) = A ( ω n 2 x ) + ω n x B ( ω n 2 x ) \begin{aligned}C(x)\ \ &=A(x^2)+xB(x^2)\\C(\omega_n^x)&=A(\omega_n^{2x})+\omega_n^xB(\omega_n^{2x})\end{aligned} C(x)  C(ωnx)=A(x2)+xB(x2)=A(ωn2x)+ωnxB(ωn2x)  考虑到 ω n 2 x = ω n / 2 x \omega_n^{2x} = \omega_{n/2}^x ωn2x=ωn/2x,则有 : : C ( ω n x ) = A ( ω n / 2 x ) + ω n x B ( ω n / 2 x ) C(\omega_n^x)=A(\omega_{n/2}^x)+\omega_n^xB(\omega_{n/2}^x) C(ωnx)=A(ωn/2x)+ωnxB(ωn/2x)  有 ω n x + n / 2 = − ω n x \omega_n^{x + n/2} = -\omega_n^x ωnx+n/2=ωnx,则有 : : C ( ω n x + n / 2 ) = A ( ω n / 2 x ) − ω n x B ( ω n / 2 x ) C(\omega_n^{x + n/2})=A(\omega_{n/2}^x)-\omega_n^xB(\omega_{n/2}^x) C(ωnx+n/2)=A(ωn/2x)ωnxB(ωn/2x)  此时,分治的去计算 A 、 B A、B AB ω n / 2 k \omega_{n/2}^k ωn/2k k = 0 , 1 , 2 , ⋯   , n 2 − 1 k = 0,1,2,\cdots,\frac n2-1 k=0,1,2,,2n1 上的值,即可在 T ( n ) = 2 T ( n 2 ) + O ( n ) T(n) =2T(\frac n2) +O(n) T(n)=2T(2n)+O(n) O ( n log ⁡ n ) O(n\log n) O(nlogn) 的复杂度下计算出 C C C ω n k \omega_n^k ωnk k = 0 , 1 , 2 , ⋯   , n − 1 k = 0,1,2,\cdots,n-1 k=0,1,2,,n1 上的值。


  在实现时,我们可以使用 C99 标准库提供的 来实现复数运行,但有一说一是真的不好用,CPP11 提供 也不大好用,凑合。

  简单的对比一下,以下代码片段中, 所使用的语法均出现在 上面。

  复数的声明 : :

double _Complex z = a + b * I;
std::complex<double> z(a, b);

  实部 | 虚部 : :

printf("z=%lf+%lfi", creal(z), cimag(z));
printf("z=%lf+%lfi", z.real(), z.imag());

  基本运算均重载了运算符,所以这里不做对比。

  好像半斤八两?


#include 

typedef std::complex<double> complex;

void FFT(complex *f, int n) {
    if (n == 1) return;
    int t = n >> 1;
    complex A[t], B[t];
    for (int i = 0; i < n; ++i)
        if (i & 1) B[i >> 1] = f[i];
        else A[i >> 1] = f[i];
    FFT(A, t), FFT(B, t);
    complex w(1), wn(cos(2 * M_PI / n), sin(2 * M_PI / n));
    for (int i = 0; i < t; ++i, w *= wn) {
        f[i] = A[i] + w * B[i];
        f[i + t] = A[i] - w * B[i];
    }
}

  如果不聊胜于无的话,我想还是挺值得一提的是 std::exp 重载了对复数的实现,我们可以直接使用 std::exp(complex(0, 2 * M_PI / n)) 来计算 ω n \omega_n ωn


傅里叶逆变换

inverse Fourier transform


预备知识

getting started


几何级数求和公式


  若一数列,从第一项 a ( a ≠ 0 ) a(a\neq0) a(a=0) 开始,以后它每一项都是它上一项乘以 r r r,即 S n = a + a r + a r 2 + ⋯ + a r n − 1 = ∑ i = 0 n − 1 a r i S_n=a+ar+ar^2+\cdots+ar^{n-1}=\sum_{i=0}^{n-1}ar^i Sn=a+ar+ar2++arn1=i=0n1ari  便称其为几何级数(或等比级数)。

  若 r ≠ 1 r \neq 1 r=1,则有 r S n − S n = a r n − a ( r − 1 ) S n = a ( r n − 1 ) S n = a ( r n − 1 ) r − 1 \begin{aligned}rS_n - S_n&= ar^n - a\\(r-1)S_n &= a(r^n - 1)\\S_n &= \frac{a(r^n - 1)}{r - 1}\end{aligned} rSnSn(r1)SnSn=arna=a(rn1)=r1a(rn1)  上式也被称为几何级数求和公式。


逆矩阵的性质


  ?好像忘了概括一下矩阵

  若矩阵 A A A 可逆,有 ( A − 1 ) T = ( A T ) − 1 (A^{-1})^T = (A^T)^{-1} (A1)T=(AT)1

  所以,若 A A A 是对称矩阵,则它的逆矩阵同样也是逆矩阵。


IDFT


  通过笔者在多项式插值的唯一性的简单证明中易得知, D F T \mathrm{DFT} DFT 在多项式 f ( x ) = { a 0 , a 1 , ⋯   , a deg ⁡ f } f(x) = \{a_0,a_1,\cdots,a_{\deg f}\} f(x)={a0,a1,,adegf} 上所做的,就是求出 [ 1 ω n 0 ω n 0 ⋯ ω n 0 1 ω n 1 ω n 2 ⋯ ω n n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ 1 ω n n − 1 ω n 2 ( n − 1 ) ⋯ ω n ( n − 1 ) 2 ] [ a 0 a 1 ⋮ a n − 1 ] = [ y 0 y 1 ⋮ y n − 1 ] \begin{bmatrix}1 & \omega_n^0 & \omega_n^0& \cdots & \omega_n^0\\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{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} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix}y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix} 111ωn0ωn1ωnn1ωn0ωn2ωn2(n1)ωn0ωnn1ωn(n1)2a0a1an1=y0y1yn1  其中 n = deg ⁡ f + 1 n = \deg f +1 n=degf+1,我们记 g = { y 0 , y 1 , ⋯   , y n − 1 } g = \{y_0,y_1,\cdots,y_{n-1}\} g={y0,y1,,yn1},容易想到求出 V ( ω n 0 , ω n 1 , ⋯   , ω n n − 1 ) − 1 V( \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1})^{-1} V(ωn0,ωn1,,ωnn1)1 左乘 g g g 傅里叶逆变换出 f f f,其中 V V V 表示范德蒙矩阵。

   V ( ω n 0 , ω n 1 , ⋯   , ω n n − 1 ) V( \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}) V(ωn0,ωn1,,ωnn1) 显然为对称矩阵,故 V ( ω n 0 , ω n 1 , ⋯   , ω n n − 1 ) − 1 = ( a i , j ) n × n V( \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1})^{-1} = (a_{i,j})_{n×n} V(ωn0,ωn1,,ωnn1)1=(ai,j)n×n 为对称矩阵, V ( ω n 0 , ω n 1 , ⋯   , ω n n − 1 ) ( a i , j ) n × n = E V( \omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1})(a_{i,j})_{n×n} = E V(ωn0,ωn1,,ωnn1)(ai,j)n×n=E,有 ∑ i = 1 n a x , i ω n y i = { 1 x = y 0 o t h e r w i s e \sum_{i=1}^na_{x,i}\omega_n^{yi}=\begin{cases}1&x=y\\0&\mathrm{otherwise}\end{cases} i=1nax,iωnyi={10x=yotherwise  根据几何级数求和公式,容易想到 ∑ i = 0 n − 1 ω n x i \sum_{i=0}^{n-1}\omega_n^{xi} i=0n1ωnxi,若 x ∈ Z x \in \mathbb Z xZ n ∤ x n \nmid x nx,则有 ( ω n x ) n − 1 ω n x − 1 = ( ω n n ) x − 1 ω n x − 1 = 0 \frac{(\omega_n^x)^n - 1}{\omega_n^x - 1} = \frac{(\omega_n^n)^x - 1}{\omega_n^x - 1} = 0 ωnx1(ωnx)n1=ωnx1(ωnn)x1=0 x = 0 x = 0 x=0,则有 ∑ i = 0 n − 1 ω n 0 = n \sum_{i=0}^{n-1}\omega_n^{0} = n i=0n1ωn0=n,故我们另 a i , j = ω n − ( i − 1 ) ( j − 1 ) a_{i,j} = \omega_n^{-(i-1)(j-1)} ai,j=ωn(i1)(j1),则有 ∑ i = 1 n ω n − x i ω n y i = { n x = y 0 o t h e r w i s e \sum_{i=1}^n\omega_n^{-xi}\omega_n^{yi}=\begin{cases}n&x=y\\0&\mathrm{otherwise}\end{cases} i=1nωnxiωnyi={n0x=yotherwise   0 ≤ x , y < n 0\leq x,y < n 0x,y<n,故 − n < x − y < n -nn<xy<n,上式从而成立。

  这里我们就能发现,对 g g g 做一遍 D F T \mathrm{DFT} DFT,变化时单位根指数取反,就能得到 n f nf nf

typedef std::complex<double> complex;

void FFT(complex *f, int n, int opt) {
    if (n == 1) return;
    int t = n >> 1;
    complex A[t], B[t];
    for (int i = 0; i < n; ++i)
        if (i & 1) B[i >> 1] = f[i];
        else A[i >> 1] = f[i];
    FFT(A, t, opt), FFT(B, t, opt);
    complex w(1), wn = std::exp(complex(0, opt * 2 * M_PI / n));
    for (int i = 0; i < t; ++i, w *= wn) {
        f[i] = A[i] + w * B[i];
        f[i + t] = A[i] - w * B[i];
    }
}

const int n = 8;

int main() {
    complex A[n]{1, 1, 4, 5, 1, 4};
    FFT(A, n, +1);
    FFT(A, n, -1);
    for (int i = 0; i < n; ++i) A[i] /= n;
    for (int i = 0; i < n; ++i)
        printf("%d", (int)(A[i].real() + 0.5));
    return 0;
}

一些优化


蝴蝶操作

bit-reversal permutation


  递归程序一种朴素的优化思路就是:递归转迭代,但要实现递归转迭代,第一件要解决的事情,就是弄清递归到 “叶节点” 时,当前的操作数与原系数序列的关系。

  蝴蝶操作可以在线性时间内,将原系数序列按分治操作最底层顺序重排,如当 n = 8 n=8 n=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 4 } { 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_4\}\{a_1,a_5\}\{a_3,a_7\} {a0,a1,a2,a3,a4,a5,a6,a7}{a0,a2,a4,a6}{a1,a3,a5,a7}{a0,a4}{a2,a4}{a1,a5}{a3,a7}  对 { a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 } \{a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7\} {a0,a1,a2,a3,a4,a5,a6,a7} 进行蝴蝶操作后,就能得到 { a 0 , a 4 , a 2 , a 4 , a 1 , a 5 , a 3 , a 7 } \{a_0,a_4,a_2,a_4,a_1,a_5,a_3,a_7\} {a0,a4,a2,a4,a1,a5,a3,a7}

  易得知 a i a_i ai A ( x ) 、 B ( x ) A(x)、B(x) A(x)B(x) 间,项数的奇偶性可能是不断变化的,比如 a 2 a_2 a2 为偶数项,则他在递归第一层的位置为 ⌈ 2 + 1 2 ⌉ \lceil\frac {2+1}2\rceil 22+1,变成了奇数项,而 a 3 a_3 a3 为奇数项,则他在递归第一层的位置为 ⌈ 3 + 1 2 ⌉ + n 2 \lceil\frac {3+1}2\rceil +\frac n2 23+1+2n,变成了偶数项。

  我们将 i i i a i a_i ai 最终的位置 i ′ i' i 按二进制拆分,容易发现,若 i i i 最低位为 1 1 1,则 i ′ i' i 的最高位为 1 1 1,否则为 0 0 0,然后对 ⌊ i 2 ⌋ \lfloor\frac i2\rfloor 2i(不要混淆了位置和序号)重复进行这个判断,最后 i ′ i' i 恰等于 i i i 在二进制表示下的反转。

  于是我们可以构造一个数列 r e v rev rev r e v i rev_i revi 表示 i i i 在二进制表示下反转的值,若我们顺序计算 r e v rev rev,容易发现计算 r e v i rev_i revi r e v i / 2 rev_{i/2} revi/2 已经被计算出来,若 i i i 为偶数。则 r e v i rev_i revi 等于 r e v i / 2 rev_{i/2} revi/2 右移一位,否则就在此基础上加上 n 2 \frac n2 2n,然后顺序交换元素就能在线性时间内完成蝴蝶操作。

#include 

typedef std::complex<double> complex;

void FFT(complex *f, int n, int opt) {
    if (n == 1) return;
    int t = n >> 1;
    FFT(f, t, opt), FFT(f + t, t, opt);
    complex w(1), wn = std::exp(complex(0, opt * 2 * M_PI / n));
    for (int i = 0; i < t; ++i, w *= wn) {
        complex tmp = w * f[i + t];
        f[i + t] = f[i] - tmp;
        f[i] += tmp;
    }
}

const int n = 8;

int rev[n];

int main() {
    for (int i = 0; i < n; ++i)
        rev[i] = rev[i >> 1] >> 1 | (i & 1) * (n >> 1);
    complex A[n]{1, 1, 4, 5, 1, 4};
    for (int i = 0; i < n; ++i)
        if (i < rev[i]) std::swap(A[i], A[rev[i]]);
    FFT(A, n, +1);
    for (int i = 0; i < n; ++i)
        if (i < rev[i]) std::swap(A[i], A[rev[i]]);
    FFT(A, n, -1);
    for (int i = 0; i < n; ++i) A[i] /= n;
    for (int i = 0; i < n; ++i)
        printf("%d", (int)(A[i].real() + 0.5));
    return 0;
}

迭代实现


  基操,没啥好说的。

#include 

typedef std::complex<double> complex;

const int N = 8;

int rev[N];

void FFT(complex *f, int lim, int opt) {
    for (int i = 0; i < lim; ++i)
        if (i < rev[i]) std::swap(f[i], f[rev[i]]);
    for (int n = 2; n <= lim; n <<= 1) {
        int t = n >> 1;
        complex wn = std::exp(complex(0, opt * 2 * M_PI / n));
        for (int i = 0; i < lim; i += n) {
            complex w(1);
            for (int j = 0; j < t; ++j, w *= wn) {
                complex tmp = w * f[i + j + t];
                f[i + j + t] = f[i + j] - tmp;
                f[i + j] = f[i + j] + tmp;
            }
        }
    }
    if (opt == -1)
        for (int i = 0; i < lim; ++i) f[i] /= lim;
}

int main() {
    for (int i = 0; i < N; ++i)
        rev[i] = rev[i >> 1] >> 1 | (i & 1) * (N >> 1);
    complex A[N]{1, 1, 4, 5, 1, 4};
    FFT(A, N, +1);
    FFT(A, N, -1);
    for (int i = 0; i < N; ++i)
        printf("%d", (int)(A[i].real() + 0.5));
    return 0;
}

IDFT 另一种思路


  设 f ( x ) = { a 0 , a 1 , ⋯   , a n − 1 } f(x) = \{a_0,a_1,\cdots,a_{n-1}\} f(x)={a0,a1,,an1} n = deg ⁡ f + 1 n = \deg f + 1 n=degf+1,构造多项式 g ( x ) = ∑ i = 0 n − 1 f ( ω n i ) x i g(x) = \sum_{i=0}^{n-1}f(\omega_n^i)x^i g(x)=i=0n1f(ωni)xi  将 ω n k \omega_{n}^k ωnk k = 0 , 1 , 2 , ⋯   , n − 1 k = 0,1,2,\cdots,n-1 k=0,1,2,,n1 依次代入式中,易知 g ( ω n k ) = ∑ i = 0 n − 1 ∑ j = 0 n − 1 a j ( ω n i ) j ( ω n i ) k = ∑ j = 0 n − 1 ∑ i = 0 n − 1 a j ( ω n j + k ) i = ∑ j = 0 n − 1 a j ∑ i = 0 n − 1 ( ω n j + k ) i \begin{aligned}g(\omega_n^k)&=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^i)^j(\omega_n^i)^k\\&=\sum_{j=0}^{n-1}\sum_{i=0}^{n-1}a_j(\omega_n^{j+k})^i\\&=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^{j+k})^i\end{aligned} g(ωnk)=i=0n1j=0n1aj(ωni)j(ωni)k=j=0n1i=0n1aj(ωnj+k)i=j=0n1aji=0n1(ωnj+k)i  根据几何级数求和公式,容易得知 g ( ω n k ) = n a n − k g(\omega_n^k) =na_{n - k} g(ωnk)=nank,计算过程大致与 I D F T \mathrm{IDFT} IDFT 相同,

  所以我们直接对 g g g 做一遍 D F T \mathrm{DFT} DFT,然后将第 2 2 2 n n n 个元素依次除 n n n 反转,整个操作的复杂度在 O ( n ) O(n) O(n),不过单位根取倒操作,单独拎出来复杂度在 O ( log ⁡ n ) O(\log n) O(logn)

  了解着玩就行了。

#include 

typedef std::complex<double> complex;

const int N = 8;

int rev[N];

void FFT(complex *f, int lim, bool inv) {
    for (int i = 0; i < lim; ++i)
        if (i < rev[i]) std::swap(f[i], f[rev[i]]);
    for (int n = 2; n <= lim; n <<= 1) {
        int t = n >> 1;
        complex wn = std::exp(complex(0, 2 * M_PI / n));
        for (int i = 0; i < lim; i += n) {
            complex w(1);
            for (int j = 0; j < t; ++j, w *= wn) {
                complex tmp = w * f[i + j + t];
                f[i + j + t] = f[i + j] - tmp;
                f[i + j] = f[i + j] + tmp;
            }
        }
    }
    if (inv) {
        for (int i = 1, j = lim - 1; i < j; ++i, --j)
            std::swap(f[i], f[j]), f[i] /= lim, f[j] /= lim;
        f[0] /= lim, f[lim >> 1] /= lim;
    }
}

int main() {
    for (int i = 0; i < N; ++i)
        rev[i] = rev[i >> 1] >> 1 | (i & 1) * (N >> 1);
    complex A[N]{1, 1, 4, 5, 1, 4};
    FFT(A, N, 0);
    FFT(A, N, 1);
    for (int i = 0; i < N; ++i)
        printf("%d", (int)(A[i].real() + 0.5));
    return 0;
}

FFT 的应用


   F F T \mathrm{FFT} FFT 的应用有很多,能力有限就不乱介绍了。


A × B Problem 升级版

P1919 【模板】A*B Problem 升级版


  给你两个正整数 a , b a,b a,b,求 a × b a \times b a×b 1 ≤ a , b ≤ 1 0 1000000 1≤a,b≤10^{ 1000000} 1a,b101000000


  设 A ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a ⌊ log ⁡ 10 a ⌋ x ⌊ log ⁡ 10 a ⌋ B ( x ) = b 0 + b 1 x + b 2 x 2 + ⋯ + b ⌊ log ⁡ 10 b ⌋ x ⌊ log ⁡ 10 b ⌋ \begin{aligned}A(x) &= a_0 + a_1x + a_2x^2 + \cdots + a_{\lfloor\log_{10}a\rfloor}x^{\lfloor\log_{10}a\rfloor}\\B(x) &= b_0 + b_1x + b_2x^2 + \cdots + b_{\lfloor\log_{10}b\rfloor}x^{\lfloor\log_{10}b\rfloor}\end{aligned} A(x)B(x)=a0+a1x+a2x2++alog10axlog10a=b0+b1x+b2x2++blog10bxlog10b  其中, a i 、 b j a_i、b_j aibj 分别为 a 、 b a、b ab 1 0 i 、 1 0 j 10^i、10^j 10i10j 位上的数,不难得知 a = A ( 10 ) 、 b = B ( 10 ) a = A(10)、b = B(10) a=A(10)b=B(10),而根据多项乘法的定义 C ( x ) = ∑ i = 0 deg ⁡ A ∑ j = 0 deg ⁡ B a i b j x i + j \displaystyle{C(x) = \sum_{i=0}^{\deg A}\sum_{j=0}^{\deg B}a_ib_jx^{i+j}} C(x)=i=0degAj=0degBaibjxi+j,不难得出 C ( 10 ) = ∑ i = 0 log ⁡ 10 a ∑ j = 0 log ⁡ 10 b a i b j 1 0 i + j = ∑ i = 0 log ⁡ 10 a a i 1 0 i ∑ j = 0 log ⁡ 10 b b j 1 0 j = a × b \begin{aligned}C(10) &= \sum_{i=0}^{\log_{10}a}\sum_{j=0}^{\log_{10}b}a_ib_j10^{i+j}\\&=\sum_{i=0}^{\log_{10}a}a_i10^i\sum_{j=0}^{\log_{10}b}b_j10^j\\&=a \times b\end{aligned} C(10)=i=0log10aj=0log10baibj10i+j=i=0log10aai10ij=0log10bbj10j=a×b  在笔者介绍多项式乘法时已经给出了一种 D F T \mathrm{DFT} DFT 做法,这里给出另一种更为高效的做法,考虑到复数 z = a + b i z = a + bi z=a+bi 的平方 z 2 = ( a + b i ) 2 = a 2 + 2 a b i + b 2 i 2 = ( a 2 − b 2 ) + 2 a b i \begin{aligned}z^2 &= (a + bi)^2\\&=a^2 +2abi +b^2i^2\\&=(a^2-b^2) +2abi\end{aligned} z2=(a+bi)2=a2+2abi+b2i2=(a2b2)+2abi  聪明的读者恐怕已经想到,我们构造一个多项式 G ( x ) = ∑ k = 0 i + j c i x i G(x) = \sum_{k=0}^{i+j}c_ix^i G(x)=k=0i+jcixi,将 a i a_i ai 放入到 c i c_i ci 的实部或虚部,将 b j b_j bj 放入到另一部,然后对 D F T N \mathrm{DFT}_N DFTN G G G 每一项系数都平方后再 I D F T \mathrm{IDFT} IDFT 回来,此时 c k c_k ck 的虚部除以二就是 C = A × B C = A × B C=A×B 真正的系数。

#include 
#include 
#include 

typedef std::complex<double> complex;

const int N = 1 << 21;

complex f[N];

char buf[N];

int rev[N];

void FFT(complex *f, int lim, int opt) {
    for (int i = 0; i < lim; ++i)
        if (i < rev[i]) std::swap(f[i], f[rev[i]]);
    for (int n = 2; n <= lim; n <<= 1) {
        int t = n >> 1;
        complex wn = std::exp(complex(0, opt * 2 * M_PI / n));
        for (int i = 0; i < lim; i += n) {
            complex w(1);
            for (int j = 0; j < t; ++j, w *= wn) {
                complex tmp = w * f[i + j + t];
                f[i + j + t] = f[i + j] - tmp;
                f[i + j] = f[i + j] + tmp;
            }
        }
    }
    if (opt == -1)
        for (int i = 0; i < lim; ++i) f[i] /= lim;
}

int n, m, len = 0, lim = 1;

int main() {
    scanf("%s", buf);
    n = strlen(buf);
    for (int i = 0, j = n; j; ++i)
        f[i] += complex(buf[--j] - '0', 0);
    scanf("%s", buf);
    m = strlen(buf);
    for (int i = 0, j = m; j; ++i)
        f[i] += complex(0, buf[--j] - '0');
    while (lim < n + m) lim <<= 1;
    for (int i = 0; i < lim; ++i)
        rev[i] = rev[i >> 1] >> 1 | (i & 1) * (lim >> 1);
    FFT(f, lim, +1);
    for (int i = 0; i < lim; ++i) f[i] *= f[i];
    FFT(f, lim, -1);
    for (int add = 0, now; len < n + m || add; ++len) {
        now = f[len].imag() / 2 + 0.5 + add;
        buf[len] = now % 10 + '0';
        add = now / 10;
    }
    if (buf[len - 1] == '0') --len;
    while (len) putchar(buf[--len]);
    return 0;
}

  优雅,永不过时。

你可能感兴趣的