快速傅里叶变换FFT

2023-10-06 14:46:10

快速傅里叶变换FFT(Fast Fourier Transformation)

本文主要讲述如何使用FFT来实现快速多项式乘法。

多项式的表示

系数表示

对于一个多项式

\[A(x)=\sum_{j=0}^{n-1} a_jx^j
\]

向量\(a=(a_0, a_1, ..., a_{n-1})\)为多项式的系数表示。

利用系数表示时，给出\(x_0\)，在\(O(n)\)的时间内可以求出\(A(x_0)\)的值，而且\(A(x)\)的加减法也可以在\(O(n)\)内求出，但乘法则需要\(O(n^2)\)的时间。

点值表示

一个多项式\(A(x)\)的点值表示为一个由\(n\)个点值对所组成的集合

\[\begin{Bmatrix}
(x_0, y_0), & (x_1, y_1), & ..., & (x_{n-1}, y_{n-1})
\end{Bmatrix}\]

使得对\(k=0, 1, 2, ..., n-1\), 所有\(x_k\)各不相同，且\(y_k=A(x_k)\).一个多项式可以有很多不同的点值表示。

次数界

如果一个多项式\(A(x)\)的最高次的非零系数为\(a_k\)，则称\(A(x)\)的次数为\(k\)，记\(degree(A)=k\)。任何严格大于一个多项式次数的整数都是该多项式的次数界。

定理：对于任意\(n\)个点值对组成的集合，其中\(x_k\)都不同，那么存在唯一的次数界为\(n\)的多项式\(A(x)\)，满足\(y_k=A(x_k), k=0, 1, 2, ..., n-1\)。

证明：由线性代数中的范德蒙行列式可证。

点值表示的优点

相对于系数表示，点值表示的加减乘法都可以在\(O(n)\)内算出(其中乘法因为当\(x\)一定时，\(y\)值相乘就是之后的\(x\)对应的\(y\)值。

因为一个多项式的点值表示可以有多种，所以可以选取一些较容易算的值，然后快速把一个多项式从系数表示转化为点值表示，运用点值表示相乘，然后再快速从点值表示转化为系数表示。

DFT

因此选择单位复数根作为\(x\)，其中\(n\)次单位复数根是满足\(\omega ^n=1\)的复数\(\omega\)，\(n\)次单位复数根有\(n\)个：\(\omega_k=e^{2\pi ik/n}, k=0, 1, 2, ..., n-1, \omega_n=e^{2\pi i/n}\)称为主\(n\)次单位根。

引理1（消去引理）：对于任何整数\(n, k \geq 0, d>0, \omega_{dn}^{dk}=\omega_{n}^{k}\)，特别地，对于任意偶数\(n>0\)，有\(\omega_{n}^{n/2}=\omega_2=-1\)

证明：\(\omega_{dn}^{dk}=(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k=\omega_n^k\)

引理2（折半引理）：如果\(n>0\)为偶数，那么\(n\)个\(n\)次单位复数根的平方的集合就是\(n/2\)个\(n/2\)次单位复数根的集合。

证明：根据消去引理，对任意非负整数\(k\)，有\((\omega_n^k)^2=\omega_{n/2}^k\)。注意，如果对于所有\(n\)词单位复数根进行平方，那么获得每个\(n/2\)次单位根正好\(2\)次，因为\((\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}\omega_n^n=\omega_n^{2k}=(\omega_n^k)^2\)

引理3（求和引理）：对任意整数\(n \geq 1\)和不能被\(n\)整除的非负整数\(k\)，有\(\sum_{j=0}^{n-1} (\omega_n^k)^j=0\)。

证明：

\[\sum_{j=0}^{n-1} (\omega_n^k)^j=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{0}{\omega_n^k-1}=0
\]

\(k\)不能被\(n\)整除保证分母不为\(0\)。

记

\[y_k=A(\omega_n^k)=\sum_{j=0}^{n-1} a_j\omega_n^{kj}
\]

则\(y=(y_0, y_1, ..., y_{n-1})\)就是系数向量\(a=(a_0, a_1, ..., a_{n-1})\)的离散傅里叶变换(DFT)，记为\(y=DFT_n(a)\).

FFT

(以下\(n\)默认为\(2\)的幂)

通过FFT，利用复数单位根的特殊性质，可以在\(O(n log n)\)的时间内算出\(DFT_n(a)\).FFT利用分治策略，采用\(A(x)\)中偶数下标的系数与奇数下标的系数，分别定义两个新的次数界为\(n/2\)的多项式\(A^{[0]}(x), A^{[1]}(x)\):

\[A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1}
\]

\[A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1}
\]

则\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)

因此\(A(x)\)在\(\omega_n^0, \omega_n^1, ..., \omega_n^{n-1}\)的值转换为：

求次数界为\(n/2\)的多项式\(A^{[0]}(x)\)和\(A^{[1]}(x)\)在点\((\omega_n^0)^2, (\omega_n^1)^2, ..., (\omega_n^{n-1})^2\)的值
根据\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)进行整合

根据折半引理，1中对点求值并不是\(n\)个不同的值，而是\(n/2\)个\(n/2\)次单位复数根。即把一个\(n\)个元素的\(DFT_n(a)\)计算划分为两个规模为\(n/2\)个元素的\(DFT_{n/2}\)计算。

当\(n=1\)时，\(y_0=a_0\omega_1^0=a_0\)

否则设\(y_k^{[0]}=A^{[0]}(\omega_{n/2}^k), y_k^{[1]}=A^{[1]}(\omega_{n/2}^k)\)，根据消去引理，有\(y_k^{[0]}=A^{[0]}(\omega_{n}^{2k}), y_k^{[1]}=A^{[1]}(\omega_{n}^{2k})\)

则

\[y_k=y_k^{[0]}+\omega_n^ky_k^{[1]}
\]

\[=A^{[0]}(\omega_{n}^{2k})+\omega_n^kA^{[1]}(\omega_{n}^{2k})
\]

\[=A(\omega_n^k)
\]

\[y_{k+(n/2)}=y_k^{[0]}-\omega_n^ky_k^{[1]}
\]

\[=y_k^{[0]}+\omega_n^{k+(n/2)}y_k^{[1]}
\]

\[=A^{[0]}(\omega_{n}^{2k})+\omega_n^{k+(n/2)}A^{[1]}(\omega_{n}^{2k})
\]

\[=A^{[0]}(\omega_{n}^{2k+n})+\omega_n^{k+(n/2)}A^{[1]}(\omega_{n}^{2k+n})
\]

\[=A(\omega_n^{k+(n/2)})
\]

逆运算

用矩阵表示\(y=V_na\)，则\(V_n\)在\((k, j)\)处元素为\(\omega_n^{kj}\)。求出\(V_n^{-1}\)，则可以求出\(a=DFT_n^{-1}(y)\)。

定理：\(V_n^{-1}\)的\((j, k)\)处元素为\(\omega_n^{-kj}/n\)。

证明：考虑\(V_n^{-1}V_n\)中(j, j')处的元素：

\[[V_n^{-1}V_n]_{jj'}=\sum_{k=0}^{n-1}(\omega_n^{-kj}/n)(\omega_n^{kj'})=\sum_{k=0}^{k(j'-j)}/n
\]

若\(j=j'\)，则此和为\(1\)；否则根据求和引理，此和为\(0\)，所以\([V_n^{-1}V_n]\)为单位矩阵。\(a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-kj}\)。

对比DFT：把\(a\)与\(y\)互换，用\(\omega_n^{-1}\)替换\(\omega_n\)，并将结果除以\(n\)。因此逆运算也可以在\(O(nlogn)\)内完成。

注意：

不够\(2\)的幂时要补零
做完FFT后顺序是乱的，应该按照编号的二进制的翻转对应的十进制进行排序。例如：\(n=8\)时

0 000 000 0

1 001 100 4

2 010 010 2

3 011 110 6

4 100 001 1

5 101 101 5

6 110 011 3

7 111 111 7

右边为对应编号

多项式乘法

#include <cstdio>

#include <cmath>

#include <algorithm>

#include <cstdlib>

#include <cstring>

#include <ctime>

#include <deque>

#include <queue>

#include <vector>

#include <map>

#include <complex>

using namespace std;

typedef complex<double> E;

const int maxn=int(1e6)+100;

const double PI=acos(-1);

int n, m;

int rev[maxn];

E a[maxn], b[maxn];

void init()

{

    scanf("%d%d", &n, &m);

    for (int i=0, tmp; i<=n; ++i)

    {

        scanf("%d", &tmp);

        a[i]=E(tmp, 0);

    }

    for (int i=0, tmp; i<=m; ++i)

    {

        scanf("%d", &tmp);

        b[i]=E(tmp, 0);

    }

    m=n+m;

    for (n=1; n<=m; n<<=1);

}

void FFT_init()

{

    for (int i=0; i<n; ++i)

        for (int j=0; 1<<j<n; ++j)

            rev[i]=(rev[i]<<1) | (i>>j & 1);

}

void FFT(E *a, int type)

{

    for (int i=0; i<=n; ++i)

        if (i<rev[i]) swap(a[i], a[rev[i]]);

    for (int i=2; i<=n; i<<=1)

        for (int j=0; j<n; j+=i)

        {

            E w(cos(2*PI/i), sin(type*2*PI/i)), wn(1, 0);

            for (int k=0; k<i>>1; ++k, wn*=w)

            {

                E tmp=a[j+k];

                a[j+k]=a[j+k]+wn*a[j+k+(i>>1)];

                a[j+k+(i>>1)]=tmp-wn*a[j+k+(i>>1)];

            }

        }

}

void solve()

{

    FFT(a, 1);

    FFT(b, 1);

    for (int i=0; i<n; ++i) a[i]=a[i]*b[i];

    FFT(a, -1);

    for (int i=0; i<=m; ++i)

        printf("%d ", int(a[i].real()/n+0.5));

}

int main()

{

    init();

    FFT_init();

    solve();

    return 0;

}

码农公寓