FFT

前置知识：复数、原根、卷积。

对于 $n$ 次多项式

$$\begin{array}{rl}f(x)=\sum _{i=0}^n{f}_i{x}^i& =f_0+f_{1}x+f_2{x}^2+\cdots +g_n{x}^n\\ g(x)=\sum _{i=0}^n{g}_i{x}^i& =g_0+g_{1}x+g_2{x}^2+\cdots +g_n{x}^n\end{array}$$

它们的乘法是 $F(x)=f(x)g(x)=\sum _{k=0}^{2n}{c}_k{x}^k$，其中

$$c_k=\sum _{i+j=k}{f}_i{g}_j$$

因此计算多项式的乘积需要 $n^2$ 次系数乘法，我们需要优化。

点值表示法

$n$ 次多项式 $f(x)$ 可以由 $n+1$ 个系数决定，也可以由 $n+1$ 个座标（点值）决定。即 $n$ 次多项式可以看作 $n+1$ 维的向量。

考虑选取 $2n+1$ 个座标来确定 $f(x)$ 和 $g(x)$。则 $F(x)$ 可以简单的通过做 $2n+1$ 次乘法得到

$$(x_k,F(x_k))=(x_k,f(x_k)g(x_k))$$

于是求多项式的乘法，可以先从系数表示法转换为点值表示法，做完乘法再变回去。

单位原根

称方程 $x^n=1$ 的 $n$ 个解为单位根，其中特殊的一个记作 ${\omega }_n=e^{\frac{2\pi i}{n}}$，称为单位原根。根据 Euler 公式，有

$${\omega }_n=e^{\frac{2\pi i}{n}}=\cos \left(\frac{2\pi }{n}\right)+i\sin \left(\frac{2\pi }{n}\right)$$

即 ${\omega }_n$ 是单位圆上的一个点，全部的 $n$ 个单位根

$$x_k={\omega }_n^k=e^{k\frac{2\pi i}{n}}=\cos \left(\frac{2\pi k}{n}\right)+i\sin \left(\frac{2\pi k}{n}\right)$$

恰全是单位元的 $n$ 等分点。因此根据 Euler 公式，单位根之间的乘法就是在单位元上转圈圈。

DFT

设多项式 $f(x)=\sum _{k=0}^{n-1}{f}_k{x}^k$，并记 $\omega ={\omega }_n$，则称向量

$${\mathrm{DFT}}_{\omega }(f)=(f(1),f({\omega }^1),\cdots ,f({\omega }^{n-1}))$$

为 $f$ 的离散 Fourier 变换（DFT）。

DFT 存在逆变换（IDFT），即从点值重新变回系数，仍是从向量到向量的变换。

IDFT 存在性质

$$({\mathrm{DFT}}_{\omega }{)}^{-1}=\frac{1}{n}({\mathrm{DFT}}_{{\omega }^{-1}})$$

篇幅所限，不放证明。现在我们可以统一的处理 DFT 和 IDFT。

分治

利用单位元根的特殊性，我们可以分治计算 DFT。比如对于 $7$ 次多项式

$$\begin{array}{rl}f(x)& =f_0+f_{1}x+f_2{x}^2+f_3{x}^3+f_4{x}^4+f_5{x}^5+f_6{x}^6+f_7{x}^7\\ & =(f_0+f_2{x}^2+f_4{x}^4+f_6{x}^6)+x(f_1+f_3{x}^2+f_5{x}^4+f_7{x}^6)\end{array}$$

奇偶分别建立函数

$$\begin{array}{rl}{f}^{[0]}(x)& =f_0+f_{2}x+f_4{x}^2+f_6{x}^3\\ f^{[1]}(x)& =f_1+f_{3}x+f_5{x}^2+f_7{x}^3\end{array}$$

则原来的函数可以表示为

$$f(x)=f^{[0]}(x^2)+x{f}^{[1]}(x^2)$$

一般的，对于度小于 $n$ 的多项式 $f(x)$，利用单位原根的性质有

$$\begin{array}{rl}f({\omega }_n^k)& =f^{[0]}({\omega }_n^k\cdot {\omega }_n^k)+{\omega }_n^k{f}^{[1]}({\omega }_n^k\cdot {\omega }_n^k)\\ & =f^{[0]}({\omega }_n^{2k})+{\omega }_n^k{f}^{[1]}({\omega }_n^{2k})\\ & =f^{[0]}({\omega }_{n/2}^k)+{\omega }_n^k{f}^{[1]}({\omega }_{n/2}^k)\end{array}$$

同理可得

$$\begin{array}{rl}f({\omega }_n^{k+n/2})& =f^{[0]}({\omega }_n^{2k+n})+{\omega }_n^{k+n/2}{f}^{[1]}({\omega }_n^{2k+n})\\ & =f^{[0]}({\omega }_{n/2}^k)-{\omega }_n^k{f}^{[1]}({\omega }_{n/2}^k)\end{array}$$

因此我们需要把多项式的系数向上补到 $2^n$ 个，方便分治。

在 DFT 中使用有

$$\begin{array}{rl}{\mathrm{DFT}}_{\omega }(f)[j]& ={\mathrm{DFT}}_{{\omega }^2}(f^{[0]})[j]+{\omega }^j{\mathrm{DFT}}_{{\omega }^2}(f^{[1]})[j]\\ {\mathrm{DFT}}_{\omega }(f)[j+n/2]& ={\mathrm{DFT}}_{{\omega }^2}(f^{[0]})[j]-{\omega }^j{\mathrm{DFT}}_{{\omega }^2}(f^{[1]})[j]\end{array}$$

至此，我们可以写出递归版的 FFT。

void DFT(Comp *f, int n, int type) {
    if (n == 1)
        return;
    for (int i = 0; i < n; i++)
        tmp[i] = f[i];
    for (int i = 0; i < n; i++) {  // 偶数放左边，奇数放右边
        if (i & 1)
            f[n / 2 + i / 2] = tmp[i];
        else
            f[i / 2] = tmp[i];
    }
    Comp *g = f, *h = f + n / 2;
    DFT(g, n / 2, type), DFT(h, n / 2, type);
    Comp step(cos(2 * PI / n), sin(2 * PI * type / n));
    Comp cur(1, 0);

    for (int k = 0; k < n / 2; k++) {
        tmp[k] = g[k] + cur * h[k];
        tmp[k + n / 2] = g[k] - cur * h[k];
        cur = cur * step;
    }
    for (int i = 0; i < n; i++)
        f[i] = tmp[i];
}

void DFT(Comp *f, int n, int type) {
    if (n == 1)
        return;
    for (int i = 0; i < n; i++)
        tmp[i] = f[i];
    for (int i = 0; i < n; i++) {  // 偶数放左边，奇数放右边
        if (i & 1)
            f[n / 2 + i / 2] = tmp[i];
        else
            f[i / 2] = tmp[i];
    }
    Comp *g = f, *h = f + n / 2;
    DFT(g, n / 2, type), DFT(h, n / 2, type);
    Comp step(cos(2 * PI / n), sin(2 * PI * type / n));
    Comp cur(1, 0);

    for (int k = 0; k < n / 2; k++) {
        tmp[k] = g[k] + cur * h[k];
        tmp[k + n / 2] = g[k] - cur * h[k];
        cur = cur * step;
    }
    for (int i = 0; i < n; i++)
        f[i] = tmp[i];
}

蝴蝶变换

递归分治总是不尽人意的，如果能一次到位就更好了。还是以 $7$ 次多项式为例

• 初始 $\{x^0,x^1,x^2,x^3,x^4,x^5,x^6,x^7\}$
• 一次 $\{x^0,x^2,x^4,x^6\},\{x^1,x^3,x^5,x^7\}$
• 两次 $\{x^0,x^4\},\{x^2,x^6\},\{x^1,x^5\},\{x^3,x^7\}$
• 结束 $\{x^0\},\{x^4\},\{x^2\},\{x^6\},\{x^1\},\{x^5\},\{x^3\},\{x^7\}$

写出二进制的形式，可以发现规律

$$\begin{array}{cccccccccc}\text{初始}& 0& 1& 2& 3& 4& 5& 6& 7& \\ \text{初始(2)}& 000& 001& 010& 011& 100& 101& 110& 111& \\ \text{结束(2)}& 000& 100& 010& 110& 001& 101& 011& 111& \\ \text{结束}& 0& 4& 2& 6& 1& 5& 3& 7& \end{array}$$

结束和开始的二进制恰好是相反的。这个变换称为蝴蝶变换，也称位逆序置换（bit-reversal permutation）。

我们可以 $O(n)$ 的预处理出变换数组。设 R(x) 是 $x$ 的变换结果，则 R(x>>1) 是已求的。

$$\begin{array}{rl}\mathtt{\text{000abcd}}& \to \mathtt{\text{dcba000}}\\ \mathtt{\text{00abcdx}}& \to \mathtt{\text{xdcba00}}\end{array}$$

即是把 R(x>>1) 右移一位再补上最高位即可。代码如下

int lim = 1, lim_2;
while (lim <= n + m)
    lim <<= 1;
lim_2 = lim >> 1;
for (int i = 0; i < lim; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1)
        rev[i] |= lim >> 1;
    // 或者合并写为
    // rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * lim_2);
}

int lim = 1, lim_2;
while (lim <= n + m)
    lim <<= 1;
lim_2 = lim >> 1;
for (int i = 0; i < lim; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1)
        rev[i] |= lim >> 1;
    // 或者合并写为
    // rev[i] = (rev[i >> 1] >> 1) | ((i & 1) * lim_2);
}

现在我们可以写出非递归版的 FFT。

void FFT(Comp *f, int n, int type) {
    for (int i = 0; i < n; ++i) {
        if (i < rev[i]) {
            swap(f[i], f[rev[i]]);
        }
    }
    for (int h = 2; h <= n; h <<= 1) {
        Comp step(cos(2 * PI / h), sin(2 * PI * type / h));
        for (int j = 0; j < n; j += h) {
            Comp cur(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                Comp f1 = f[k], f2 = f[k + h / 2];
                f[k] = f1 + cur * f2;
                f[k + h / 2] = f1 - cur * f2;
                cur = cur * step;
            }
        }
    }
    if (type == 1)
        return;
    for (int i = 0; i < n; i++)
        f[i].x /= n, f[i].y /= n;
}

void FFT(Comp *f, int n, int type) {
    for (int i = 0; i < n; ++i) {
        if (i < rev[i]) {
            swap(f[i], f[rev[i]]);
        }
    }
    for (int h = 2; h <= n; h <<= 1) {
        Comp step(cos(2 * PI / h), sin(2 * PI * type / h));
        for (int j = 0; j < n; j += h) {
            Comp cur(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                Comp f1 = f[k], f2 = f[k + h / 2];
                f[k] = f1 + cur * f2;
                f[k + h / 2] = f1 - cur * f2;
                cur = cur * step;
            }
        }
    }
    if (type == 1)
        return;
    for (int i = 0; i < n; i++)
        f[i].x /= n, f[i].y /= n;
}

三步并做两步

实战一下：P3803 多项式乘法。

我们的计算步骤是：

FFT(F, lim, 1);
FFT(G, lim, 1);
for (int i = 0; i <= lim; i++)
    F[i] = F[i] * G[i];
FFT(F, lim, -1);

FFT(F, lim, 1);
FFT(G, lim, 1);
for (int i = 0; i <= lim; i++)
    F[i] = F[i] * G[i];
FFT(F, lim, -1);

实际上，我们并不用三次 FFT，两次足以。注意到若把 $F(x)$ 放在实部而 $G(x)$ 放在虚部

$$(F+iG{)}^2=(F^2-G^2)+2iFG$$

平方之后虚部恰是答案。

无 REV

实际上 REV 是可以消掉的，完全的推导比较复杂，这里只是以板子的形式给出。

template <class iter>
void ntt(iter f, int n) {
    pre_w(n), ntt_size += n;
    for (int l = n / 2; l; l >>= 1)
        for (int i = 0; i < n; i += l * 2)
            for (int j = 0; j < l; j++) {
                int x = f[i + j], y = f[i + j + l];
                f[i + j + l] = 1ll * (x - y + P) * w[j + l] % P;
                f[i + j] = mo(x + y);
            }
}

template <class iter>
void intt(iter f, int n) {
    pre_w(n), ntt_size += n;
    for (int l = 1; l < n; l <<= 1)
        for (int i = 0; i < n; i += l * 2)
            for (int j = 0; j < l; j++) {
                int x = f[i + j];
                int y = 1ll * w[j + l] * f[i + j + l] % P;
                f[i + j] = mo(x + y), f[i + j + l] = mo(x - y + P);
            }
    const int iv = qpow(n);
    for (int i = 0; i < n; i++)
        f[i] = 1ll * f[i] * iv % P;
    reverse(f + 1, f + n);
}

template <class iter>
void ntt(iter f, int n) {
    pre_w(n), ntt_size += n;
    for (int l = n / 2; l; l >>= 1)
        for (int i = 0; i < n; i += l * 2)
            for (int j = 0; j < l; j++) {
                int x = f[i + j], y = f[i + j + l];
                f[i + j + l] = 1ll * (x - y + P) * w[j + l] % P;
                f[i + j] = mo(x + y);
            }
}

template <class iter>
void intt(iter f, int n) {
    pre_w(n), ntt_size += n;
    for (int l = 1; l < n; l <<= 1)
        for (int i = 0; i < n; i += l * 2)
            for (int j = 0; j < l; j++) {
                int x = f[i + j];
                int y = 1ll * w[j + l] * f[i + j + l] % P;
                f[i + j] = mo(x + y), f[i + j + l] = mo(x - y + P);
            }
    const int iv = qpow(n);
    for (int i = 0; i < n; i++)
        f[i] = 1ll * f[i] * iv % P;
    reverse(f + 1, f + n);
}

应用

FFT 是卷积计算的基础工具。更多应用见 卷积。

练习

• P3803 多项式乘法

参考资料

• 快速傅里叶变换 - OI Wiki