FFT板子题

# include <bits/stdc++.h>
using namespace std;
const int N = 3e6;
const double PI = acos(-1);
struct node {
    double x, y;
    node operator + (const node& t) const {
        return {x + t.x, y + t.y};
    }
    node operator - (const node& t) const {
        return {x - t.x, y - t.y};
    }
    node operator * (const node& t) const {
        return {x * t.x - y * t.y, x * t.y + y * t.x};
    }
};
node A[N], B[N];
char s1[N], s2[N];
int R[N], ans[N];
void FFT(node A[], int n, int op) {
    for (int i = 0; i < n; i++) {
        R[i] = R[i / 2] / 2 + ((i & 1) ? n / 2 : 0);
    }
    for (int i = 0; i < n; i++) {
        if (i < R[i]) {
            swap(A[i], A[R[i]]);
        }
    }
    for (int i = 2; i <= n; i <<= 1) {
        node w1 = {cos(2 * PI / i), sin(2 * PI / i) * op};
        for (int j = 0; j < n; j += i) {
            node wk = {1, 0};
            for (int k = j; k < j + i / 2; k++) {
                node x = A[k], y = A[k + i / 2] * wk;
                A[k] = x + y; A[k + i / 2] = x - y;
                wk = wk * w1;
            }
        }
    }
}
int main() {
    scanf("%s%s", s1, s2);
    int n = strlen(s1) - 1, m = strlen(s2) - 1;
    for (int i = 0; i <= n; i++) {
        A[i].x = s1[n - i] - '0';
    }
    for (int i = 0; i <= m; i++) {
        B[i].x = s2[m - i] - '0';
    }
    for (m = n + m, n = 1 ;n <= m; n <<= 1);
    FFT(A, n, 1);
    FFT(B, n, 1);
    for (int i = 0; i < n; i++) {
        A[i] = A[i] * B[i];
    }
    FFT(A, n, -1);
    int k = 0;
    for (int i = 0, t = 0; i < n || t; i++) {
        t += A[i].x / n + 0.5;
        ans[k++] = t % 10;
        t /= 10;
    }
    while (k > 1 && !ans[k - 1]) {
        k--;
    }
    for (int i = k - 1; i >= 0; i--) {
        printf("%d", ans[i]);
    }
    return 0;
}

编写它的目的是备忘。

---

前置知识

• 我们把（最多）n次多项式 $F(x)$ 表示成 $\sum_{d=0}^n{a_dx^d}$ 的方法称作系数表示法。
• 我们把（最多）n次多项式 $F(x)$ 用 $F(x_0),F(x_1),F(x_2),\dots F(x_n)$ 表示的方法叫做点值表示法，点值表示和系数表示可以相互转化。
• 基本的任意角三角函数知识(即自变量的取值范围是任意实数）。
• 基本的向量知识：了解坐标向量的和，差，数乘，向量积即可。
• 基本的复数知识：$a+bi$。a，b是实数，$i^2=-1$。
• 基本的向量表示复数知识：我们把 $a+bi$ 看作 $(a,b)$ 向量，将这个向量的模长看作复数的模长（可见模长是非负实数），将这个向量与x轴正半轴所成的角叫做幅角（在x轴上为正，在下为负，可加减2π），复数相乘，模长相乘，幅角相加。
• 基本的反演知识，推荐炫酷反演魔术（THE XUANKU INVERSION MAGIC），如果你希望自己推出来加深记忆的话。
• 基本的求模知识：我们把 $a$ 除以 $b$ 的余数记作 $a\bmod b$ ,将 $a$ 除以 $c$ 的余数与 $b$ 除以 $c$ 的余数相同记作 $a=b\pmod c$。
• 一个证明技巧：若一组数 $x_1,x_2,\dots x_n$ 互不相等，一组数 $y_1,y_2,\dots y_n$ 互不相等，且对任意 $x_i$ 存在 $y_j=x_i$ ，则 $x_1,x_2,\dots x_n$ 与 $y_1,y_2,\dots y_n$ 一一对应，由此可推出和，乘积和一些奇怪的东西相等。

---

简述基本过程

我依稀记得我到最后才理解它大概操作的悲惨故事，因此我首先要讲它为什么可以加速多项式乘法。

• 对于一个 n 次多项式 $F(x)$ 和一个 m 次多项式 $G(x)$，乘积最多不过n+m次，因此如果可以迅速将 $F(x)$,$G(x)$ 转换为点值表示法($x_1,x_2,...,x_n$ 相同以便于乘起来)，再用O(n)求其乘积 $H(x)$ 的点值表示，再快速把 $H(x)$ 转为系数表示，就可以了。
• O(n)进行点值操作基本原理： $F(x_0)G(x_0),F(x_1)G(x_1),F(x_2)G(x_2),\dots F(x_n)G(x_n)$便是$H(x)$的点值表示。
• 唯一的问题是：普通的系数转点值是一个个代入求值，需要 O(n2），系数转点值更是涉及到高斯消元，O(n3)，我会拉格朗日插值O(n2） 。
• 不过，我们发现$x_0,x_1,x_2,...,x_n$的选取几乎没有任何要求，只要它们两两不等即可。
• 基于这个思路，我们有了 FFT，它通过选取极度特殊的点，使得点值系数互相转化的复杂度成为 $O(n\log n)$。
• 这个特殊的点值是什么呢，下面会提到：本原单位根，很佛系的名字。

---

本原单位根

• 定义：$x^n=1$ 在复数域上的所有根，称作 $n$ 次本原单位根。
• 根据模长相乘幅角相加，又因模长为非负实数，所以这个数模长为1，而幅角，则一定是 $\frac{2\pi}n$ 的整数倍（排除重复的，有n个互不相等的复数根）。
• 模长为1，幅角为 $\frac{2\pi k}n$ 的一个数（k是整数），计作 $\omega_n^k$，它是 $x^n=1$ 的一个根。
• 我们当然有必要疯狂探寻它的性质。
• 明显数值上 $\omega_n^k=\cos(\frac{2\pi k}n)+i\sin(\frac{2\pi k}n)$
• $\omega_n^k=\omega_n^{k+n},\omega_n^k=\omega_{nm}^{km}$ 和 $\omega_n^k=(\omega_n^1)^k$ 也是很明显的。
• 此外 $\omega_n^0=1,\omega_{2n}^n=-1$ 因此有 $\omega_{2n}^{n+k}=-\omega_{2n}^k$ 。
• 尤为重要的：$\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$ 互不相等，这意味着它可以作为点值表示的自变量。

---

FFT是什么？

FFT 是一种以分治算法为基础，运用本原单位根的特点进行反演来用。$O(n\log n)$进行点值和系数的互相转化的方法，它包含系数转点值(DFT)，点值转系数（IDFT）。

为(分治）方便起见，下面提到的n是二的正整数次幂。

---

DFT的具体做法

手抓一个萌萌的系数表示的多项式：$F(x)=\sum_{d=0}^{n-1}{a_dx^d}$

（求和号顶上的那个必须是$n-1$，这样点值和系数个数都是二的整数次幂，利于分治）

按次数的奇偶性分治：

设

$$\begin{cases} F_0(x)=\sum_{d=0}^{(n-1)/2}{a_{2d}x^d}\\ F_1(x)=\sum_{d=0}^{(n-1)/2}{a_{2d+1}x^d} \end{cases}$$

则

$$F(x)=F_0(x^2)+xF_1(x^2)$$

对于 $\omega_n^0,\omega_n^1,...,\omega_n^{n-1}$，在这里把它分成两类：$\omega_n^k,\omega_n^{n/2+k}$($k<\frac{n}2$)

$$F(\omega_n^k)=F_0((\omega_n^k)^2)+\omega_n^kF_1((\omega_n^k)^2)=F_0(\omega_n^{2k})+\omega_n^kF_1(\omega_n^{2k}) $$$$=F_0(\omega_{n/2}^k)+\omega_n^kF_1(\omega_{n/2}^k) $$$$F(\omega_n^{n/2+k})=F_0((\omega_n^{n/2+k})^2)+\omega_n^{n/2+k}F_1((\omega_n^{n/2+k})^2) $$$$=F_0(\omega_n^{n+2k})-\omega_n^kF_1(\omega_n^{n+2k})=F_0(\omega_{n/2}^k)-\omega_n^kF_1(\omega_{n/2}^k) $$

因此我们只要处理总共n个 $F_1(x),F_2(x)$ 的值即可O(n)求出 $F(x)$ 的n个值。

递归下去有logn层。时间复杂度O(nlogn)。

---

IDFT的大概做法

• 似乎有人用逆矩阵搞出来了，但我才初中啊。

反演是怎么看出来的：

我们已经知道了通过 $a_0,a_1,\dots,a_{n-1}$ 可表示 $F(\omega_n^0),F(\omega_n^1),\dots F(\omega_n^{n-1})$，想要通过 $F(\omega_n^0),F(\omega_n^1),\dots F(\omega_n^{n-1})$ 求 $a_0,a_1,\dots,a_{n-1}$。这不是反演是什么？

大概流程：

• 一个先要条件：

$$F(\omega_n^k)=\sum_{d=0}^{n-1}{a_d\omega_n^{dk}}$$

• 一句大白话：（$[a]$是当a为真时值为1，a为假时值为0。）

$$a_k=\sum_{d=0}^{n-1}{[d=k]a_d}$$

不过，更好的形式是：

$$a_k=\sum_{d=0}^{n-1}{[(k-d)\bmod n=0]a_d}$$

• 一个等式：萌萌哒=$[(k-d)\bmod n=0]$。
• 接下来，我们就可以反演了。

不过，我们需要寻找萌萌哒（一个与题面密切相关的式子）。

---

寻找萌萌哒

寻找萌萌哒首先需要探寻更多的本原单位根性质。

考虑 $\sum_{d=0}^{n-1}\omega_n^d$ :

$$\sum_{d=0}^{n-1}\omega_n^d=\sum_{d=0}^{n-1}(\omega_n^1)^d=\frac1{\omega_n^1-1}(\omega_n^1-1)\sum_{d=0}^{n-1}(\omega_n^1)^d $$$$=\frac1{\omega_n^1-1}(\sum_{d=1}^{n}(\omega_n^1)^d-\sum_{d=0}^{n-1}(\omega_n^1)^d)=\frac1{\omega_n^1-1}((\omega_n^1)^n-1)=\frac{1-1}{\omega_n^1-1}=0 $$

注意，当n=1时，$\omega_n^1=1$ ，此时变换不成立。$\sum_{d=0}^{n-1}\omega_n^d=1$

综上：

$$\sum_{d=0}^{n-1}\omega_n^d=[n=1]$$

考虑高级一点的式子：

$$G(n,m)=\sum_{d=0}^{n-1}\omega_n^{dm}$$

还记得 $\omega_n^k=\omega_n^{k+n}$吗

这让我们把它和模运算联系起来：$\omega_n^k=\omega_n^{k\mod n}$ ,这对下面很重要。

考虑两种情况(默认$\gcd(a,0)=a$)：

• $\gcd(n,m)=1$

此时 若 $d_1\neq d_2\pmod n$ ，则 $md_1\neq md_2\pmod n$

也就是 $0,m,2m,3m,\dots(n-1)m \mod n$ 两两不等。

也就是说 $0,m,2m,3m,\dots(n-1)m \mod n$ 与 $0,1,2,...,n-1$ 一一对应。

$$G(n,m)=\sum_{d=0}^{n-1}\omega_n^d=[n=1]$$

• $\gcd(n,m)\neq1$

直觉上：$0,m,2m,3m,\dots(n-1)m \mod n$ 与 $0,\gcd(n,m),2\gcd(n,m),...,n-\gcd(n,m)$有某种对应关系。

对于$0,1,\dots n-1$中的一个数 $d$ ,它可以表示成

$$p\frac{n}{\gcd(n,m)}+q(p<\gcd(n,m),q<\frac{n}{\gcd(n,m)}) $$

则

$$dm=pn\frac{m}{\gcd(n,m)}+qm=qm\pmod n$$

因此如果 $d_1=d_2\pmod{\frac{n}{\gcd(n,m)}}$ ，则$md_1=md_2\pmod n$

那么如果 $d_1\neq d_2\pmod{\frac{n}{\gcd(n,m)}}$ 呢

由于 $\gcd(\frac n{\gcd(n,m)},m)=1$：

$$md_1\neq md_2\pmod{ \frac n{\gcd(n,m)}}$$

即

$$md_1\neq md_2\pmod n$$

也就是说 $0,m,2m,3m,\dots(\frac n{\gcd(n,m)}-1)m \bmod n$两两不等

也就是说 $0,m,2m,3m,\dots(\frac n{\gcd(n,m)}-1)m \bmod n$ 与 $0,\gcd(n,m),2\gcd(n,m),...,n-\gcd(n,m)$ 一一对应。

也就是说 $0,m,2m,3m,\dots(n-1)m \mod n$ 与 $\gcd(n,m)$组 $0,\gcd(n,m),2\gcd(n,m),...,n-\gcd(n,m)$ 不计顺序是完全相同的。

因此可进行如下变形。

$$G(n,m)=\gcd(n,m)\sum_{d=0}^{\frac n{\gcd(n,m)}-1}\omega_n^{d\cdot\gcd(n,m)}=\gcd(n,m)\sum_{d=0}^{\frac n{\gcd(n,m)}-1}\omega_{\frac n{\gcd(n,m)}}^d $$$$=\gcd(n,m)[\frac n{\gcd(n,m)}=1]=\gcd(n,m)[m\bmod n=0] $$

即当$m\bmod n=0$时，$G(n,m)=n$ ，否则$G(n,m)=0$

• 综上， $G(n,m)=n[m\bmod n=0]$

也就是说，萌萌哒就是：

$$\frac{G(n,m)}n=[m\bmod n=0]=\frac1n\sum_{d=0}^{n-1}\omega_n^{dm} $$

所有前置准备已经完了，接下来我们开始反演。

---

IDFT是怎样炼成的

• 一个先要条件：

$$F(\omega_n^k)=\sum_{d=0}^{n-1}{a_d\omega_n^{dk}}$$

• 一句大白话：

$$a_k=\sum_{d=0}^{n-1}{[(k-d)\mod n=0]a_d}$$

• 一个等式：

$$[m\bmod n=0]=\frac1n\sum_{d=0}^{n-1}\omega_n^{dm}$$

• 代入变形：

$$a_k=\sum_{d=0}^{n-1}{(\frac1n\sum_{m=0}^{n-1}\omega_n^{m(k-d)})a_d}=\frac1n\sum_{d=0}^{n-1}{(\sum_{m=0}^{n-1}\omega_n^{m(k-d)})a_d} $$$$=\frac1n\sum_{d=0}^{n-1}{(\omega_n^{-md}\sum_{m=0}^{n-1}\omega_n^{mk})a_d}=\frac1n\sum_{m=0}^{n-1}{(\omega_n^{mk}\sum_{d=0}^{n-1}a_d\omega_n^{d(-m)})} $$

• 捕捉小朋友：

看看这一坨：

$$\sum_{d=0}^{n-1}a_d\omega_n^{d(-m)}$$

再看看先要条件：

$$F(\omega_n^k)=\sum_{d=0}^{n-1}{a_d\omega_n^{dk}}$$

是它，就是它，我们的朋友——多项式！

$$a_k=\frac1n\sum_{m=0}^{n-1}{(\omega_n^{mk}\sum_{d=0}^{n-1}a_d\omega_n^{d(-m)})}=\frac1n\sum_{m=0}^{n-1}{\omega_n^{mk}F(\omega_n^{-m})} $$

• 还看不出来怎么求吗？

用 $A(\omega_n^k)$ 表示 $a_k$ ，用 $f_k$ 表示 $\frac{F(\omega_n^{-k})}n$ ，这个式子长这样

$$A(x)=\sum_{d=0}^{n-1}f_dx^d$$

这不就巧了吗，它其实是 DFT 。因此 IDFT 的复杂度也是$O(n\log n)$

理论部分在此就告一段落了。

---

可怕的事情

完全按照我说的打的程序长这样：

改了下，在云剪贴板。

然后就毫无疑问地超时了。我需要优化。

---

优化1

我们每次都要新开一个数组，非常耗空间，我们看看每次变换时，各个系数到底去哪了（以下是$n=8$的情况）。

我们可以看到，对于第$4a+2b+c(0\leq a,b,c<1)$个数,它去到了$a+2b+4c$的位置(相当于被二进制翻转了),我们可以编一个数组来存下每次对应的位置变换，再向上合并即可。 这个数组的具体推法是什么呢？

$$\begin{cases} op(x)=0&x=0\\ op(x)=op(x/2)/2&2\mid x\\ op(x)=op((x-1)/2)/2+n/2&2\nmid x \end{cases}$$

然后空间就省到O(n)了。

---

优化2

引自NaCly_Fish。

• 我们可以将3次FFT变成2次FFT。
• 对于一个n次多项式$F(x)$和一个m次多项式$G(x)$，设$H(x)=F(x)+iG(x)$，则$\frac{H^2(x)}2$的虚部即为答案。
• 理由：

$$\frac{H^2(x)}2=\frac{(F(x)+iG(x))^2}2=\frac{F^2(x)-G^2(x)+2F(x)G(x)i}2 $$$$=\frac{F^2(x)-G^2(x)}2+F(x)G(x)i$$

• 但似乎精度会更低。
• 代码实现。

---

总结

最后再看一看$\omega_n^k$的特殊性质。

• $\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}$两两不等。
• $\omega_n^k=(\omega_n^1)^k$
• $\omega_{2n}^k+\omega_{2n}^{n+k}=0$
• $\omega_n^k=\omega_n^{k\bmod n}$
• $\omega_n^0=1$

实际上，是这些性质保证了它可以进行快速变换(其它性质可通过这些证明）。

问题来了，除了本原单位根，难道就没有其他东西有这样的性质吗？

据说这种东西叫原根。

• 另外，数学家们已经证明，在复数域内，只有本原单位根满足这些性质。

---

原根

• 后来发现各大题解其实讲得很清楚，不知道我当时为什么不知道原根怎么用。
• 对于一个数 $m$ ，在 $1,2,\dots,m-1$ 中与 $m$ 互质的数的个数记作 $\varphi(m)$ 。
• 在这 $\varphi(m)$ 个数中，容易证明 $a^{\varphi(m)}=1$ （可搜索词条欧拉定理），那么如果 $a^b=1$ 在 $0<b<\varphi(m)$ 时永远不成立，则 $a$ 是模 $m$ 意义下的一个原根。
• 这里不再详细讲求法了。

---

998244353的原根

• 首先998244353是质数，满足对任意 $0<a<998244353$ 有 $a^{998244352}=1\pmod {998244353}$
• 其次，$998244353=7\times17\times2^{23}+1$
• 最后，3 是模 998244353 意义下的一个原根，也就是说：$3^0,3,\dots3^{998244351}$ 模998244353互不相同
• 这足以激发人的一点灵感，设 $\omega_n^k=3^{\frac{998244352k}n}\pmod {998244353}$ （注意n在这里是二的整数次幂且$\log_2 n\leq23$）

则$\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}\pmod {998244353}$两两不等。

$\omega_n^k=(\omega_n^1)^k\pmod {998244353}$

$\omega_{n}^k=-\omega_{n}^{n/2+k}\pmod {998244353}$

因为$3^{499122176}=-1\pmod {998244353}$

$\omega_n^k=\omega_n^{k\bmod n}\pmod {998244353}$

$\omega_n^0=1\pmod {998244353}$

• 就这么搞完了。
• 代码实现。
• 实测两次 FFT 更快。

吐槽与总结

• 吐槽：以前初三的代码好锅啊！而且只做完模板题就溜，根本不给我留一点后路啊！
• 从这篇中，你们可以看到一类卷积问题广泛的应用前景。
• 而我们历经辛苦，其实就是推出了这么个式子：

$$\omega_n^k=\cos(\frac{2\pi k}n)+i\sin(\frac{2\pi k}n) $$$$a_k=\frac1n\sum_{m=0}^{n-1}{(\omega_n^{mk}\sum_{d=0}^{n-1}a_d\omega_n^{d(-m)})}=\frac1n\sum_{m=0}^{n-1}{\omega_n^{mk}F(\omega_n^{-m})} $$

• 它说明逆变换就是把正变换后的系数反过来然后结果除个 $n$。
• 然后我们利用分治兜兜转转推出了 $O(n\log n)$ 的变换方式来计算卷积。
• 我们辛苦自主推出的结论其实不过是范德蒙德矩阵的逆，出于实用的角度谴责这一行为。
• 同样，出于实用的角度，本题的代码有大幅度优化空间且例题很少，这里会补充部分例题。
• 计算卷积长度：

L=1<<(__lg(n+m)+1);

• 蝴蝶变换转化:

for(int i=1;i<L;++i)
		sw[i]=i&1?sw[i>>1]:(sw[i>>1]|(L>>1));

• 正变换：
• 逆变换：
• FFT，NTT。
• 现在想来，十分简单，但是如果你其它数学不过关，这些也只是模板而已，是不会发挥其应有的效果的。

例题 1

• 题目。
• 先做个除法约分，然后问题就变成卷积的形式了，多么简单啊……（靠这话你也说的出口）
• 代码，多测不清空，爆零两行泪。

例题 2

• 既然标题都起成这样了，考虑分治。
• 假设已经求出了 $0$ 到 $n/2-1$ 考虑对右边的贡献其实也是个卷积的形式，可以 $O(n\log n)$ 暴力算贡献。
• 最终容易得到复杂度为 $O(n\log^2n)$，代码。

Ruby yyds

a=gets.to_i
b=gets.to_i
print a*b