FFT快速傅里叶变换（多项式乘法）

本文仅供作者本人练习推导FFT，如有侵权请私信我

FFT算法过程

FFT的精髓在于将 $n-1({\log }_{2}n\in \mathbb{N})$ 次的 $n$ 项多项式 $A(x)$ 在 $O(n\log n)$ 的时间（而不是硬算的 $O(n^2)$ ）内求出其在 $n$ 个 $n-1$ 次单位根 ${\omega }_{n-1}^i(i=0,1,\dots ,n-1)$ 上的取值，也就是将任意多项式O(nlogn)从系数表示法转化为点值表示法

推式子

$$\begin{gathered} 设A(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_{n-1}{x}^{n-1} \\ 令A_0(x)=a_0+a_{2}x+a_4{x}^2+\cdots +a^{\frac{n}{2}-1}, \\ {A}_1(x)=a_1+a_{3}x+a_5{x}^2+\cdots +a^{\frac{n}{2}-1}. \\ 则显然\begin{array}{c}A(x)=A_0(x^2)+x{A}_1(x^2)\end{array} \\ 将x={\omega }_{n-1}^i代入(1)有：\begin{array}{c}A({\omega }_{n-1}^i)=A_0({\omega }_{n-1}^{2i})+{\omega }_{n-1}^i{A}_1({\omega }_{n-1}^{2i})(\end{array} \\ \Rightarrow \begin{array}{c}A({\omega }_{n-1}^i)=A_0({\omega }_{(n-1)/2}^i)+{\omega }_{n-1}^i{A}_1({\omega }_{(n-1)/2}^i)\end{array} \\ 将x={\omega }_{n-1}^{i+(n-1)/2}代入(1)有：\begin{array}{c}A({\omega }_{n-1}^i)=A_0({\omega }_{(n-1)/2}^{2i+(n-1)})+{\omega }_{n-1}^{i+(n-1)/2}{A}_1({\omega }_{n-1}^{2i+(n-1)})(\end{array} \\ \Rightarrow \begin{array}{c}A({\omega }_{n-1}^i)=A_0({\omega }_{(n-1)/2}^i)-{\omega }_{n-1}^i{A}_1({\omega }_{(n-1)/2}^i)\end{array} \end{gathered}$$
注意到
$$\{\begin{array}{l}A({\omega }_{n-1}^i)=A_0({\omega }_{(n-1)/2}^i)+{\omega }_{n-1}^i{A}_1({\omega }_{(n-1)/2}^i)\\ A({\omega }_{n-1}^i)=A_0({\omega }_{(n-1)/2}^i)-{\omega }_{n-1}^i{A}_1({\omega }_{(n-1)/2}^i)\end{array}$$
所以只需要求出 $A_0(x)$ 和 $A_1(x)$ 的各 $n/2$ 个点值就可以得到 $A(x)$ 的 $n$ 个点值。看似求的点值没变，但是只需要进行 $\log n$ 轮操作，所以时间复杂度是低的
附记：此处推式子应该用 $n$ 而不是 $n-1$ 的，会更好看……

代码实现

实现代码时，建议使用非递归写法来节省空间和行数
【例题：】O(nlogn)高精度乘法

#include <bits/stdc++.h>
using namespace std;
const double PI = acos(-1);
struct Co {
    double Re, Im;
};
Co operator+(Co A, Co B) {return Co{A.Re + B.Re, A.Im + B.Im};}
Co operator-(Co A, Co B) {return Co{A.Re - B.Re, A.Im - B.Im};}
Co operator*(Co A, Co B) {return Co{A.Re*B.Re-A.Im*B.Im, A.Im*B.Re+A.Re*B.Im};}
void FFT(int n, vector<Co> &a, int type);
void init();
string s1, s2, o;
int n, m, lim = 1, lg;
vector<Co> a, b, ans;
vector<int> rev;
int main()
{
    // freopen("D.in", "r", stdin);
    // freopen("D.out", "w", stdout);
    // ios::sync_with_stdio(0);
    // cin.tie(0);
    cin >> s1 >> s2;
    n = s1.size();
    m = s2.size();
    for (int i = n - 1; i >= 0; --i) a.push_back(Co{double(s1[i] - '0'), 0});
    for (int i = m - 1; i >= 0; --i) b.push_back(Co{double(s2[i] - '0'), 0});

    init();
    FFT(lim, a, 1);
    FFT(lim, b, 1);

    for (int i = 0; i < lim; ++i) {
        ans[i] = a[i] * b[i];
    }
    FFT(lim, ans, -1);
    int carry = 0;
    for (int i = 0; i < lim; ++i) {
        int t = (int)(ans[i].Re / lim + 0.5);
        t += carry;
        carry = t / 10;
        t %= 10;
        o += char(t + '0');
    }
    int end = o.size() - 1;
    while (o[end] == '0' && end > 0) --end;
    while (~end) {
        putchar(o[end]);
        --end;
    }
    return 0;
}
void init() {
    while (lim <= n + m) {
        lim <<= 1;
        ++lg;
    }
    rev.resize(lim, 0);
    for (int i = 1; i < lim; ++i) rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (lg - 1));
    a.resize(lim, Co{0, 0});
    b.resize(lim, Co{0, 0});
    ans.resize(lim, Co{0, 0});
}
void FFT(int n, vector<Co> &a, int type) {
    // 二进制重排
    for (int i = 1; i < n; ++i) {
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    }
    
    for (int mid = 1; mid < n; mid <<= 1) {
        Co Wn{cos(PI / mid), type * sin(PI / mid)};
        for (int l = 0; l < n; l += (mid << 1)) {
            Co w{1, 0};
            for (int i = 0; i < mid; ++i) {
                Co tmp = w * a[l + mid + i];
                a[l + mid + i] = a[l + i] - tmp; // 新右=左-w*右
                a[l + i] = a[l + i] + tmp; // 新左=左+w*右
                w = w * Wn;
            }
        }
    }
}