快速傅里叶变换

FFT 原理

令 \(a, b \in \mathbb R ^ {n}\)，记 \(a\) 和 \(b\) 循环卷积的结果 \(c\)：

\[c_k = \sum\nolimits_{i + j \bmod n = k} a_i b_j \]

我们想要的变换是：

\[M(Aa \odot Bb) = c \]

我们想要好一点的性质，如交换和结合，我们令 \(A = B = M ^ {-1} = D\)，则有：

\[Da \odot Db = Dc \]

此时，考虑 \(e_i\) 取值，有 \(D_{k, i}D_{k, j} = D_{k, i + j \bmod n}\)，解得 \(D_{k, i} = D_{k, 1} ^ i, D_{k, 1} ^ n = 1\)，为了满秩，令 \(D_{k, 1} = w_n ^ k\)，则得到了范德蒙德矩阵：

\[D = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & w_n & w_n^2 & \cdots & w_n^{n-1} \\ 1 & w_n^2 & w_n^4 & \cdots & w_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_n^{n-1} & w_n^{2(n-1)} & \cdots & w_n^{(n - 1)(n-1)} \end{bmatrix}\]

再考虑 \(D ^ {-1}(Da \odot Db) = c\) 在 \(e_i\) 处的点值：

\[\sum_{p = 1} ^ r D ^ {-1}_{i + j \bmod n, p} D_{p, i} D_{p, j} = \sum_{p = 1} ^ r D ^ {-1}_{k, p} w_n ^ {kp} = 1 \]

容易发现 \(D^{-1}_{i, j} = \dfrac{1}{n}w_n ^ {-ij}\)，而实践中更常常实现为 \(D^{-1} = \dfrac{1}{n} \operatorname{perm}(0, n - 1, \cdots, 2, 1)D\)，其中：

\[\operatorname{perm}(0, n - 1, \cdots, 2, 1) := \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 1 \\ 0 & 0 & 0 & \cdots & 1 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 1 & \cdots & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 \end{bmatrix}\]

在实现时多项式乘法时，计算出 \(Da, Db\)，将二者点积得到 \(Da \odot Db\)，再计算 \(D(Da \odot Db)\)，然后翻转区间 \([1, n - 1]\)，最后整体除以 \(n\) 即可。

FFT

我们需要的是快速计算 \(Da\)，为了方便，我们设 \(n\) 是 \(2\) 的幂。

展开式子，就是给定 \(x_i\)，求：

\[X_k = \sum_{i = 0} ^ {n - 1} x_i w_n ^ {ik} \]

称快速进行这种操作的算法为 FFT。

DIT - FFT

考虑分治，设 \(y_r = x_{2r}, z_r = x_{2r + 1}, m = n / 2\)，则有：

\[\begin{aligned} X_k = \sum_{i = 0} ^ {n - 1} x_i w_n ^ {ik} & = \sum_{i = 0} ^ {m - 1} y_i w_m ^ {ik} + w_n ^ k \sum_{i = 0} ^ {m - 1} z_i w_m ^ {ik} = Y_k + w_n ^ k Z_k \\ X_{k + m} = \sum_{i = 0} ^ {n - 1} x_i w_n ^ {i(k + m)} & = \sum_{i = 0} ^ {m - 1} y_i w_m^{ik} - w_n ^ k \sum_{i = 0} ^ {m - 1} z_i w_m ^ {ik} = Y_k - w_n ^ k Z_k \end{aligned}\]

有时间复杂度 \(T(n) = O(n) + 2T(n / 2) = O(n \log n)\)

注意到，初始时直接将序列排成左半边是 \(y\)，右半边是 \(z\) 的形式，则只需在原地做变换，这种变换称为蝴蝶变换。

记蝴蝶变换的矩阵为 \(P\)，DIT-FFT 对应的矩阵为 \(Q\)，则有：

\[QP = D \]

每次变换的系数矩阵为：\(\displaystyle \begin{bmatrix} 1 & w_n ^ k \\ 1 & -w_n ^ k \end{bmatrix}\)

DIF - FFT

注意到，\(D, P\) 都是自逆的，故 \(D = QP = PQ^\top\)

发现做蝴蝶变换矩阵 \(P ^ 2 = I\)，有 \(D ^ 2 = Q Q^\top\)，可以节省掉蝴蝶变换！

由于 \((AB)^\top = B^\top A^\top\)，故则其应当从大区间至小区间枚举。

每次变换的系数矩阵为：\(\displaystyle \begin{bmatrix} 1 & 1 \\ w_n ^ k & -w_n ^ k \end{bmatrix}\)。

这称为 转置原理。

若将计算过程画成一张图，以边权表示乘的系数，则如下：

代码

void InitOmega(int N) {
    static int tb[23]; long long rt = 31;
    for(int i = 22; i >= 0; --i)
        tb[i] = (int)rt, rt = rt * rt % Mod;
    for(int L = 1, p = 0; L < N; L <<= 1) {
        long long t = 1; const int q = tb[p++];
        for(int i = 0; i < L; ++i)
            w[L + i] = (int)t, t = t * q % Mod;
    }
}
void NTT(int *A, int N) {
    const auto &trans = [&](int i, int j, int L) {
        const int u = A[i + j], v = A[i + j + L];
        A[i + j] = u + v >= Mod ? u + v - Mod : u + v;
        A[i + j + L] = 1ll * (u + Mod - v) * w[L + j] % Mod;
    };
    for(int L = N >> 1; L; L >>= 1)
        for(int i = 0; i < N; i += (L << 1))
            for(int j = 0; j < L; ++j) trans(i, j, L);
}
void INTT(int *A, int N) {
    const auto &trans = [&](int i, int j, int L) {
        const int u = A[i + j], v = 1ll * w[L + j] * A[i + j + L] % Mod;
        A[i + j] = u + v >= Mod ? u + v - Mod : u + v;
        A[i + j + L] = u - v >= 0 ? u - v : u + Mod - v;
    };
    for(int L = 1; L < N; L <<= 1)
        for(int i = 0; i < N; i += (L << 1))
            for(int j = 0; j < L; ++j) trans(i, j, L);
    reverse(A + 1, A + N); const int iv = Mod - (Mod - 1) / N;
    for(int i = 0; i < N; ++i) A[i] = (int)(1ll * A[i] * iv % Mod);
}