FFT 学习笔记

模板及讲解

简介

$DFT$:离散傅里叶变换
$IDFT$:离散傅里叶逆变换
$FFT$:快速傅里叶变换,计算多项式乘法 (卷积)
$FNTT/NTT$:快速傅里叶变换的优化版,优化常数及误差

FFT

多项式

系数表示法

设$A(x)$表示一个$n−1$次多项式

则$A(x)=\sum\limits_{i=0}^n a_ix^i$

利用这种方法计算多项式乘法复杂度为$O(n^2)$
(第一个多项式中每个系数都需要与第二个多项式的每个系数相乘)

点值表示法

将$n$互不相同的$x$带入多项式,会得到$n$个不同的取值$y$

则该多项式被这$n$个点$(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)$唯一确定

其中$y_i=\sum\limits_{j=0}^{n-1} a_jx_i^j$

利用这种方法计算多项式乘法复杂度为$O(n^2)$

点值表示法用处

我们可以把系数表示转化成点值表示点值表示下的多项式就是选相同的$x_i$,把对应的$y_i$相乘
两个长度为$n$的多项式相乘得到的是长度为$2n−1$的多项式,需要$2n−1$个点值才能唯一表示,因此一开始两个多项式也要选$2n−1$个点表示。

Markdown

随意选$O(n)$个点,计算$y$是$O(n)$的,总时间还是$O(n^2)$的,所以我们要介绍后面的利器,单位根

复数

定义

设$a,b$为实数,$i^2=−1$,形如$a+bi$的数叫复数,其中$i$被称为虚数单位

在复平面中,$x$代表实数,$y$轴(除原点外的点)代表虚数,从原点$(0,0)$到$(a,b)$的向量表示复数$a+bi$
模长:从原点$(0,0)$到点$(a,b)$的距离,即$\sqrt{a^2+b^2}$
幅角:假设以逆时针正方向,从$x$轴正半轴到已知向量的转角的有向角叫做幅角

运算法则

加法:$(a,b)+(x,y)=(a+x,b+y)$
乘法:$(a,b) \times (x,y)=(ac-bd, bc+ad)$ (可从定义式中推出)

单位根

定义

$n$次单位根(记为$\omega$)
它是$n$个复数的集合,每一个的$n$次方都等于$1$。其中的一个是$e^{\frac {2\pi i} n}$,记为$\omega_n$。

欧拉公式

$$
e^{ix}=\cos x+(\sin x)i
$$

$x$就是这个复数向量的旋转角

易得
$$
\omega_{n}^{k}=\cos k \frac{2\pi}{n}+i\sin k\frac{2\pi}{n}
$$

几何意义

复平面上,以原点为圆心,$1$为半径作圆,所得的圆叫单位圆。以圆点为起点,圆的$n$等分点为终点,做$n$个向量,设幅角为正最小的向量对应的复数为$\omega_n$,称为$n$次单位根

下面是算导的图片,这$n$个复数向量在单位圆上呈放射状
Markdown

性质

1、消去引理 (从定义式可推,可从图中观察):
$$
\omega _{2n}^{2k}=\omega _{n}^{k}
$$

2、(从定义式可推,可从图中观察)
$$
\omega _{n}^{k+\frac{n}{2}}=-\omega _{n}^{k}
$$

3、(从定义式可推,可从图中观察)
$$
\omega_n^0=\omega_n^n=1
$$

4、(可从图中观察)
$$
\omega^{\frac n 2}_n=-1
$$

$DFT$和$IDFT$

$DFT$

我们把$\omega ^k_n(k \in [n−1])$分别带入多项式求点值

假设$n$为偶数,那么我们可以把它的奇偶项分开,用两个多项式表示它
$$
A^{[0]}(x)=a_0+a_2x+a_4x^2+…+a_{n−2}x^{\frac n2−1} \\
A^{[1]}(x)=a_1+a_3x+a_5x^2+…+a_{n−1}x^{\frac n2−1}
$$
则 $A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$

带入单位根$((k < \frac n 2))$
$$
A(\omega^k_n)=A^{[0]}(\omega^{2k}_n)+\omega^k_nA^{[1]}(\omega^{2k}_n) \\
=A^{[0]}(\omega^k_{\frac n 2})+\omega^k_nA^{[1]}(\omega^k_{\frac n 2})(k < \frac n 2)
$$

$((k \geq \frac n 2))$时
$$
A(\omega^{\frac n 2+k}_n)=A^{[0]}(\omega^{n+2k}_n)+\omega^{\frac n 2+k}_nA^{[1]}(\omega^{n+2k}_n)\\
=A^{[0]}(\omega^{2k}_n)+\omega^{\frac n 2}_n\omega^k_nA^{[1]}(\omega^{2k}_n) \\
=A^{[0]}(\omega^k_{\frac n 2})-\omega^k_nA^{[1]}(\omega^k_{\frac n 2})
$$

可以由算法写出$DFT$的复杂度$T(n)=2T(\frac n 2)+O(n)=O(n\log n)$

$IDFT$

还原成系数表示的过程叫做$IDFT​$。

只需要记住,$IDFT$的$\omega_n$是$e^{-\frac{2\pi i}n}$,最后的结果除以$n$,其它过程同$DFT$,可以写在一个函数里。

蝴蝶操作和$Rader$排序

在操作过程中,取出了a0[k]a1[k]的值,通过计算又求出了a[k]a[k+n]的值。我们把这样的一次运算叫做「蝴蝶操作」。

以$n=8$为例,看看递归过程的结构
Markdown

我们完全可以从底向上递推求解
剩下的问题就是把初始的数组变成最后一层的样子了。
观察最后序列的编号的二进制表示000,100,010,110,001,101,011,111与原来000,001,010,011,100,101,110,111相比,每个位置上的二进制位都反过来了。这样的变化叫做$Rader$排序

我们可以$O(n)$将$Rader$排序的映射关系求出。设$iRader$排序后的数为$r_i$,我们可以通过$r_{\frac i2}$递推求出

这个递推式相当于偶数位$i​$的值就是$\frac i2​$值的一半,而奇数位则在前一个偶数位上$\operatorname{or} 2^{l-1}​$,其中$l​$是$2^l​$超过多项式最终长度的最小$l​$

例题:【模板】多项式乘法(FFT)

#include<cstdio> 
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<map>
#include<set>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db long double
#define fir first
#define sec second
#define mp make_pair
using namespace std;

namespace flyinthesky {

    const int MAXN = 4200000 + 5;
    const double PI = acos(-1);

    struct C { // 手写 complex 类
        double r, i;
        C() {r = i = 0;}
        C(double x, double y){r = x, i = y;}
        C operator + (C &x) {return C(r + x.r, i + x.i);}
        C operator - (C &x) {return C(r - x.r, i - x.i);}
        C operator * (C &x) {return C(r * x.r - i * x.i, r * x.i + i * x.r);}
        void operator += (C &x) {r += x.r, i += x.i;}
        void operator *= (C &x) {
            double t = r; 
            r = r * x.r - i * x.i, i = t * x.i + i * x.r;
        }
    }f[MAXN], g[MAXN];

    int n, m, r[MAXN];

    void FFT(C *a, int op) { // op = 1 DFT; op = -1 IDFT; 
        C W, w, t, *a0, *a1;
        for (int i = 0; i < n; ++i) // 根据映射关系交换元素
            if (i < r[i]) swap(a[i], a[r[i]]);
        for (int i = 1; i < n; i <<= 1) { // 第几层
            W = C(cos(PI / i), sin(PI / i) * op);
            for (int j = 0; j < n; j += (i << 1)) { // 一层中的子问题
                w = C(1, 0), a0 = a + j, a1 = a0 + i;
                for (int k = 0; k < i; ++k) { // 处理子问题
                    t = *a1 * w, *a1 = *a0 - t, *a0 += t; // 蝴蝶操作
                    w *= W, ++a0, ++a1;
                }
            }
        }
    }

    void clean() {
    }
    int solve() {

        clean();
        scanf("%d%d", &n, &m);
        for (int x, i = 0; i <= n; ++i) scanf("%d", &x), f[i].r = (db)x;
        for (int x, i = 0; i <= m; ++i) scanf("%d", &x), g[i].r = (db)x;

        int l = 0;
        for (m += n, n = 1; n <= m; n <<= 1, ++l);
        for (int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); //递推求r

        FFT(f, 1), FFT(g, 1);
        for (int i = 0; i < n; ++i) f[i] *= g[i]; // 对应 y 相乘
        FFT(f, -1);
        for (int i = 0; i <= m; ++i) printf("%.0f ",fabs(f[i].r) / n); // IDFT 要除以 n

        return 0;
    }
}
int main() {
    flyinthesky::solve();
    return 0;
}
/*
1 6
5 6
1 2 3 4 5 6 7
*/

参考文章

------ 本文结束 ------