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+x_3\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,e_1,\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+a_2x^2+\cdots+a_{n-1}x^{n-1},B(x)=b_0+\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}=(\omega_n^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)=c_0+\cdots+c_{2n-2}\cdot x^{2n-2}\)。

现在我们再弄两个多项式出来:

\[C_1(x)=c_0+c_2\cdot x+\cdots+c_{2n-2}\cdot x^{n-1} \]

\[C_2(x)=c_1+c_3\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}) \]

\[=C_1(\omega_\frac{n}{2}^{i})+\omega_n^i\cdot C_2(\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}) \]

\[=C_1(\omega_{\frac{n}{2}}^{k+{\frac{n}{2}}})+\omega_n^{k+\frac{n}{2}}\cdot C_2(\omega_{\frac{n}{2}}^{k+{\frac{n}{2}}}) \]

\[=C_1(\omega_{\frac{n}{2}}^{k}\cdot \omega_n^n)+\omega_n^{k+\frac{n}{2}}\cdot C_2(\omega_{\frac{n}{2}}^{k}\cdot \omega_n^n) \]

\[=C_1(\omega_{\frac{n}{2}}^{k})+(\omega_n^{k}\cdot \omega_n^{\frac{n}{2}})\cdot C_2(\omega_{\frac{n}{2}}^{k}) \]

\[=C_1(\omega_{\frac{n}{2}}^{k})+(\omega_n^{k}\cdot -1)\cdot C_2(\omega_{\frac{n}{2}}^{k}) \]

\[=C_1(\omega_{\frac{n}{2}}^{k})- \omega_n^{k}\cdot C_2(\omega_{\frac{n}{2}}^{k}) \]

我们对比一下两种情况下 \(C(\omega_n^i)\) 最终化出来的式子：

\(C(\omega_n^i)=C_1(\omega_\frac{n}{2}^{i})+\omega_n^i\cdot C_2(\omega_\frac{n}{2}^{i})\)
\(C(\omega_n^i)=C_1(\omega_{\frac{n}{2}}^{k})- \omega_n^{k}\cdot C_2(\omega_{\frac{n}{2}}^{k})\)

发现没有，我们通过式子推导，使得我们只需要求函数 \(C_1,C_2\) 在 \((\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<n)\) ，如果用递归实现，递归边界为显然为 \(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,d_1,\cdots,d_{n-1})\) 。

我们将 \(C(x)\) 的离散傅里叶变换 \((d_0,\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+d_1x+\cdots+d_{n-1}x^{n-1} \]

根据 \(d_i\) 的定义，显然有 \(d_i=C(\omega_n^i)=\sum\limits_{j=0}^{n-1}c_j\cdot (\omega_n^i)^j\)，其中 \(c_j\) 为多项式 \(C\) 的 \(j\) 次项系数。

我们再将 \(D(x)\) 的离散傅里叶变换表示成 \((e_0,\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(\omega_n^{-k})=\sum\limits_{i=0}^{n-1}d_i\cdot (\omega_n^{-k})^i\) 。

那么对于每一个 \(e_i\) ，我们可以尝试一下对它进行式子推导，看看能否找出一种通过 \(e_i\) 快速求出 \(c_i\) 的办法：

\[e_k=D(\omega_n^{-k})=\sum\limits_{i=0}^{n-1}d_i\cdot (\omega_n^{-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 (\omega_n^{-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 \omega_n^{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}c_j\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}(\omega_n^{j-k})^i=\sum\limits_{i=0}^{n-1}(\omega_n^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}(\omega_n^{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}(\omega_n^{j-k})^i\) 才会有值，且此时\(\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\) 的值为 \(n\)。

那么我们就可以进一步转化 \(e_k\) ：

\[e_k=\sum\limits_{j=0}^{n-1}c_j\cdot \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i \]

我们按上面分类讨论的两种情况进行求值：

\[e_k=\left(\sum\limits_{j=0,j\neq k}^{n-1}c_j \cdot \sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i\right)+c_k\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;
			}
		}
	}
}

posted @ 2023-04-26 20:13 vickyの橙子林阅读(76) 评论(0) 编辑收藏举报

刷新页面返回顶部

vickyqwq

别为失败做打算，那比庸常的打算还愚蠢。