「快速傅里叶变换」总结

前言

作为 NOI 大纲里的十级算法，FFT拥有着逼格十分高的名称，也十分令人头疼？（雾
通俗地讲，FFT旨在以优异的时间复杂度，解决两个多项式相乘的问题，即：
$c_k=\sum_{i=0}^{k}a_ib_{k-i}$
FFT的核心思想在于：将一个多项式转化成点值表示法，再由点值表示法推出原多项式。

~~如果觉得抽象可以看图~~ 在这里插入图片描述

前置知识

多项式的点值表示法和系数表示法
向量
复数

复数可以自行查阅 Oi-wiki。
这里简略介绍一下多项式表示法（不妨令 $A (x)$ 表示一个 $n - 1$ 次的多项式）。
系数表示法：若 $A(x)=\sum_{i=0}^{n}a_i\times x^{i}$ 。
点值表示法：不妨将 $n$ 个互不相同的 $x$ 代入多项式得到 $n$ 个 $y$ ，则 $A$ 可被这 $n$ 个值 唯一确定，其中 $y_i=\sum_{j=0}^{n-1}a_j\times x_{i}^{j}$ 。

单位根

请注意单位根的概念是建立在 复平面 上的，因此需要熟知复数的运算法则。

定义：以原点为圆心作单位圆，并钦定圆心为起点，圆点的 $n$ 等分点为终点，作 $n$ 条向量，这 $n$ 条向量即为 $n$ 次单位根，通常将幅角为正且最小的向量对应的复数称为 $\omega_n$ 。

一般地， $n$ 次单位根可由 欧拉公式 直接得出： $\omega_{n}^{k}=cosk\frac{2\pi}{n}+isink\frac{2\pi}{n}$

$\omega_n$ 显然满足以下性质：

每乘一次 $\omega_n$ 便逆时针转动 $\frac{2\pi}{n}$ 角度
$\omega_{n}^{n} = \omega_{n}$
$\omega_{2n}^{2k}=\omega_{n}^{k}$

快速傅里叶变换（DFT）

下皆令 $n=2^k(k\in Z)$
不妨令一个多项式 $A(x)=a_0+a_1x+a_2x^{2}+a_3x^{3}+······+a_{n-1}x^{n-1}$ 。
将其按照下标奇偶分类，可得：
$A_1(x)=a_0+a_2x+a_4x^{2}+······+a_{n-2}x^{\frac{n}{2}}$ $A_2(x)=a_1+a_3x+a_5x^{2}+······+a_{n-1}x^{\frac{n}{2}}$
显然有， $A(x)=A_1(x^2)+xA_{2}(x^2)$ 。
注意到，多项式 $A (x)$ 可被 $n$ 个值唯一确定，所以，如果我们可以利用上述式子，代入特殊的 $x\prime$ 值，以高效的时间复杂度得到 $A(x\prime)$ 对应的值，那么便可以用 点值表示法 表示 $A (x)$ 了。

代入 $\omega_{n}^{k}$ ，则有：
$A(\omega_{n}^{k})=A_1(\omega_{n}^{2k})+\omega_{n}^{k}A_2(\omega_{n}^{2k})$ $A(\omega_{n}^{k})=A_1(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k})$
代入 $\omega_{n}^{k+\frac{n}{2}}$ ，同理可得：
$A(\omega_{n}^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}A_2(\omega_{\frac{n}{2}}^{k})$

通过大眼观察法可以得出，当我们计算出 $A(\omega_{n}^{k})$ 时，我们同时也得到了 $A(\omega_{n}^{k+\frac{n}{2}})$ 的值，这便是又名的 蝴蝶操作。
于是，我们把原问题分割成了两个等价的子问题。这也就意味着，通过不断递归，我们将 $n$ 个单位根 $\omega_{n}^{k}$ 依次代入，便可以实现以 $O (n l o g n)$ 的时间复杂度用 点值表示法 表示多项式 $A (x)$ ！

快速傅里叶逆变换（IDFT）

然而，题目通常要求我们用 系数表示法 表示一个多项式，而非 点值表示法。因此，我们需要将点值转化为系数。

结论：IDFT 求解的矩阵即为 DFT 的逆矩阵，则我们可得：搬自YLY的课件

形式化地讲，令 $B (x) = I D FT (A (x))$ ，则有：
$b_k=\frac{1}{n}\sum_{i=0}a_i\omega_n^{-ik}=\frac{1}{n}A(\omega_n^{-k})$

由此，我们只需要将 DFT 实现过程中的 $\omega_n$ 修改为 $\omega_n^{-1}$ ，最后将所得结果乘上 $\frac{1}{n}$ 即可。

真正的快速傅里叶变换（FFT）

结合 DFT 与 IDFT，我们实现了 $O (n l o g n)$ 的多项式乘法。
然而，多次递归导致程序的实际效率并不高，~~甚至过不了模板题~~。所以，我们考虑迭代实现 FFT。

正难则反，如果我们可以正向分割问题处理，为何不可以倒着合并得出答案呢？
~~同样地，如果觉得抽象可以看图。~~
在这里插入图片描述
不妨令 $a$ 表示初始数组， $b$ 表示分治完后的数组， $rev_i$ 表示 $i$ 二进制位反转后所得到的数。
通过大眼观察法，可以得到： $a_i=b_{rev_i}$ ，因此，我们可以通过 $b$ 数组逆推得到答案。

NTT

某些题目要求我们在模意义下输出答案，这个时候，FFT便可能会产生较大的精度误差。为了解决这类模数问题，NTT 应运而生。

FFT 无法在计算中途取模，归咎于我们代入了复数意义下的单位根。因此，我们可以规避掉单位根，转而代入模意义下与单位根有着相似性质的数，这便是 NTT。
我们通常将这个数称为 原根( $g$ )。

注意到 FFT 处理的数组的长度一定是 $2$ 的整次幂，因此，我们的模数必须满足 $p=k2^x+1$ 的苛刻条件。
常见地， $998244353, 469762049, 1004535809$ 的原根是 $3$ 。

于是，我们只需要将 FFT 中的 $\omega_n$ 替换为原根即可。具体地，用 $g_n^1$ 满足 $[g_n^1 \equiv g^{\frac{p-1}{n}} \pmod{p} ]$ 替换 $\omega_n$ ，作逆运算的时候用逆元即可。

代码

给出我常用的 FFT，NTT模板。

#include<bits/stdc++.h>
#define db double 
// #define int long long
using namespace std;
const int MAXN = 1e6 + 5;
const db Pi = acos(-1.0);
namespace poly{
    const int MOD = 998244353;
    int to[MAXN] , G0 = 3 , G1 = 332748118 , wn[MAXN][2];
    int qpow(int base , int k) {
        int res = 1;
        while(k) {
            if (k & 1) res = res * base % MOD;
            base = base * base % MOD;
            k >>= 1;
        }
        return res;
    }
    int inv(int p) {return qpow(p , MOD - 2);}
    struct Complex{
        db x , y;
        Complex() {}
        Complex(db _x , db _y) {x = _x , y = _y;}
        friend Complex operator + (Complex a , Complex b) {return Complex(a.x + b.x , a.y + b.y);};
        friend Complex operator - (Complex a , Complex b) {return Complex(a.x - b.x , a.y - b.y);};
        friend Complex operator * (Complex a , Complex b) {return Complex(a.x * b.x - a.y * b.y , a.x * b.y + a.y * b.x);};
        friend Complex operator / (Complex a , db b) {return Complex(a.x / b , a.y / b);};
        friend Complex operator * (Complex a , db b) {return Complex(a.x * b , a.y * b);};
    }a[MAXN] , b[MAXN] , q1[MAXN] , q2[MAXN] , q3[MAXN];
    void init(int len) {
        int lg = __lg(len) + 1 ,limit = (1 << lg);
        for (int i = 0 ; i < limit ; i ++) to[i] = (to[i >> 1] >> 1) | ((i & 1) << (lg - 1));
        wn[lg][1] = qpow(G1 , (MOD - 1) / limit);
        wn[lg][0] = qpow(G0 , (MOD - 1) / limit);
        for (int i = lg ; i >= 1 ; i --) {
            wn[i - 1][0] = wn[i][0] * wn[i][0] % MOD;
            wn[i - 1][1] = wn[i][1] * wn[i][1] % MOD;
        }
    }
    void FFT(Complex *A , int up , int type){
        for (int i = 0 ; i < up ; i ++) {
            if (i > to[i]) swap(A[i] , A[to[i]]);
        }
        for (int mid = 1 ; mid < up ; mid <<= 1) {
            Complex Wn(cos(Pi / mid) , type * sin(Pi / mid));
            for (int len = mid << 1 , j = 0 ; j < up ; j += len) {
                Complex W(1 , 0);
                for (int i = 0 ; i < mid ; i ++ , W = W * Wn) {
                    Complex tmpx = A[j + i] , tmpy = W * A[j + mid + i];
                    A[j + i] = tmpx + tmpy , A[j + mid + i] = tmpx - tmpy;
                }
            }
        }
    }
    void NTT(int *A , int up , int type) {
        for (int i = 0 ; i < up ; i ++) {
            if (i > to[i]) swap(A[i] , A[to[i]]);
        }
        for (int mid = 1 , lg = 1 ; mid < up ; mid <<= 1 , lg ++) {
            int Wn = wn[lg][type];
            for (int len = mid << 1 , j = 0 ; j < up ; j += len) {
                int W = 1;
                for (int i = 0 ; i < mid ; i ++ , W = W * Wn % MOD) {
                    int tmpx = A[j + i] , tmpy = W * A[j + mid + i] % MOD;
                    A[j + i] = (tmpx + tmpy) % MOD , A[j + mid + i] = (tmpx - tmpy + MOD) % MOD;
                }
            }
        }
        if (!type) return;
        int inv = qpow(up , MOD - 2);   
        for (int i = 0 ; i < up ; i ++) A[i] = A[i] * inv % MOD;
    }
    void FFT_mul(int n , db *F , int m , db *G , db *res) {
        init(n + m);
        int up = __lg(n + m) + 1;
        up = (1 << up);
        for (int i = 0 ; i <= up ; i ++) a[i].x = a[i].y = b[i].x = b[i].y = 0;
        for (int i = 0 ; i <= n ; i ++) a[i].x = F[i];
        for (int i = 0 ; i <= m ; i ++) b[i].x = G[i];
        FFT(a , up , 1) , FFT(b , up , 1);
        for (int i = 0 ; i <= up ; i ++) a[i] = a[i] * b[i];
        FFT(a , up , -1); 
        for (int i = 0 ; i <= n + m ; i ++) res[i] = a[i].x;
    }
    void NTT_mul(int n , int *F , int m , int *G , int *res) {
        init(n + m);
        int up = __lg(n + m) + 1;
        up = (1 << up);
        NTT(F , up , 0) , NTT(G , up , 0);
        for (int i = 0 ; i < up ; i ++) F[i] = F[i] * G[i] % MOD;
        NTT(F , up , 1);
        for (int i = 0 ; i < up ; i ++) res[i] = F[i];
    }
}
int n , m;
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr) , cout.tie(nullptr);
    return 0;
}