【Coel.学习笔记】【调整一下状态】快速傅里叶变换（Fast Fourier Transform）

学了好几天的网络流，现在满脑子都是源点、汇点、残留网络、最小割……
今天剩下了挺多时间，调整一下学学 FFT。想当年越级学这玩意还学得半懂不懂……
警告：以下内容用到了大量数学公式，请做好心理准备。此外，在阅读本篇博客时，最好确保您已经学习过北师大版高中必修二数学课本的一至五章，包括选学内容，否则阅读本博客会存在困难。

引入

快速傅里叶变换（$\mathcal{Fast\, Fourier\, Transform}$）是用来快速求两个多项式的乘积，即卷积的算法。直接展开求解的复杂度是 $O(n^2)$，而快速傅里叶变换可以把结果优化到 $O(n\log n)$，得到了极大的优化。

前置知识

多项式的两种表示方法

对于一个 $n-1$ 次多项式 $A(x)=a_0+a_1x+a_2x^2+...+a_nx^n$，可以用 $n$ 个不同点唯一确定，被称为多项式的点值表示法。对应地，上述表示法叫做系数表示法。
可以看作是一个 $n-1$ 次函数，这就和初中时学到的函数一样，用待定系数法求出多项式了。

点值表示法有什么用呢？假设还是求 $A(x)\times B(x)$，其中 $A(x)$ 为 $n$ 次,$B(x)$ 为 $m$ 次，则取 $n+m+1$ 个不同点，那么乘积得到的新多项式可以一步到位求出，时间复杂度是 $O(n+m)$。

那么我们的目标就是尽可能快地将一个多项式在系数表示法和点值表示法间相互转换。

单位根

在高中数学学过，一个复数可以表示成 $a+bi$，$a,b$ 分别是复数的实部和虚部；复数和复平面上的点一一对应，复数的模长等于 $\sqrt{a^2+b^2}$。复数对应的向量与实轴正半轴的夹角叫做辐角。

在复数域上取一个单位圆并且等分为 $n$ 份，每一份对应的复数叫做 $n$ 次单位根，第 $k$ 个记作 $\omega^k_n$。另一种单位根的定义为：方程 $x^n=1$ 在复数范围内的解，第 $k$ 个解记作 $\omega^k_n$。根据欧拉公式 $e^{ix}=\cos x + i\sin x$，可以得到单位根的通项表示 $\omega_n=e^{\frac{2\pi i}{n}}=\cos\dfrac{2\pi}{n}+i\sin\dfrac{2\pi}{n}$。

单位根的性质：

1. $\omega^k_n=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}$。本质上是复数的三角表示法的转换。
2. $\omega^n_n=\omega^0_n=1$。这和单位圆的周期性有关，也可以从上面式子推得。
3. $\omega^{2k}_{2n}=\omega^k_n$。这很显然，因为分成 $n$ 份取 $k$ 份和分成 $2n$ 份取 $2k$ 份本质是一样的。
4. $\omega^{k+\frac{n}{2}}_n=-\omega^k_n$。加上 $\dfrac{n}{2}$ 其实就是旋转了 $\dfrac{2n\pi}{2n}=\pi$。还记得诱导公式吗？

记忆性质时请务必理解记忆。

算法推导

利用单位根的性质，我们对 $n-1$ 次多项式取一组 $n-1$ 次单位根，并将多项式按奇偶排列：

$$A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$$

拆解一下，令 $x$ 表示原来的 $x^2$：

$$A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1}$$

提取 $x$ 再用 $x$ 表示原来的 $x^2$：

$$A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1}$$

这样就会有 $A(x)=A_1(x^2)+xA_2(x^2)$。
那么，当 $k\in [0,\dfrac{n}{2}-1)$ 时，有 $A(\omega^k_n)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega^{2k}_n)$。这里利用到了复数三角表示法的运算法则。
]这里打个右括号是因为博客园把左括号看成链接了= =

再利用性质 3 化简，原式等于

$$A_1(\omega^k_\frac{n}{2})+\omega^k_nA_2(\omega^k_\frac{n}{2})$$

看看 $\dfrac{n}{2}$ 到 $n-1$ 的部分。
$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega^{2k}_{n})+\omega_n^{k+\frac{n}{2}}A_2(\omega^{2k}_{n})$
$=A_1(\omega^{k}_{\frac{n}{2}})-\omega_n^{k}A_2(\omega^{k}_{\frac{n}{2}})$（运用性质 3,4）

发现了什么？这两个式子正好满足前一项相等，后一项相反！

也就是说，我们可以通过分治的方式推出值。由于这个分治算法能处理的长度只能为 $2$ 的正整数次幂，所以要在变换之前补零占位。

现在看看逆变换。不加证明地给出结论：记 $(\omega^k_n,A(\omega^k_n))$ 为 $(x_k, y_k)$，则多项式相乘得到的第 $k$ 项系数

$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$

接下来继续推导：

$$c_k=\sum_{i=0}^{n-1}（\sum^{n-1}_{j=0}a_j\omega_n^i)(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}\sum^{n-1}_{j=0}a_j(\omega_n^{j-k})^i$$

$$=\sum_{i=0}^{n-1}a_i\sum^{n-1}_{j=0}(\omega_n^{j-k})^i$$

$$=na_k$$

最后一步有点复杂，细说一下：
构造一个多项式 $S(x)= 1+x^1+x^2+...+x^{n-1}$。当 $k\ne0$ 时，$S(\omega^k_n)=\omega^0_n+\omega^k_n+\omega^{2k}_n+...+\omega^{(n-1)k}_n$。提取一个 $\omega^k_n$，有

$$\omega^k_nS(\omega^k_n)=\omega^k_n+\omega^{2k}_n+...+\omega^0_n$$

这两个式子式是相等的，所以可以相减得 $S=0$。

当 $k=0$ 时，$S(\omega^k_n)=S(1)=n$。代入原式情况，就有结果为 $na_k$了。

实现时如果递归实现会很慢，所以继续优化成迭代。

对于一组系数：

$$a_0,a_1,a_2,a_3,a_4,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_6,a_1,a_5,a_3,a_7$$

这样分组本质上就是递归分治的过程。

现在观察一下它们下标的二进制变化：

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

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

可以发现二进制下标实际上翻转了一遍。

而在实际求解中，我们可以以 $O(n)$ 的复杂度递推出翻转后的值，这也就保证了快速傅里叶变换的效率。
不加证明地给出结果：
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));

顺带一提，这个二进制反转的过程叫做位逆序置换，或者蝴蝶变换。

代码实现

【模板】多项式乘法（FFT）

// Problem: P3803 【模板】多项式乘法（FFT）
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3803
// Memory Limit: 500 MB
// Time Limit: 2000 ms
// Author: Coel
// 
// Powered by CP Editor (https://cpeditor.org)

#include <cmath>
#include <cstring>
#include <iostream>

using namespace std;


const int maxn = 3e6 + 10;
const double pi = acos(-1.0);

struct complex {
    double x, y;
    complex operator+(const complex &t) const { return {x + t.x, y + t.y}; }
    complex operator-(const complex &t) const { return {x - t.x, y - t.y}; }
    complex operator*(const complex &t) const { return {x * t.x - y * t.y, x * t.y + y * t.x}; }
} a[maxn], b[maxn];

int rev[maxn], bit, tot;

void init_bit() { //预处理翻转
    while ((1 << bit) < n + m + 1) bit++;
    tot = 1 << bit; //二进制移位
    for (int i = 0; i < tot; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
}

void fft(complex a[], int inv) { //inv 表示变换/逆变换，由于两种变换只有辐角的区别所以可以这么写
    for (int i = 0; i < tot; i++)
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int mid = 1; mid < tot; mid <<= 1) { //按照推导过程迭代
        auto w1 = complex({cos(pi / mid), inv * sin(pi / mid)});
        for (int i = 0; i < tot; i += mid * 2) {
            auto wk = complex({1, 0});
            for (int j = 0; j < mid; j++, wk = wk * w1) {
                auto x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

int main(void) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n, m;
    cin >> n >> m;
    for (int i = 0; i <= n; i++)
        cin >> a[i].x;
    for (int i = 0; i <= m; i++)
        cin >> b[i].x;
    init_bit();
    fft(a, 1), fft(b, 1);
    for (int i = 0; i < tot; i++)
        a[i] = a[i] * b[i];
    fft(a, -1);
    for (int i = 0; i <= n + m; i++)
        cout << (int)(a[i].x / tot + 0.5) << ' '; //四舍五入
    return 0;
}

【模板】A*B Problem 升级版（FFT 快速傅里叶变换）

高精度乘法的 fft 需要把高精度乘法转化为卷积形式。
对于一个大整数 $A$，设第 $i$ 位为 $a_i$，显然

$$A=a_n\times 10^n+a_{n-1}\times 10^{n-2}+...+a_0\times 10^0$$

把 $x$ 用 $10$ 代换，就可以轻松得到一个卷积了。

// Problem: P1919 【模板】A*B Problem 升级版（FFT 快速傅里叶变换）
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P1919
// Memory Limit: 256 MB
// Time Limit: 1500 ms
// Author: Coel
// 
// Powered by CP Editor (https://cpeditor.org)

#include <cmath>
#include <cstring>
#include <iostream>

using namespace std;

const int maxn = 5e6 + 10;
const double pi = acos(-1.0);

struct complex {
    double x, y;
    complex operator+(const complex &t) const { return {x + t.x, y + t.y}; }
    complex operator-(const complex &t) const { return {x - t.x, y - t.y}; }
    complex operator*(const complex &t) const { return {x * t.x - y * t.y, x * t.y + y * t.x}; }
} a[maxn], b[maxn];

int n, m, k;
int ans[maxn];
int rev[maxn], bit, tot;
char s1[maxn], s2[maxn];

void init_bit() {
    while ((1 << bit) < n + m + 1) bit++;
    tot = 1 << bit;
    for (int i = 0; i < tot; i++)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1));
}

void fft(complex a[], int inv) {
    for (int i = 0; i < tot; i++)
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int mid = 1; mid < tot; mid <<= 1) {
        auto w1 = complex({cos(pi / mid), inv * sin(pi / mid)});
        for (int i = 0; i < tot; i += mid * 2) {
            auto wk = complex({1, 0});
            for (int j = 0; j < mid; j++, wk = wk * w1) {
                auto x = a[i + j], y = wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

int main(void) {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin >> s1 >> s2;
    n = strlen(s1) - 1, m = strlen(s2) - 1;
    for (int i = 0; i <= n; i++)//值转化为多项式
        a[i].x = s1[n - i] - '0';
    for (int i = 0; i <= m; i++)
        b[i].x = s2[m - i] - '0';
    init_bit();
    fft(a, 1), fft(b, 1);
    for (int i = 0; i < tot; i++)
        a[i] = a[i] * b[i];
    fft(a, -1);
    for (int i = 0, t = 0; i < tot || t; i++) { // 多项式转化为值
        t += a[i].x / tot + 0.5;
        ans[k++] = t % 10;
        t /= 10;
    }
    while (k > 1 && !ans[k - 1]) k--; //删除前导零
    for (int i = k - 1; i >= 0; i--)
        cout << ans[i];
    return 0;
}

三转二优化

上面的卷积需要做三次变换（$A$ 转点值，$B$ 转点值， $A\times B$ 转系数），可不可以做个优化呢？
答案是肯定的。我们把 $B$ 的值放在 $A$ 的虚部上，让 $A$ 自我相乘，再转换回来，就只需要两次变换。这时答案在虚部上，并且要除上 $2$。
另外我们还可以把三角函数的值预处理出来，进一步减小常数。代码略。