「快速傅里叶变换」总结

前言

作为 NOI 大纲里的十级算法，FFT拥有着逼格十分高的名称，也十分令人头疼？（雾
通俗地讲，FFT旨在以优异的时间复杂度，解决两个多项式相乘的问题，即：
$$c_k=\sum _{i=0}^k{a}_i{b}_{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_{1}x+a_2{x}^2+a_3{x}^3+\cdot \cdot \cdot \cdot \cdot \cdot +a_{n-1}{x}^{n-1}$。
将其按照下标奇偶分类，可得：
$$A_1(x)=a_0+a_{2}x+a_4{x}^2+\cdot \cdot \cdot \cdot \cdot \cdot +a_{n-2}{x}^{\frac{n}{2}}$$$$A_2(x)=a_1+a_{3}x+a_5{x}^2+\cdot \cdot \cdot \cdot \cdot \cdot +a_{n-1}{x}^{\frac{n}{2}}$$
显然有，$A(x)=A_1(x^2)+x{A}_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(nlogn)$ 的时间复杂度用 点值表示法 表示多项式 $A(x)$！

快速傅里叶逆变换（IDFT）

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

结论：IDFT 求解的矩阵即为 DFT 的逆矩阵，则我们可得：

形式化地讲，令 $B(x)=IDFT(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(nlogn)$ 的多项式乘法。
然而，多次递归导致程序的实际效率并不高，甚至过不了模板题。所以，我们考虑迭代实现 FFT。

正难则反，如果我们可以正向分割问题处理，为何不可以倒着合并得出答案呢？
同样地，如果觉得抽象可以看图。

不妨令 $a$ 表示初始数组，$b$ 表示分治完后的数组，$re{v}_i$ 表示 $i$ 二进制位反转后所得到的数。
通过大眼观察法，可以得到：$a_i=b_{re{v}_i}$，因此，我们可以通过 $b$ 数组逆推得到答案。

NTT

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

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

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

于是，我们只需要将 FFT 中的 ${\omega }_n$ 替换为原根即可。具体地，用 $g_n^1$ 满足 $[g_n^1\equiv g^{\frac{p-1}{n}}(\mathrm{mod}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;
}