简明的FFT分治法学习记录

目的

多项式$F(x)={\sum }_{i=0}^n{f}_i{x}^i$​ $,G(x)={\sum }_{i=0}^m{g}_i{x}^i$若要实现乘法，朴素情况下显然可以发现是$O(n^2)$级别的。

然而我们知道，一个次数为n的多项式，又可以被n个点值构成的集合表示。

通过巧妙选取点值（FFT），将多项式乘法转化为点值乘法，就可以将复杂度压缩至$O(nlogn)$级别。

简明理论背景

${\omega }_n^1$称为n次单位根，将圆n等分，第k个点就记为${\omega }_n^k$

$${\omega }_n^1=cos\frac{1}{n}2\pi +isin\frac{1}{n}2\pi$$

$${\omega }_n^k=cos\frac{k}{n}2\pi +isin\frac{k}{n}2\pi =({\omega }_n^1{)}^k$$

还要用到三个前置性质

$${\omega }_n^k={\omega }_{2n}^{2k}$$

$$\begin{gathered} {\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k \\ {\omega }_n^{k+n}={\omega }_n^k \end{gathered}$$

代入上式即可证，第二个式子可以理解为表示的向量等大反向（n/2即圆n等分的一半，即向量转了半圈），第三个就是转了一圈

FFT分治实现

$$\begin{gathered} F(x)=\sum _{i=0}^{n-1}{a}_i{x}^i=a_0+a_{1}x+a_2{x}^2+...+a_{n-1}{x}^{n-1} \\ =(a_0+a_2{x}^2+...+a_{n-2}{x}^{n-2})+x(a_1+a_3{x}^2+...+a_{n-1}{x}^{n-2}) \end{gathered}$$

$$\begin{gathered} 令F_L(x)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i}{x}^i,F_R(x)=\sum _{i=0}^{\frac{n}{2}-1}{a}_{2i+1}{x}^i \\ 则F(x)=F_L(x^2)+x{F}_R(x^2) \end{gathered}$$

利用单位根的性质进行代入

$$\begin{gathered} F({\omega }_n^i)=F_L({\omega }_n^{2i})+{\omega }_n^i{F}_R({\omega }_n^{2i})=F_L({\omega }_{\frac{n}{2}}^i)+{\omega }_n^i{F}_R({\omega }_{\frac{n}{2}}^i),0\le i\le \frac{n}{2}-1,i\in Z \\ F({\omega }_n^{i+\frac{n}{2}})=F_L({\omega }_{\frac{n}{2}}^i)-{\omega }_n^i{F}_R({\omega }_{\frac{n}{2}}^i),\frac{n}{2}\le i\le n-1,i\in Z \end{gathered}$$

发现$F_L({\omega }_{\frac{n}{2}}^i),F_R({\omega }_{\frac{n}{2}}^i)$只有前面的系数有所差异。即我们只要得到前n/2的点就可以推后n/2个点(知道前半圆就可以知道一整个圆)

这又显然符合分治结构。而求解$F_L({\omega }_{\frac{n}{2}}^i),F_R({\omega }_{\frac{n}{2}}^i)$后的合并是$O(n)$的。则有$T(n)=2T(\frac{n}{2})+\Theta (n)=\Theta (n\log n)$

递归实现

参考了http://www.longluo.me/blog/2022/04/02/FFT

struct Complex {
    double x, y;
    Complex(double _x = 0.0, double _y = 0.0) {
        x = _x, y = _y;
    }
    Complex operator -(const Complex &a) const {
        return Complex(x - a.x, y - a.y);
    }
    Complex operator +(const Complex &a) const {
        return Complex(x + a.x, y + a.y);
    }
    Complex operator *(const Complex &a) const {
        return Complex(x * a.x - y * a.y, x * a.y + y * a.x);
    }
}; // 封装了复数

const double PI = acos(-1.0);

vector<Complex> FFT(vector<Complex> &a) {
    // 实部为系数
    int n = a.size();
    if (n == 1) return a;
    vector<Complex> FL(n / 2), FR(n / 2);
    for (int i = 0; 2 * i < n; i++) { // 拆分奇偶
        FL[i] = a[2 * i];
        FR[i] = a[2 * i + 1];
    }
    vector<Complex> GL = FFT(FL), GR = FFT(FR); // 分治
    vector<Complex> G(n);
    double angle = 2 * PI / n; // omega^1_n
    Complex omega = Complex(cos(angle), sin(angle));
    Complex curRoot = Complex(1, 0);
    for (int i = 0; 2 * i < n; i++) { // 上述递推式
        G[i] = GL[i] + curRoot * GR[i];
        G[i + n / 2] = GL[i] - curRoot * GR[i];
        curRoot = curRoot * omega;
    }
    return G;
}   

实例应用

大数乘法https://leetcode.cn/problems/multiply-strings/

#include <iostream>
#include <algorithm>
#include <cmath>
#include <vector>
using namespace std;

struct Complex {
    double x, y;
    Complex(double _x = 0.0, double _y = 0.0) {
        x = _x, y = _y;
    }
    Complex operator -(const Complex &a) const {
        return Complex(x - a.x, y - a.y);
    }
    Complex operator +(const Complex &a) const {
        return Complex(x + a.x, y + a.y);
    }
    Complex operator *(const Complex &a) const {
        return Complex(x * a.x - y * a.y, x * a.y + y * a.x);
    }
};

const double PI = acos(-1.0);

vector<Complex> FFT(vector<Complex> &a, bool flag) {
    // 实部为系数
    int n = a.size();
    if (n == 1) return a;
    vector<Complex> FL(n / 2), FR(n / 2);
    for (int i = 0; 2 * i < n; i++) {
        FL[i] = a[2 * i];
        FR[i] = a[2 * i + 1];
    }
    vector<Complex> GL = FFT(FL, flag), GR = FFT(FR, flag);
    vector<Complex> G(n);
    double angle = 2 * PI / n * (flag ? -1 : 1);
    Complex omega(cos(angle), sin(angle));
    Complex curRoot(1, 0);
    for (int i = 0; 2 * i < n; i++) {
        G[i] = GL[i] + curRoot * GR[i];
        G[i + n / 2] = GL[i] - curRoot * GR[i];
        curRoot = curRoot * omega;
    }
    return G;
}   

string multiply(string num1, string num2) {
    if (num1 == "0" || num2 == "0") return "0";
    int len1 = num1.size(), len2 = num2.size();
    int n = 1;
    while (n < len1 + len2) n <<= 1;
    vector<Complex> a(n), b(n);
    for (int i = len1 - 1; i >= 0; i--) {
        a[i] = Complex(num1[len1 - 1 - i] - '0', 0);
    }
    for (int i = len2 - 1; i >= 0; i--) {
        b[i] = Complex(num2[len2 - 1 - i] - '0', 0);
    }
    a = FFT(a, 0), b = FFT(b, 0);
    for (int i = 0; i < n; i++)
        a[i] = a[i] * b[i];
    a = FFT(a, 1);
    string ans;
    int carry = 0;
    for (int i = 0; i < n; i++) {
        int sum = round(round(a[i].x) / n) + carry;
        carry = sum / 10;
        ans += sum % 10 + '0';
    }
    if (carry > 0) {
        ans += carry % 10 + '0';
    }
    int ind = ans.size() - 1;
    while (ans[ind] == '0' && ind > 0) {
        ind --;
    }
    ans = ans.substr(0, ind + 1);
    reverse(ans.begin(), ans.end());
    return ans;
}

int main() {
    string a, b;
    cin >> a >> b;
    cout << multiply(a, b) << endl;
    return 0;
}

参考

1. http://www.longluo.me/blog/2022/04/02/FFT
2. https://www.cnblogs.com/JustinRochester/p/17036805.html
3. https://blog.csdn.net/enjoy_pascal/article/details/81478582/