FFT可爱QwQ！！！
\(\text{upd 2023.11.23:}\)
发现自己 FFT 老忘，所以决定在文章开头弄个速通版 FFT 讲解。
给定多项式 \(A,B\) 的系数表示，求 \(C=A\times B\) 的系数表示，FFT 大致操作：
分别求 \(A,B\) 的点值表示，求出来之后将对应位相乘得到 \(C\) 的点值表示：
- 拆分多项式：\(C(x)=C_1(x^2)+x\cdot C(x^2),C_1(x)=x_0+x_1\cdot x+\dots,C_2(x)=x_1+x3\cdot x+\dots\)，对于 \(C{1⁄2}(x)\) 的系数拆分，可运用二进制位翻转来直接进行位置互换。
贴个二进制位翻转的代码。
int sum_len = len_a+len_b, len = 0;
while (n <= sum_len)
++len, n <<= 1;
for (int i = 1; i < n; i++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
for (int i = 0; i < n; i++)
if (i < r[i]) swap(a[i], a[r[i]]);
* 然后再运用式子进行\\(C(x)\\) 的离散傅里叶变换求值：
* \\(i\in[0,\frac{n}{2}):C(\omega_{n}^{i})=C_1(\omega_{\frac{n}{2}}^{i})+\omega_{n}^{i}\cdot C_2(\omega_{\frac{n}{2}}^{i})\\)。
* \\(i\in[\frac n2,n-1]:C(\omega_{n}^{i})=C_1(\omega_{\frac{n}{2}}^{i-\frac n2})-\omega_{n}^{i-\frac n2}\cdot C_2(\omega_{\frac{n}{2}}^{i-\frac n2})\\)。
* 根据 \\(A,B\\) 的对应位相乘得到 \\(C\\) 的点值表示，时间复杂度 \\(O(nlog_2 n)\\)。
将 \(C\) 的点值表示作为多项式 \(D\) 的系数，求多项式 \(D\) 的离散傅里叶变换 \((e_0,e1,\dots,e{n-1})\)，最后用关系式 \(c_i=\frac{e_i}{n}\) 进行系数计算即可。
关于复数运算的常数优化：单位根在每次准备开始迭代时就提前算好，不要每次需要调用时才算，不然会增大常数。
这部分顺便贴个代码吧。
for (int len = 2; len <= n; len <<= 1) { cp wn(cos(P * 2 / len), inv * sin(P * 2 / len)); for (int st = 0; st + len <= n; st += len) { cp ax(1, 0); for (int j = st; j < st + (len >> 1); j++, ax = ax * wn) { cp a1 = a[j], a2 = a[j + (len >> 1)] * ax; a[j] = a1 + a2, a[j + (len >> 1)] = a1 - a2; } } }
三次变两次：直接将 \(B\) 的系数放虚部，跑离散傅里叶变换，平方，跑离散傅里叶变换，除 \(2n\)，即可得到系数。
\(\text{Reference meterial}\)
-
小学生都能看懂的FFT！！！——\(\text{Written by }\)胡小兔
教练PPT（优质PPT你值得拥有！）
\(\text{O(n log}_2\text{n)}\) 的多项式乘法：\(\text{FFT}\)
\(\text{Part 0.}\) 点值表示&系数表示
\(\text{0.1}\) 点值表示
对于次数为 \(n\) 的多项式，我们很容易发现，只要我们代入 \(n+1\) 个不同的 \(x\) 值求得 \(n+1\) 个函数值，这 \(n+1\) 个函数值完全可以代替原来多项式的系数表示，让我们能够通过拉格朗日插值法 \(O(n)\) 求出任意一个函数值。
也就是说想要通过函数值表示一个次数为 \(n\) 的多项式我们只需要 \(n+1\) 个 \(x\) 的取值不同的函数值就可以表示出这个多项式。
此时我们用 \(n+1\) 个点值表示一个次数为 \(n\) 的多项式，故称其为点值表示。
PS：拉插不是 \(FFT\) 的重点，所以这里为了日后vicky复习 \(FFT\)
时不犯迷糊着想，就不对拉插进行过多阐述。想学的话自己搜blog吧，这玩意还是蛮好学的。
\(\text{0.2}\) 系数表示
我写由题意得应该不会被骂吧。- 具体例子就是： \(A(x)=a_0+a_1x+x_2x^2+\cdots+a_nx^n\)，这就是系数表示。
\(\text{Part 1.FFT}\) 思路概述
现在我们有三个多项式:
\(A(x)=a_0+a_1x+a2x^2+\cdots+a{n-1}x^{n-1},B(x)=b0+\cdots+b{n-1}x^{n-1}\)。
\(C(x)=A(x)*B(x)=\sum\limits_{i=0}^{2n-2}c_ix^i\) 。
对于多项式 \(C(x)=A(x)*B(x)\) 而言，我们如果用点值表示这个多项式，因为它的次数为 \(2n-2\) ，所以我们只需要求出 \(2n-1\) 个 \(C(x)\) 的函数值，就可以在 \(O(n)\) 的时间复杂度内求出多项式 \(A(x)\) 和 \(B(x)\) 的乘积。
那么此时求两个多项式相乘的思路其实就很明确了：
- 将多项式 \(A(x)\) 和 \(B(x)\) 由系数表示转化为点值表示，然后然后将它们俩相乘。【每求解一个 \(C(x)\) 所需要的时间均为 \(O(n)\) ，\(C(x)\) 若要使用点值表示则需要求 \(2n+1\) 个 \(C(x)\) 的点值，所以这里如果使用朴素算法转化点值表示的话时间复杂度就是 \(O(n^2)\) 】
- 将 \(C(x)\) 的点值表示转化为系数表示，然后我们就可以非常自然地求得任意 \(C(x)\) 了。（这里用拉插的话时间复杂度又是 \(O(n^2)\) ）
- 那么， \(FFT\) 要做的其实就是用一些神奇技巧来将两个转化过程简化。
\(\text{Part 2.}\) 关于复数
这里坚持住的话后面的 \(FFT\) 理解起来其实就很简单了。
首先我们了解一些概念。
对于一个复数 \(a+bi\)，我们有两种理解方式：（因为这些都是vicky边口胡边自学的，所以有错记得跟我说一声
- 一个起点为原点的向量。
- 平面直角坐标系上的点 \((a,b)\)。
两个复数相乘的规则如下：模长相乘，幅角相加。
用代数式进行表示的话就是：
\[(a+bi)(c+di)=ac+bdi^2+adi+cbi=(ac-bd)+(ad+cb)i \]
模长就是点 \((a,b)\) 到原点的距离，这里也可以理解为向量的模长。幅角为该复数所对应的向量与 \(x\) 轴的逆时针夹角。
如下图：
图中 \(P_6\) 的幅角为 \(270\) 度，\((P_3)^2=P_6\)。
幅角相加的概念要理解到位。
介绍了基本虚数的概念之后，我们此时再引入一个新的概念：\(\omega _n^i\)。
我们以原点为圆心作一个圆，取该圆与 \(x\) 轴正半轴的交点所表示的那个复数作为第一个点，我们将其命名为 \(\omega_n^0\) 。
我们将这个圆平均分为 \(n\) 份，每一份都代表着一个在圆上的复数，即 \(\omega_n^i\) ，我们将复数 \(\omega_n^0\) 看作一个向量，对它进行若干次旋转，每次旋转都转 \(\frac{360}{n}\)度，那么当我们进行了 \(i\) 次旋转之后我们就会得到新的复数 \(\omega_n^{i\%n}\) ，其幅角为 \(\frac{360}{n}\cdot i\) 度。
如上图，\(P_0\) 为 \(\omega_8^0\) ，\(P_5\) 为 \(\omega_8^5\) ，即 \(\omega_8^0\) 经过 \(5\) 次旋转之后得到的幅角为 \(\frac{360}{8}*5=225\) 度的复数。
总结一下，\(\omega_n^i\) 的意义即为：
- 将一个以原点为圆心的圆等分 为 \(n\) 份，每一份对应圆上的一个代表复数 的点。
- 定义该圆与x轴正半轴的交点为 \(\omega_n^0\)。
- 该复数的幅角度数为 \(\frac{360}{n}\cdot i\) 度，也可以理解为它在 \(\omega_n^0\) 幅角的基础上加了 \(i\) 份度数为 \(\frac{360}{n}\) 度的幅角（即 \(\omega_n^0\) 等角度地旋转 \(i\) 次得到 \(\omega_n^i\) ）。
那么，这个 \(\omega_n^k\) 该怎么计算呢？
套公式即可：
\[\omega_n^k=\cos(2\pi\cdot \frac{k}{n})+i\cdot \sin(2\pi\cdot \frac{k}{n}) \]
具体为啥是这个公式，只要稍微了解一下三角函数的弧度制即可理解该公式。
且根据复数相乘幅角相加 的规则以及 \(\omega_n^i\) 的定义，我们很容易得到不同复数间的一些关系：
\(\omega_n^{2i}=(\omegan^i)^2=\omega{\frac{n}{2}}^i\)
\(\omega_n^i=\omega_n^{i\%n}\)
\(\omega_n^{i+\frac{n}{2}}=-\omega_n^i\)
其中最后一条关系还牵及到一个概念：共轭复数 。
我们称\(\omega_n^{i+\frac{n}{2}}\) 和 \(-\omega_n^i\) 互为共轭复数。正式一点地描述，即实部相等，虚部互为相反数的两复数互为共轭复数。从几何意义上来说，两个共轭复数所对应的向量是关于 \(x\) 轴对称的。
上面的两个式子直接套 \(\omega_n^i\) 上下标的定义就可以很容易证得了吧。QwQ
离散傅里叶变换
我知道你很慌但你先别慌。
离散傅里叶变换其实只是一个概念而已，没啥难理解的地方。
上面说过了，对于一个次数为 \(n\) 的多项式 \(f(x)\) ，若我们要对其用点值表示 ，则需要取 \(n+1\) 个不同的 \(x\) 值代入得到 \(n+1\) 个函数值来表示这个函数 \(f(x)\) 。
那么这 \(n+1\) 个不同的 \(x\) 值我们为何不去取 \((\omega{n+1}^0,\cdots,\omega{n+1}^{n})\) 上的值呢？
因为正经人谁会想到这种奇怪的取法啊喂！ QAQ
离散傅里叶变换的定义其实就是对于一个次数为 \(n\) 的函数 \(f(x)\) ，我们取函数 \(f(x)\) 在 \((\omega{n+1}^0,\cdots,\omega{n+1}^{n})\) 上的取值作为点值表示。
简而论之，一个函数在 \(n+1\) 个特殊位置上取值并将其作为该函数的点值表示，这个点值表示就是离散傅里叶变换 。
再说一遍！一个函数的离散傅里叶变换是该函数的点值表示 而不是单独一个点值！（写给我自己看的，没有任何侮辱各位智商的意思在。
\(\text{Part 3.}\) 复数在 \(\text{FFT}\) 的两个变换过程中的运用
搞定了复数！我们其实离搞懂 \(\text{FFT}\) 就差推几个式子的功夫啦！QwQ
下面的式子推导对我这样的数学弱智来说都还是很友好的，所以说其他人完全不用怕搞不懂 \(\text{FFT}\) 的式子推导就是了。
实在推不出来你把结论记下来感觉其实也没啥问题。
\(\text{Part 3.0}\) 回顾
咱们先回顾一下 \(\text{FFT}\) 的两个转化过程：
- 系数表示 - > 点值表示。（具体来说是利用多项式 \(A(x)\) 和 \(B(x)\) 的系数表示求得 \(C(x)\) 的点值表示。）
- 点值表示 - > 系数表示。（将 \(C(x)\) 的点值表示转化为系数表示。）
\(\text{Part 3.1}\) 系数表示 -> 点值表示
我们现有多项式 \(C(x)=A(x)*B(x)=c0+\cdots+c{2n-2}\cdot x^{2n-2}\)。
现在我们再弄两个多项式出来:
\[C_1(x)=c_0+c2\cdot x+\cdots+c{2n-2}\cdot x^{n-1} \]
\[C_2(x)=c_1+c3\cdot x+\cdots+c{2n-1}\cdot x^{n-1} \]
此时我们就可以将 \(C(x)\) 表示如下：
\[C(x)=C_1(x^2)+x\cdot C_2(x^2) \]
那么若我们将 \(C(x)\) 的系数表示转化为 \(C(x)\) 对应的离散傅里叶变换。
即系数表示 - > 点值表示。
为了下面的式子推导看上去更简洁，我们将 \(C(x)\) 的次数设为 \(n-1\) ，即函数 \(C(x)\) 的离散傅里叶变换在 \((\omega_n^0,\cdots,\omega_n^{n-1})\) 处取值。（不然满屏的 \(\omega_{n+1}^i\)谁受得了啊喂！QAQ
那么对于任意 \(C(x)\) 在 \(\omega_n^i\) 处的点值求值，我们可以将点值求值过程分类讨论后进行转化：
温馨提示： \(\omega_n^0=\omega_n^n=1,\omega_n^\frac{n}{2}=-1\)
- \(0\leq i<\frac{n}{2}\)：
\[C(\omega_n^i)=C_1((\omega_n^i)^2)+\omega_n^i\cdot C_2((\omega_n^i)^2) \]
\[=C_1(\omega_n^{2i})+\omega_n^i\cdot C_2(\omega_n^{2i}) \]
\[=C1(\omega\frac{n}{2}^{i})+\omega_n^i\cdot C2(\omega\frac{n}{2}^{i}) \]
- \(\frac{n}{2}\leq i\leq n\)，设 \(i=k+\frac{n}{2}(0\leq k< \frac{n}{2})\) ：
\[C(\omega_n^i)=C(\omega_n^{k+\frac{n}{2}}) \]
\[=C_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}\cdot C_2(\omega_n^{2k+n}) \]
\[=C1(\omega{\frac{n}{2}}^{k+{\frac{n}{2}}})+\omega_n^{k+\frac{n}{2}}\cdot C2(\omega{\frac{n}{2}}^{k+{\frac{n}{2}}}) \]
\[=C1(\omega{\frac{n}{2}}^{k}\cdot \omega_n^n)+\omega_n^{k+\frac{n}{2}}\cdot C2(\omega{\frac{n}{2}}^{k}\cdot \omega_n^n) \]
\[=C1(\omega{\frac{n}{2}}^{k})+(\omega_n^{k}\cdot \omega_n^{\frac{n}{2}})\cdot C2(\omega{\frac{n}{2}}^{k}) \]
\[=C1(\omega{\frac{n}{2}}^{k})+(\omega_n^{k}\cdot -1)\cdot C2(\omega{\frac{n}{2}}^{k}) \]
\[=C1(\omega{\frac{n}{2}}^{k})- \omega_n^{k}\cdot C2(\omega{\frac{n}{2}}^{k}) \]
我们对比一下两种情况下 \(C(\omega_n^i)\) 最终化出来的式子：
- \(C(\omega_n^i)=C1(\omega\frac{n}{2}^{i})+\omega_n^i\cdot C2(\omega\frac{n}{2}^{i})\)
- \(C(\omega_n^i)=C1(\omega{\frac{n}{2}}^{k})- \omega_n^{k}\cdot C2(\omega{\frac{n}{2}}^{k})\)
发现没有，我们通过式子推导，使得我们只需要求函数 \(C_1,C2\) 在 \((\omega{\frac{n}{2}}^0,\cdots,\omega_{\frac{n}{2}}^{\frac{n}{2}-1})\)处的取值就可以求得函数 \(C\) 的点值了。
现在问题就变得很简单了，我们只需递归实现即可，碰到递归边界 \(n=1\) 时我们直接套 \(C(\omega_n^i)=A(\omega_n^i)*B(\omega_n^i)\) 求值返回即可。
现在我们看回时间复杂度。
对于每一个 \(C(\omega_n^i)(0\leq i)\) ，如果用递归实现，递归边界为显然为 \(C(\omega_n^i)\) 中的 \(n\) 的值为 \(1\) 。
每次 \(C(\omega_n^i)\) 中的 \(n\) 都会缩小一半，不难得到我们要进行 \(log_2n\) 次递归，即每次求 \(C(\omega_n^i)\) 的点值的时间复杂度为 \(\text{O(log}_2\text{n)}\)，我们一共要求的点值总数与 \(n\) 成正比，所以求 \(C(x)\) 的点值表示所需的时间复杂度为 \(\text{O(n log n)}\) 。
但是，我们不能忽视的是，\(\text{O(nlog}_2\text{n)}\) 仅仅只是递归实现 \(\text{FFT}\) 的渐进时间复杂度而已，也就是说，递归实现的 \(\text{FFT}\) 自带的常数在某些凉心畜题人的数据下还是很容易被卡成 \(TLE\) 的。
事实上我们有更快的非递归 \(\text{FFT}\) 实现方式，且vicky个人认为 非递归版本并不比递归版本难写多少。
都学到这里了不把 \(\text{FFT}\) 的优化学完不就前功尽弃了嘛！QAQ
但是为了日后vicky看这篇帖子不会被冲昏头脑，所以我们先把 \(\text{FFT}\) 的基本思想讲完再进一步对 \(\text{FFT}\) 的代码实现进行讲解。QwQ
\(\text{Part 3.2}\) 点值表示 -> 系数表示
这部分主要也是推式子，而且相对上一部分来说这部分的式子推导还是友好很多的。
但是上一部分的其实也不难，直接套 \(\omega\) 的定义就能推出来了对吧。
说实话其实这部分直接记结论也完全没有问题的。
在上一部分中我们求的是 \(C(x)\) 在 \((\omega_n^0,\cdots,\omega_n^{n-1})\) 处取值的点值表示，也就是它的离散傅里叶变换，那么我们是否可以这些在点值上进行一些瞎搞去求得 \(C(x)\) 的系数表示呢？
答案是可以的。
为了下面的式子推导可读性更强，我们先将 \(C(x)\) 的离散傅里叶变换 \((C(\omega_n^0),\cdots,C(\omega_n^{n-1}))\) 分别设为数组 \((d_0,d1,\cdots,d{n-1})\) 。
我们将 \(C(x)\) 的离散傅里叶变换 \((d0,\cdots,d{n-1})\) 作为另一个次数为 \(n-1\) 的多项式 \(D(x)\) 的系数表示，即：
\[D(x)=C(\omega_n^0)+C(\omega_n^1)\cdot x+\cdots+C(\omega_n^{n-1})\cdot x^{n-1} \]
\[=d_0+d1x+\cdots+d{n-1}x^{n-1} \]
根据 \(d_i\) 的定义，显然有 \(d_i=C(\omegan^i)=\sum\limits{j=0}^{n-1}c_j\cdot (\omega_n^i)^j\)，其中 \(c_j\) 为多项式 \(C\) 的 \(j\) 次项系数。
我们再将 \(D(x)\) 的离散傅里叶变换表示成 \((e0,\cdots,e{n-1})\)，但是这里需要注意的是，我们不再是在 \((\omega_n^0,\cdots,\omega_n^{n-1})\) 处取值，而是在 \((\omega_n^{0},\omega_n^{-1},\cdots,\omega_n^{1-n})\) 处进行取值，即 \(e_k=D(\omegan^{-k})=\sum\limits{i=0}^{n-1}d_i\cdot (\omega_n^{-k})^i\) 。
那么对于每一个 \(e_i\) ，我们可以尝试一下对它进行式子推导，看看能否找出一种通过 \(e_i\) 快速求出 \(c_i\) 的办法：
\[e_k=D(\omegan^{-k})=\sum\limits{i=0}^{n-1}d_i\cdot (\omegan^{-k})^i=\sum\limits{i=0}^{n-1}D(\omega_n^{i})\cdot (\omega_n^{-k})^i \]
\[=\sum\limits{i=0}^{n-1}\left(\sum\limits{j=0}^{n-1}c_j\cdot (\omega_n^i)^j\right)\cdot (\omega_n^{-k})^i \]
\[=\sum\limits{i=0}^{n-1}\sum\limits{j=0}^{n-1}c_j\cdot (\omega_n^i)^j\cdot (\omegan^{-k})^i=\sum\limits{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_j\cdot \omega_n^{i\cdot j}\cdot \omega_n^{-k\cdot i} \]
\[=\sum\limits{i=0}^{n-1}\sum\limits{j=0}^{n-1}c_j\cdot \omegan^{i\cdot (j-k)}=\sum\limits{i=0}^{n-1}\sum\limits_{j=0}^{n-1}c_j\cdot (\omega_n^{j-k})^i \]
\[=\sum\limits_{j=0}^{n-1}cj\cdot \sum\limits{i=0}^{n-1}(\omega_n^{j-k})^i \]
\(k\) 显然可以理解成一个定值，我们枚举到每一个 \(j\) 的时候，\(j\) 其实也算作是一个定值，此时不难发现 \(\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\) 其实就是个等比数列求和。
因为 \(j\) 的值是在 \([0,n-1]\) 范围内进行枚举的，因此我们可以对 \(j-k\) 进行分类讨论，然后再求 \(\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\) 的值：
- \(j-k=0\)，即 \(j=k\) ：
此时显然有：
\[\sum\limits_{i=0}^{n-1}(\omegan^{j-k})^i=\sum\limits{i=0}^{n-1}(\omegan^0)^i=\sum\limits{i=0}^{n-1}1^i=n \]
- \(j-k\neq 0\)：
根据等比数列求和公式 \(sum=\dfrac{a_1\cdot (1-q^n)}{1-q}\) ，我们易得：
PS ：公式中的 \(a_1\) 为数列首项，\(q\) 为公比，\(n\) 为项数。
为了下面的式子推导看起来更简洁一些，我们设 \(W=j-k\) 。
注意，下面式子中的 \(j\) 与 \(k\) 我们均看作定值。
\[\sum\limits_{i=0}^{n-1}(\omegan^{j-k})^i=\sum\limits{i=0}^{n-1}(\omega_n^W)^i \]
\[=\dfrac{1\cdot (1-(\omega_n^W)^n)}{1-\omega_n^W}=\dfrac{1-(\omega_n^n)^W}{1-\omega_n^W} \]
\[=\dfrac{1-1^W}{1-\omega_n^W}=\dfrac{1-1}{1-\omega_n^W}=0 \]
也就是说，当且仅当 \(\sum\limits{j=0}^{n-1}\) 枚举到 \(k\) 时，\(\sum\limits{i=0}^{n-1}(\omegan^{j-k})^i\) 才会有值，且此时\(\sum\limits{i=0}^{n-1}(\omega_n^{j-k})^i\) 的值为 \(n\)。
那么我们就可以进一步转化 \(e_k\) ：
\[ek=\sum\limits{j=0}^{n-1}cj\cdot \sum\limits{i=0}^{n-1}(\omega_n^{j-k})^i \]
我们按上面分类讨论的两种情况进行求值：
\[ek=\left(\sum\limits{j=0,j\neq k}^{n-1}cj \cdot \sum\limits{i=0}^{n-1}(\omega_n^{j-k})^i\right)+ck\cdot \sum\limits{i=0}^{n-1}(\omega_n^0)^i \]
\[=\left(\sum\limits_{j=0,j\neq k}^{n-1}c_j \cdot 0\right)+c_k\cdot n \]
\[=c_k\cdot n \]
这个时候要怎么求 \(c_k\) 就不用我说了吧。 QwQ
\[c_k=\dfrac{e_k}{n} \]
此时我们运用上一部分中的系数表示- >点值表示的方法即可在 \(\text{O(log}_2\text{n)}\) 的时间复杂度下求每个 \(e_k\) 也就是 \(D(\omega_n^{-k})\) 的值，每求出一个 \(e_k\) ，我们就可以运用 \(c_k=\dfrac{e_k}{n}\) 这一关系 \(\text{O(1)}\) 求出 \(C(x)\) 的第 \(k\) 项系数。
而 \(C(x)\) 的项数与 \(n\) 成正比，所以我们运用这一方法实现点值表示- >系数表示的渐进时间复杂度即为 \(\text{O(nlog}_2\text{n)}\)。
\(\text{Part 3.3}\) 回顾 \(\text{FFT}\) 的全过程
\(\text{FFT}\) 是一种可以在 \(\text{O(n log}_2\text{n)}\) 的渐进时间复杂度内求次数都为 \(n\) 的两个多项式 \(A(x)\) 和 \(B(x)\) 的卷积 \(C(x)\) 的值的算法。
\(\text{Part 3.1}\)
- 设 \(C(x)=A(x)*B(x)\)，分治求\(A(x),B(x)\) 的离散傅里叶变换，然后求对应的 \(C(x)\) 的离散傅里叶变换，渐进时间复杂度为 \(\text{O(nlog}_2\text{n)}\)。
\(\text{Part 3.2}\)
- 以 \(C(x)\) 的离散傅里叶变换作为另一个与 \(C(x)\) 同次的多项式 \(D(x)\) 的系数，每次用 \(\text{Part 3.1}\) 中的办法 \(\text{O(log}2\text{n)}\) 求每个 \(D(x)\) 在 \((\omega{2n-1}^{0},\omega{2n-1}^{-1},\cdots,\omega{2n-1}^{-(2n-2)})\) 处取值的得到点值 \(e_i\)，根据 \(c_i=\dfrac{e_i}{n}\) 的关系求 \(C(x)\) 的每一项系数，求 \(C(x)\) 的 \(2n-2\) 项系数的渐进时间复杂度为 \(\text{O(n log}_2\text{n)}\) 。
- 求出 \(C(x)\) 的每一项系数之后，因为 \(C(x)=A(x)*B(x)\) ，所以对于多项式 \(A(x)\) 和 \(B(x)\) 的乘积，我们直接将 \(x\) 代入系数表示下 \(C(x)\) 在 \(\text{O(n)}\) 渐进时间复杂度下进行求解即可。
所以 \(\text{FFT}\) 的渐进时间复杂度为 \(\text{O(nlog}_2\text{n)}\) 。
放一下板题（P3803 【模板】多项式乘法（FFT））递归写法的代码：
#include<bits/stdc++.h>
#define N 4000005
#define P acos(-1.0)
using namespace std;
int aa,bb;
struct cp{
double x,y;
cp (double xx=0,double yy=0){ x=xx,y=yy; }
}a[N],b[N];
cp operator +(cp a,cp b){ return cp(a.x+b.x,a.y+b.y);}
cp operator -(cp a,cp b){ return cp(a.x-b.x,a.y-b.y);}
cp operator *(cp a,cp b){ return cp(a.x*b.x-a.y*b.y,a.y*b.x+a.x*b.y);}
void fft(cp*,int,int);
int main(){
scanf("%d%d",&aa,&bb);
for(int i=0;i<=aa;i++)
scanf("%lf",&a[i].x);
for(int i=0;i<=bb;i++)
scanf("%lf",&b[i].x);
int n=1;
while(n<=aa+bb) 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);
for(int i=0;i<=aa+bb;i++)
printf("%d ",(int)(a[i].x/n+0.5));
return 0;
}
void fft(cp *a,int n,int inv){
if(n==1) return ;
cp a1[(n>>1)+5],a2[(n>>1)+5];
for(int i=0;i<n;i+=2)//这里是i+=2！！！
a1[i>>1]=a[i],a2[i>>1]=a[i+1];
fft(a1,n>>1,inv),fft(a2,n>>1,inv);
for(int i=0;i<(n>>1);i++){
cp x(cos(2.0*P*i/n),inv*sin(2.0*P*i/n)),ax=x*a2[i];
a[i]=a1[i]+ax,a[i+(n>>1)]=a1[i]-ax;
}
}
\(\text{Part 3.4}\) \(\text{FFT}\) 优化
优化一：递归->迭代
迭代即用循环模拟递归的过程，对 FFT 的常数有不小的优化作用。
我们观察一下在做FFT的时候多项式 \(C(x)\) 的每一项系数的位置变化。
0 1 2 3 4 5 6 7
0 2 4 6|1 3 5 7
0 4|2 6|1 5|3 7
0|4|2|6|1|5|3|7
这里有一个很奇妙的性质：对于每一个系数，它的最终下标为它原来下标的二进制翻转。
如原来处在位置 \(6(=(110)_2)\) 的系数最终到了位置 \(3(=(011)_2)\) 。
那么此时我们其实就可以把要转换的多项式 \(C(x)\) 的系数一开始就调换到它的最终位置，然后再通过循环模拟递归的过程。
代码也很简单：
int sum_len = len_a+len_b, len = 0;
while (n <= sum_len)
++len, n <<= 1;
for (int i = 1; i < n; i++)
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
for (int i = 0; i < n; i++)
if (i < r[i]) swap(a[i], a[r[i]]);
for (int len = 2; len <= n; len <<= 1) {
cp wn(cos(P * 2 / len), inv * sin(P * 2 / len));
for (int st = 0; st + len <= n; st += len) {
cp ax(1, 0);
for (int j = st; j < st + (len >> 1); j++, ax = ax * wn) {
cp a1 = a[j], a2 = a[j + (len >> 1)] * ax;
a[j] = a1 + a2, a[j + (len >> 1)] = a1 - a2;
}
}
}
优化二：减少复数乘法运算次数
而且实现FFT时还有一点十分重要，就是转化过程中的复数乘法运算。
有图有真相：D
让我们来一起看看两份代码有什么区别吧：D
\(5.42s\)：
for(int len=2;len<=n;len<<=1){
for(int st=0;st+len-1<n;st+=len){
for(int j=st;j<=st+(len>>1)-1;j++){
cp x(cos(2.0*P*(j-st)/len),inv*sin(2.0*P*(j-st)/len)),ax=x*a[j+(len>>1)],a1=a[j];
a[j]=a1+ax,a[j+(len>>1)]=a1-ax;
}
}
}
\(2.83s\)：
for(int len=2;len<=n;len<<=1){
complex wn(cos(2*P/len),inv*sin(2*P/len));
for(int st=0;st+len<=n;st+=len){
complex a1,a2,ax(1,0);
for(int j=st;j<st+(len>>1);j++,ax=ax*wn)
a1=a[j],a2=a[j+(len>>1)]*ax,
a[j]=a1+a2,a[j+(len>>1)]=a1-a2;
}
}
两份代码的其余部分均相同。
两份代码仅仅差了几次的复数运算，最终的运行结果却相差甚远。
所以说，复数乘法这种东西，咱们还是能少算几次就少算几次 。
在第二份代码中我们就少算了很多次 (cos(2*P/len),inv*sin(2*P/len)) 的值，故常数优化了不少。
因为这里是我感性理解的，并没有找理性证明的博客进行学习，所以如果我说错了的话麻烦在评论区指出一下我的错误 www。awa
优化三：三次FFT变为两次FFT
在我们求 \(A(x)\) 和 \(B(x)\) 的卷积 \(C(x)\) 时，我们会对它们分别做 \(\text{FFT}\) 转化为点值表示，然后对应每项乘起来得到 \(C(x)\) 的点值表示，最后再做一次 \(\text{FFT}\) 将 \(C(x)\) 由点值表示转化为系数表示。
这样的话实际上我们一共进行了三次 \(\text{FFT}\) ，于是我们考虑能否减少做 \(\text{FFT}\) 的次数。
当然可以啊，不然我怎么可能写这个板块啊。
我们弄两个系数为复数的多项式，按正常 \(\text{FFT}\) 的做法，我们应该分别把多项式 \(A(x)\) 和 \(B(x)\) 的系数放到两个多项式系数的实部上面去。
但是在三次变两次 \(\text{FFT}\) 优化中，我们只需要弄一个系数为复数的多项式，并将多项式 \(A(x)\) 和 \(B(x)\) 的系数分别放到该多项式系数的实部和虚部中，然后对这个多项式跑一边 \(\text{FFT}\) ，然后对它做平方，再做一遍 \(\text{FFT}\) 逆变换回来,最后取它每一位系数虚部上的数字再 \(/2n\) 即可（该多项式次数为 \(n\)）。
证明也很简单，了解完全平方公式以及复数概念即可证明。
\[(a+bi)^2=a^2+(b^2i^2)+2abi=a^2-b^2+2abi \]
此时我们的常数就可以优化到约为原来的 \(\frac{2}{3}\) 了。
这是三个优化都加上了的代码跑出来的：D
最后放个加了优化的代码QwQ
#include<bits/stdc++.h>
#define P acos(-1.0)
#define N 1000005
using namespace std;
int a1,b1,n=1,r[N<<2],ans[N<<2],aa1;
struct cp{
double x,y;
cp(double xx=0,double yy=0){ x=xx,y=yy;}
}a[N<<2],b[N<<2];//这里的复数我用的是结构体表示，实际上复数运算也可以用STL
cp operator +(cp a,cp b){ return cp(a.x+b.x,a.y+b.y);}
cp operator -(cp a,cp b){ return cp(a.x-b.x,a.y-b.y);}
cp operator *(cp a,cp b){ return cp(a.x*b.x-b.y*a.y,a.y*b.x+b.y*a.x);}
void pre(),fft(cp*,int);
int main(){
scanf("%d%d",&a1,&b1);
for(int i=0;i<=a1;i++) scanf("%lf",&a[i].x);
for(int i=0;i<=b1;i++) scanf("%lf",&a[i].y);
pre(),fft(a,1);
for(int i=0;i<n;i++) a[i]=a[i]*a[i];
fft(a,-1);
for(int i=0;i<=a1+b1;i++) printf("%d ",(int)(a[i].y/n/2+0.5));//取虚部求最终的系数表示
return 0;
}
void pre(){
aa1=max(a1,b1)<<1;
while(n<=aa1) n<<=1;
int len=1;
while((1<<len)<=aa1) ++len;
for(int i=1;i<n;i++) r[i]=((r[i>>1]>>1|(i&1)<<(len-1)));//预处理每一个系数的最终位置（二进制翻转）
}
void fft(cp *a,int inv){
for(int i=0;i<n;i++)
if(i<r[i]) swap(a[r[i]],a[i]);
//用循环模拟递归的过程
for(int len=2;len<=n;len<<=1){
cp wn(cos(P*2/len),inv*sin(P*2/len));//提前算好单位复数根，减少重复运算
for(int st=0;st+len-1<n;st+=len){
cp ax(1,0);
for(int j=st;j<st+(len>>1);j++,ax=ax*wn){
cp a1=a[j],a2=ax*a[j+(len>>1)];
a[j]=a1+a2,a[j+(len>>1)]=a1-a2;
}
}
}
}
/xin-xi-ji-zhu-fu-wu/fftke-ai-qwq-7352.html