FFT学习笔记

恶心人的多项式

快速傅里叶变换

多项式

定义

不严谨定义：形如 $f(x) = \sum \limits _{i=0}^{n} a_ix^i$ 的式子为多项式。

定义（from OIWiki）：对于求和式 $\sum a_nx^n$，如果是有限项相加，称为多项式，记作 $f(x)=\sum \limits_{n=0}^m a_nx^n$。

次数：对于多项式 $F(x) = \sum \limits _{i=0}^{n} a_ix^i$，该多项式的次数为 $n$。表示为 $\text{degree} (F) = n$

次数界：任何一个大于多项式次数的整数。

基本运算

系数表示法：

$$\mathbf{a} = [a_1, a_2,\cdots, a_n] ^\top$$

加法

$$\mathbf{c} = \mathbf{a} + \mathbf{b}$$

$$c_i = \sum \limits_{i = 0}^{n} (a_i + b_i)$$

时间复杂度 $O(n)$。

乘法

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

$$\text{degree} (C) = \text{degree} (A) + \text{degree} (B)$$

$$c_i = \sum \limits_{j = 0}^i a_j b_{i - j}$$

记为 $\mathbf{c} = \mathbf{a} \otimes \mathbf{b}$

时间复杂度 $O(n^2)$。

求值

秦九韶算法：

$$\begin{aligned} f(x) ={}& a_nx^n+a_{n-1}x^{n-1}+\cdots a_1x+a_0 \\ ={}& (a_nx^{n-1}+a_{n-1}x^{n-2}+\cdots a_1)x+a_0 \\ ={}& ((\cdots(a_nx+a_{n-1})x+a_{n-2})x+\cdots))x+a_0 \end{aligned}$$

时间复杂度 $O(n)$。

点值表示法：

$$\{(x_0, A(x_0)), (x_1, A(x_1)), \cdots(x_n, A(x_n))\}$$

加法

$$C(x_i) = A(x_i) + B(x_i)$$

时间复杂度 $O(n)$。

乘法

$$C(x_i) = A(x_i)B(x_i)$$

时间复杂度 $O(n)$。

求值

拉格朗日插值法：

$$A(x)=\sum \limits_{i=1}^n A(x_i) \frac{\prod \limits _{j\ne i}^{n}(x- x_j)}{\prod \limits_{j \ne i}^{n}(x_i-x_j) }$$

时间复杂度 $O(n^2)$

代码实现：

double res = 0;
for(int i = 1; i <= n; i ++ )
{
    double s1 = 1, s2 = 1;
    for(int j = 1; j <= n; j ++ )
        if(j != i)
        {
            s1 = s1 * (k - x[j]);
            s2 = s2 * (x[i] - x[j]);
        }

    res = res + y[i] * s1 / s2;
}

复数

基本运算

高中数学基础，不多介绍

$$(a+bi) + (c+di) = (a+c) + (b + d)i$$

$$(a + bi) - (c + di) = (a - c) + (b - d)i$$

$$(a + bi)(c + di) = (ac - bd) + (ad + bc)i$$

$$\frac{a + bi}{c + di} = \frac{ac+bd}{c^2 + d^2} + \frac{bc-ad}{c^2+d^2}i$$

欧拉定理

欧拉定理：$e^{i\theta} = \cos \theta + i\sin \theta$

根据这个定理，我们可以把复数转化成在复平面上的一堆点来进行一些公式的推导。

算术基本定理

$n$ 次代数方程 $f(x)=a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 = 0(a_n \ne 0 且 n > 0)$ 在复数域上恰好有 $n$ 个根。

根据这个定理，我们可以找出一个特殊的式子：$x^n=1$ 的根，在复数域中恰好有 $n$ 个，设 $\omega _n = e^{\frac{2\pi}{n}i}$，则 $x^n = 1$ 的解集可以记作 $\{ \omega _n^k|k=0,1,2,\cdots,n-1 \}$。

单位根

单位根具有一些性质，在 FFT 的时候会有很大用处。

折半引理

$$\omega _{2n}^{k+n} = -\omega _{2n}^{k}$$

这个结论在复平面中很好证明，当然也可以用公式来证明：

$$\begin{aligned} \omega _{2n}^{k+n} &= e^{i\frac{2\pi}{2n}(k+n)} \\ &=e^{i\frac{2\pi}{2n}k+i\pi} \\ &= e^{i\frac{2\pi}{2n}k}\times(-1) \\ &= -\omega _{2n}^k \end{aligned}$$

消去引理

$$\omega _{dn}^{dk} = \omega _{n}^{k}$$

证明：

$$\begin{aligned} \omega _{dn}^{dk} &= e^{i\frac{2\pi}{dn}dk } \\ &=e^{i\frac{2\pi}{n} k} \\ &= \omega ^{k}_{n} \end{aligned}$$

求和引理

$$\sum \limits_{i=0}^{n-1} \omega_n^i= 0$$

证明：

$$\begin{aligned} \sum \limits_{i=0}^{n-1} \omega_n^i &= \frac{(\omega_n^k)^n -1}{\omega _n^k-1} \\ &=\frac{(\omega _n^n)^k-1}{\omega _n^k-1} \\ &= 0 \end{aligned}$$

FFT

在计算多项式乘法时，系数表示法的复杂度为 $O(n^2)$，而点值表示法的复杂度仅为 $O(n)$，我们可以尝试将系数表示法转化为点值表示法来进行运算。

我们可以对多项式进行如下变换，将奇数项和偶数项分开：

（前提：$n$ 为奇数且 $n+1$ 为 $2$ 的幂）

$$\begin{aligned} A(x) &= a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0\\ &= a_nx^n+a_{n-2}x^{n-2}+\cdots+a_1x+a_{n-1}x^{n-1}+a_{n-3}x^{n-3}+\cdots+a_0\\ \end{aligned}$$

将奇数项的系数构成的多项式记为 $A^{[1]}$，偶数项记为 $A^{[0]}$，就可以得到下面的式子。

$$\begin{aligned} A(x)&= a_nx^n+a_{n-2}x^{n-2}+\cdots+a_1x+a_{n-1}x^{n-1}+a_{n-3}x^{n-3}+\cdots+a_0\\ &= A^{[0]}(x^2)+xA^{[1]}(x^2) \end{aligned}$$

这时候我们想到，可以对两边的多项式分治处理，但是实数中没有可以把这个平方消掉的数，但是我们可以从复数中找到。

将 $\omega _{2n}^k$ 带入 $A$ 中，就可以得到：

$$\begin{aligned} A(\omega _{2n}^k)&= A^{[0]}(\omega _{2n}^{2k})+\omega _{2n}^kA^{[1]}(\omega _{2n}^{2k}) \\ &= A^{[0]}(\omega _n^k)+\omega _{2n}^kA^{[1]}(\omega _{n}^{k}) \end{aligned}$$

我们就可以对 $A^{[0]}$ 和 $A^{[1]}$ 进行递归处理。

同时，在复平面中，$\omega _{2n}^{n+k}$ 与 $\omega _{2n}^{k}$ 在多项式 $A^{[0]}(x^2)$ 和 $A^{[1]}(x^2)$ 的值相同，因此我们也可以算出该点的值。

$$\begin{aligned} A(\omega _{2n}^{n+k})&= A^{[0]}(\omega _n^k)+\omega _{2n}^{n+k}A^{[1]}(\omega _{n}^{k})\\ &=A^{[0]}(\omega _n^k)-\omega _{2n}^{k}A^{[1]}(\omega ^k_n) \end{aligned}$$

这就是 FFT 的实现原理。做完 FFT 之后，我们就得到了一个多项式在复数域中 $2n$ 个点的点值，这时再将两个多项式乘起来，就可以得到乘积的点值表示。

但是如果我们需要的是系数表示，我们则需要做一遍快速傅里叶逆变换，将点值转化为系数。

我们可以列出刚才求点值时的矩阵表示：

$$\begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}=\begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix}\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}$$

我们现在只需要求出逆矩阵即可。由于这个矩阵元素非常特殊（没记错应该叫范德蒙德矩阵），它的逆矩阵即为每一项的倒数再除以变换的长度。

为了取倒数，我们可以做变换：

$$\begin{aligned} \omega _k^{-1}&= e^{-\frac{2\pi}{k}i } \\ &= \cos (-\frac{2\pi}{k})+\sin(-\frac{2\pi}{k})\\ &=\cos(\frac{2\pi}{k})-\sin(\frac{2\pi}{k}) \end{aligned}$$

我们可以看到与 FFT 时仅差了一个负号，因此我们可以将 FFT 和 IDFT 结合为一题，传参时传入一个 $opt$，为 $1$ 时执行 FFT，为 $-1$ 时执行 IDFT。

代码实现：

#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
const int N = 3e6 + 10;
struct Complex
{
    double a, b;
    Complex() {a = 0, b = 0;}

    Complex(double real, double imag): a(real), b(imag) {}

    Complex operator + (const Complex& x) const
    {
        return Complex(a + x.a, b + x.b);
    }

    Complex operator - (const Complex& x) const
    {
        return Complex(a - x.a, b - x.b);
    }

    Complex operator * (const Complex& x) const
    {
        return Complex(a * x.a - b * x.b, a * x.b + b * x.a);
    }
}F[N], G[N];

void FFT(Complex a[], int lim, int opt)
{
    if(lim == 1) return;
    Complex a1[lim >> 1], a2[lim >> 1];
    for(int i = 0; i < lim; i += 2)
        a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    FFT(a1, lim >> 1, opt);
    FFT(a2, lim >> 1, opt);
    Complex wn = Complex(cos(2.0 * PI / lim), opt * sin(2.0 * PI / lim));
    Complex w = Complex(1, 0);
    for(int k = 0; k < (lim >> 1); k ++ )
    {
        a[k] = a1[k] + w * a2[k];
        a[k + lim / 2] = a1[k] - w * a2[k];
        w = w * wn;
    }
}

int n, m;

int main()
{
    n = read(), m = read();
    for(int i = 0; i <= n; i ++ ) scanf("%lf", &F[i].a);
    for(int i = 0; i <= m; i ++ ) scanf("%lf", &G[i].a);

    int lim = 1;
    while(lim <= n + m) lim <<= 1;
    FFT(F, lim, 1), FFT(G, lim, 1);

    for(int i = 0; i <= lim; i ++ ) F[i] = F[i] * G[i];
    FFT(F, lim, -1);

    for(int i = 0; i <= n + m; i ++ ) cout << (int)(F[i].a / lim + 0.5) << ' ';

    return 0;
}

优化

BitReverse

上面我们使用了递归实现 FFT，我们可以思考一下能否用迭代来实现 FFT。

我们可以观察一下每次递归，都是哪些元素被分到了一组。接下来拿 $0\sim 7$ 来举例。

$$\begin{aligned}(0,1,2,3,4,5,6,7)\\(0,2,4,6)(1,3,5,7)\\(0,4)(2,6)(1,5)(3,7)\\0,4,2,6,1,5,3,7\end{aligned}$$

乍一看好像毫无规律，但是我们把最后分的一组中二进制表示出来：

$$000,100,010,110,001,101,011,111$$

（以下所说的数字全为二进制）

我们可以发现，$000$ 处在了第 $000$ 号位，$100$ 处在了第 $001$ 号位，$010$ 处在了第 $010$ 号位，$110$ 处在了第 $011$ 号位…

结论已经显然易见了，把二进制位翻转一下即为 FFT 时的排列顺序，因此我们可以通过这个操作来增加效率，我们把这个操作称为位逆序置换。

不难发现其中有 $4$ 条规律：

1. 奇数位二进制首位为 $1$，偶数位二进制首位为 $0$。
2. 每两个二进制数除了首位都相同。
3. 每两个数除首位外是上一级子问题的解
4. 前一半数除末位外是上一级子问题的解

因此我们可以写出下面的递推代码：rev[i] = (rev[i >> 1] >> 1) | (i & 1) << (len - 1)

蝴蝶操作

在代码中，有一段操作：

a[k] = a1[k] + w * a2[k];
a[k + lim / 2] = a1[k] - w * a2[k];

其中都包含了共同的一项 w * a2[k]，我们可以将这一项单独存储，减小常数。

同时我们也可以将 a1，a2 数组省去，直接在原数组上进行操作。

代码实现：

#define LOCAL
#include <bits/stdc++.h>
using namespace std;

const double PI = acos(-1);
int lim, len;
int rev[N];

void FFT(Complex a[], int opt)
{
    for(int i = 0; i < lim; i ++ ) 
        if(i < rev[i])
            swap(a[i], a[rev[i]]);
    int up = log2(lim);
    for(int dep = 1; dep <= up; dep ++ )
    {
        int m = 1 << dep;
        Complex wn = Complex(cos(2.0 * PI / m), opt * sin(2.0 * PI / m));
        for(int k = 0; k < lim; k += m)
        {
            Complex w = Complex(1, 0);
            for(int j = 0; j < m / 2; j ++ )
            {
                Complex t = w * a[k + j + m / 2];
                Complex u = a[k + j];
                a[k + j] = u + t;
                a[k + j + m / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if(opt == -1)
    {
        for(int i = 0; i < lim; i ++ ) a[i].a /= lim;
    }
}

int n, m;

int main()
{
    n = read(), m = read();
    for(int i = 0; i <= n; i ++ ) scanf("%lf", &F[i].a);
    for(int i = 0; i <= m; i ++ ) scanf("%lf", &G[i].a);

    lim = 1;
    while(lim <= n + m) lim <<= 1, len ++;
    for(int i = 0; i < lim; i ++ )  
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len - 1));

    FFT(F, 1), FFT(G, 1);

    for(int i = 0; i <= lim; i ++ ) F[i] = F[i] * G[i];
    FFT(F, -1);

    for(int i = 0; i <= n + m; i ++ ) cout << (int)(F[i].a + 0.5) << ' ';

    return 0;
}