前言（不想听的可以跳到下面）

OK。蒟蒻又来口胡了。
自从ZJOI2019上Day的数论课上的多项式听到懵逼了，所以我就下定决心要学好多项式。感觉自己以前学的多项式都是假的。
但是一直在咕咕，现在是中午，一个早上的努力就完成了FFT的学习，其实并没有想象中的那么难。
文笔较渣，想到什么就写什么，可能逻辑性比较差，来回看个几遍差不多就懂了。

介绍

先简单介绍一下FFT（Fast Fourier Transformation） ，中文全名叫做快速傅里叶变换。
应用在加速多项式的乘法，或者是高精度加速（就是以$10$10为底的多项式）

前置芝士

一、欧拉公式

这个东西啊，被称为数学界最优美的公式。我也不知道为什么（雾）
对于任意的$x$x，都有以下的性质
$e^{ix}=cos(x)+i \times sin(x)$eix=cos(x)+i×sin(x)

我想来证明一下：显然
跟上自己的渣渣证明：【传送门】
这个东西太变态了，好像还有什么别的什么微积分证明，不管了，不会。。

二、复数

复数的定义

形如：$a+bi$a+bi的数，叫做复数。
那么其中的$i$i是$\sqrt{-1}$−1​，$i$i就是在复数虚部的基本单位。
引入了虚数的概念，我们就从数轴到了复数的平面。
模长：这个定义和向量中的定义差不多，也是表示线段的长度。
幅角：以$x$x轴的正半轴为$0°$0°，到向量的夹角为幅角。

复数的计算法则

• 加法：虚实部分别相加
• 减法：虚实部分别相减
• 乘法：其实和代数的乘法差不多，我们可以把这个东西看成两个$2$2的代数式的乘法。

稍微推导一下虚数的乘法

(b+ai)*(d+ci)//展开
=ai*ci+b*ci+ai*d+b*d//因为i*i=-1那么就可以继续化简合并同类项
=(b*d-a*c)+(b*c+a*d)i

以上这个东西是不是非常个简单，是不是觉得自己是一个天才。（滑稽脸）
虽然C++的STL中有complex的库，但是还是手写的比较放心，没准被卡掉了就哭死了qwq
给一下自己手写的complex的结构体，注意需要大写首字母或者是改变几个字母，不然会变量名冲突导致CE。

struct Complex {
    db x, y;
    Complex (db xx = 0.0, db yy = 0.0) { x = xx, y = yy; }
    Complex operator + (Complex B) { return Complex(x + B.x, y + B.y); }
    Complex operator - (Complex B) { return Complex(x - B.x, y - B.y); }
    Complex operator * (Complex B) { return Complex(x * B.x - y * B.y, x * B.y + y * B.x); }
}

三、向量

我们还稍微需要知道一点向量的知识，可能可以抽象成这样的模型。
不详细讲解。
平行四边形法则：

$AB+AD=AC$AB+AD=AC

这个东西适用于向量之间的加法，但是也可以适用于复数之间的加法。
怎么讲这句话。
因为在复数的平面内，复数可以看成一个一个向量。
还有一个定理是向量的乘法是模长相乘，幅角相加
上面提过了模长和幅角，大致理解一下就可以了。

单位根

简单介绍

这个东西极其的重要，如果这个东西不理解，FFT是很难学好的。
默认$n$n为$2$2的若干次幂。
在复平面内，以原点为圆心，$1$1为半径，做一个圆，叫做单位圆。
单位根，就是将一个单位圆平均分成$n$n分，形成了$n$n个向量，然后每一个向量所代表的的复数值就是我们的$\omega$ω。
那么单位根就是$\omega^1_n$ωn1​。

以上就是把一个单位圆分成了$8$8份，也就是$8$8个向量的情况。
我们从$x$x轴正半轴上的第一个点开始标号，从$0$0开始。
那么可以得到$\omega^{0\cdots n}_n$ωn0⋯n​，其中这个$n$n和$0$0是相等向量，方向长度都相等。

在代数中，如果$\omega^{n}=1$ωn=1，那么称$\omega$ω为$n$n次单位根。和我们FFT没有太大关系

我们先记$\omega^k_n$ωnk​表示的是编号为$k$k的点的复数值，那么根据模长相乘，幅角相加的定理得到$\omega^k_m=$ωmk​=
各个单位根的复数值，乍一看还是非常不好求的，那么我们就用欧拉公式来求解。

---

插入一个东西

弧度制，等于半径长的圆弧所对的圆心角叫做1弧度的角，用符号rad表示，读作弧度。用弧度作单位来度量角的制度叫做弧度制。

---

回归正文。
将$\omega$ω代入欧拉公式的每一个点的权值为$\omega^k_n=cos(k \times \frac{2\pi}{n})+i\times sin(k \times \frac{2\pi}{n})$ωnk​=cos(k×n2π​)+i×sin(k×n2π​)这里用到了弧度制的知识，高中必修应该有的吧。。。

单位根的性质

这个东西有一点神奇。
首先$\omega^0_n=\omega^n_n=0$ωn0​=ωnn​=0，上文应该有提到过吧。
下面两个性质，也是FFT能够加速成$nlogn的原因$nlogn的原因

性质1

$\omega^{2k}_{2n}=\omega^k_n$ω2n2k​=ωnk​

证明
$\omega^{2k}_{2n}=cos \ 2k\times \frac{2\pi}{2n}+i \ sin 2k \times \frac{2\pi}{2n}$ω2n2k​=cos 2k×2n2π​+i sin2k×2n2π​

$=\omega^k_n$=ωnk​

性质2

$\omega^{k+\frac{n}{2}}_n=-\omega^k_n$ωnk+2n​​=−ωnk​

证明
$\omega^{\frac n2}_{n}=cos \ \frac{n}{2}\times \frac{2\pi}{n}+i \ sin \frac n2\times \frac{2\pi}{n}$ωn2n​​=cos 2n​×n2π​+i sin2n​×n2π​

$=cos \pi +i \ sin \ \pi$=cosπ+i sin π

$= -1$=−1

好像还有人不知道为什么这个$=cos \pi +i \ sin \ \pi=-1$=cosπ+i sin π=−1，其实很简单，回到上面那个图，然后你就会发现$\pi$π的弧度制，表示的刚好是半个圆，那么就是$-1$−1。不需要证明。

四、卷积

在泛函分析中，卷积、旋积或摺积(英语：Convolution)是通过两个函数f 和g 生成第三个函数的一种数学算子，表征函数f 与g经过翻转和平移的重叠部分的面积。by baidubaike

这样讲有一点空泛，那么来举一个事例。
比如说多项式之间的乘法就是卷积。
个人理解卷积就是两个多项式或者是函数通过某一种计算法则，变成另一个多项式或者是函数。
因为有系数的区别，那么会有平移或者是翻转的情况。

---

区分一些概念

FFT是快速傅里叶变换
IFFT是快速傅里叶逆变换
DFT是离散傅里叶变换
IDFT是离散傅里叶逆变换
其实比较好理解，逆变换和变换之间的关系就是变换是把一个式子变成了点值表达的形式，逆变换就是把这个东西变回来。（不对好像剧透了什么，算了这个博客来回看个两遍应该就全都懂了。。。）

正文

多项式的两种表示方法

方法1-系数表达法

也就是我们的定义，任意的$n$n次多项式$f(x)$f(x)都可以表示成$\sum^{n}_{i=0} a_ix^i$∑i=0n​ai​xi的形式。
差不多就是长以下这个样子：
$f(x)=a$f(x)=a
很久之前我们学过高精度，同理可以用$O(n^2)$O(n2)的时间暴力求出我们的答案。
大致框架如下

for (int i = 1; i <= n; i ++) {
	for (int j = 1; j <= n; j ++) ans[i * j - 1] += a[i] * b[j];
}

那么这个东西也没有什么可以挖下去的可能了，我们就放弃这个$O(n^2)$O(n2)算法吧。

方法2-点值表达法

给你$n+1$n+1个点，你可以确定一个$n$n次多项式

给你一个非常简单的证明：小学生都知道$n$n个未知数的方程，如果有$n+1$n+1个方程，就可以求出各个根的值很多时候是太烦了所以不想算（大雾
同理我们就可以把一个多项式表示成$n+1$n+1个点的形式
就是以下的形式
$f(x)=((x_1,y_1),(x_2,y_2),(x_3,y_3)\cdots)$f(x)=((x1​,y1​),(x2​,y2​),(x3​,y3​)⋯)

---

DFT（离散傅里叶变换）

那么可以回到正题了，上文我们说到了点值表达法，那么我们来研究一下点值表达法的一些特点。
比较容易可以发现，如果我们已经得到了点值的表达式，那么我们可以在$O(1)$O(1)时间内卷积。
比如说：多项式$f(x)$f(x)和$g(x)$g(x)，都是$n$n次的多项式。
$f(x)=((x_0,f(x_0)),(x_1,f(x_1))\cdots)$f(x)=((x0​,f(x0​)),(x1​,f(x1​))⋯)

$g(x)=((x_0,g(x_0)),(x_1,g(x_1))\cdots)$g(x)=((x0​,g(x0​)),(x1​,g(x1​))⋯)

那么卷积生成的多项式设为$h(x)$h(x)。
$h(x)=((x_0,f(x_0)\times g(x_0)),(x_1,f(x_1)\times g(x_1)),\cdots)$h(x)=((x0​,f(x0​)×g(x0​)),(x1​,f(x1​)×g(x1​)),⋯)

也就是$x$x坐标不变，$y$y坐标相乘。
注意这里的时间复杂度好像只有$O(n)$O(n)，需要注意我们系数转点值和点值转系数的过程如果暴力依旧是$O(n^2)$O(n2)。

---

那么我们继续。
前置芝士中提过了，对于任意一个$n$n次多项式，如果给定$n+1$n+1个点，那么就能够确定。
但是很明显，难道这些点是随机取的吗？当然不是。还记不记得我们刚才有一个叫做单位根的东西。
那么傅里叶就提出了，取这些单位圆上的点就可以解决问题。
傅里叶是这样想的：如果代入的值的若干次方是$1$1，那么就不需要计算所有的值了。
那么我们就把这些单位根$\omega$ω代入多项式中，然后求解。
但是很可惜傅里叶不是一个信息学家，他又不关这些时间复杂度。DFT的时间复杂度依旧是$O(n^2)$O(n2)。（傅里叶：这锅我不背，时辰出来背锅）

IDFT（离散傅里叶逆变换）

至于点值变成系数的过程，也很显然易见，通过单位根转回去的时间复杂度依旧是$O(n^2)$O(n2)。
那么多项式的操作也就只能停步在$O(n^2)$O(n2)了吗？
前方高能，我已经把车门锁紧了，坐稳了。

FFT（快速傅里叶变换）

先给出一点提示，FFT的时间复杂度是$O(nlogn)$O(nlogn)，那么很明显，我们可以分治或者是二分。
我们先假设一个多项式是$A(x)=\sum^{n-1}_{i=0}a_i \times x^i$A(x)=∑i=0n−1​ai​×xi。展开来就是$=a_0+a_1x+a_2x^2+\cdots +a_{n-1}x^{n-1}$=a0​+a1​x+a2​x2+⋯+an−1​xn−1。
按照多项式下标的奇偶分一个组。
$A(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1})$A(x)=(a0​+a2​x2+⋯+an−2​xn−2)+(a1​x+a3​x3+⋯+an−1​xn−1)

右边全部提出来一个$x$x
$A(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2})$A(x)=(a0​+a2​x2+⋯+an−2​xn−2)+x(a1​+a3​x2+⋯+an−1​xn−2)

是不是稍微看出了一点什么东西。
然后我们设$A_1(x)=a_0+a_2x^2+\cdots+a_{n-2}x^{n-2}$A1​(x)=a0​+a2​x2+⋯+an−2​xn−2，和$A_2=a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}$A2​=a1​+a3​x2+⋯+an−1​xn−2，是不是有得到了两个多项式
代回去得到了
$A(x)=A_1(x^2)+xA_2(x^2)$A(x)=A1​(x2)+xA2​(x2)

代入我们点值表达式的单位根
因为分上下半圆，我们就先代入上半圆$\omega^k_n(k&lt;\frac{n}{2})$ωnk​(k<2n​)
$A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2)$A(ωnk​)=A1​((ωnk​)2)+ωnk​A2​((ωnk​)2)

$A(\omega^k_n)=A_1(\omega^{2k}_n)+\omega^k_nA2(\omega^{2k}_n)=A_1(\omega^k_{\frac{n}{2}})+\omega^k_nA_2(\omega^k_{\frac{n}{2}})$A(ωnk​)=A1​(ωn2k​)+ωnk​A2(ωn2k​)=A1​(ω2n​k​)+ωnk​A2​(ω2n​k​)
那么代入下半圆，也就是代入$\omega^{k+\frac{n}{2}}_n(k&lt;\frac{n}{2})$ωnk+2n​​(k<2n​)
$A(\omega^{k+\frac{n}{2}}_n)=A_1((\omega^{k+\frac{n}{2}}_n)^2)+\omega^{k+\frac n2 }_nA_2((\omega^{k+\frac{n}{2}}_n)^2)$A(ωnk+2n​​)=A1​((ωnk+2n​​)2)+ωnk+2n​​A2​((ωnk+2n​​)2)

$=A_1(\omega^{2k+n}_n)+\omega^{k+\frac{n}{2}}_nA2(\omega^{2k+n}_n)$=A1​(ωn2k+n​)+ωnk+2n​​A2(ωn2k+n​)

$=A_1(\omega^{2k}_n\times \omega^n_n)+\omega^{ \frac n2}_n \omega^k_nA2(\omega^{2k}_n\times\omega^n_n)$=A1​(ωn2k​×ωnn​)+ωn2n​​ωnk​A2(ωn2k​×ωnn​)

$=A_1(\omega^{2k}_n)-\omega^k_nA2(\omega^{2k}_n)$=A1​(ωn2k​)−ωnk​A2(ωn2k​)

$=A_1(\omega^k_{\frac{n}{2}})-\omega^k_nA_2(\omega^k_{\frac{n}{2}})$=A1​(ω2n​k​)−ωnk​A2​(ω2n​k​)

发现了什么这两堆东西，之后后面的符号不同。
那么也就是说，在我们算出了$A_1(\omega^k_{\frac{n}{2}})$A1​(ω2n​k​)和$A_2(\omega^k_{\frac{n}{2}})$A2​(ω2n​k​)的时候，我们就可以得到$A_1(\omega^{k+\frac{n}{2}}_{\frac{n}{2}})$A1​(ω2n​k+2n​​)和$A_2(\omega^{k+\frac{n}{2}}_{\frac{n}{2}})$A2​(ω2n​k+2n​​)的值了，全部计算完，也就可以直接算$\omega^k_n$ωnk​和$\omega^{k+\frac{n}{2}}_n$ωnk+2n​​的值了。
也就是说我们在询问区间$[0,\frac{n}{2})$[0,2n​)的时候，我们已经把区间$[\frac{n}{2},n)$[2n​,n)已经算出来了。
哇塞，好像我们的问题就直接折半了，在用递归求解，那么在只有一项的时候返回常数项就可以了。
是不是非常牛逼。。。
这样写就直接把时间复杂度降到$O(nlogn)$O(nlogn)了。

IFFT（快速傅里叶逆变换）

好像问题还没有结束，因为我们把系数表示法，转成了点值表示法，但是转不回去。。。
那么我们就需要把点值转回系数。
下面的推导来自【自为风月马前卒大佬的博客】
$(y_0,y_1,y_2,…,y_{n−1})为(a_0,a_1,a_2,…,a_{n−1})$(y0​,y1​,y2​,…,yn−1​)为(a0​,a1​,a2​,…,an−1​)的傅里叶变换（即点值表示）
$c_k=∑_{i=0}^{n−1}y_i(ω^{−k}_n)^i$ck​=i=0∑n−1​yi​(ωn−k​)i

$=∑^{i=0} _{n−1}(∑^{j=0}_{n−1}a_j(ω^i_n)^j)(ω^{−k}_n)^i$=n−1∑i=0​(n−1∑j=0​aj​(ωni​)j)(ωn−k​)i

$=∑^{i=0}_{n−1}(∑^{j=0}_{n−1}a_j(ω^j_n)^i)(ω^{−k}_n)^i$=n−1∑i=0​(n−1∑j=0​aj​(ωnj​)i)(ωn−k​)i

$=∑_{i=0}^{n−1}(∑_{j=0}^{n−1}a_j(ω^j_n)^i(ω^{−k}_n)^i)$=i=0∑n−1​(j=0∑n−1​aj​(ωnj​)i(ωn−k​)i)

$=∑_{i=0}^{n−1}∑_{j=0}^{n−1}a_j(ω^j_n)^i(ω^{−k}_n)^i$=i=0∑n−1​j=0∑n−1​aj​(ωnj​)i(ωn−k​)i

$=∑_{i=0}^{n−1}∑_{j=0}^{n−1}a_j(ω^{j−k}_n)^i$=i=0∑n−1​j=0∑n−1​aj​(ωnj−k​)i

$=∑_{j=0}^{n−1}a_j(∑_{i=0}^{n−1}(ω^{j−k}_n)^i)$=j=0∑n−1​aj​(i=0∑n−1​(ωnj−k​)i)

设$S(x)=∑^{n−1}_{i=0}x^i$S(x)=∑i=0n−1​xi
$S(ω^k_n)=1+(ω^k_n)+(ω^k_n)^2+\cdots+(ω^k_n)^{n−1}$S(ωnk​)=1+(ωnk​)+(ωnk​)2+⋯+(ωnk​)n−1

当$k!=0$k!=0时，等式两边同乘$ω^k_n$ωnk​得
$ω^k_nS(ω^k_n)=ω^k_n+(ω^k_n)^2+(ω^k_n)^3+\cdots(ω^k_n)^n$ωnk​S(ωnk​)=ωnk​+(ωnk​)2+(ωnk​)3+⋯(ωnk​)n

两式相减得
$ω^k_nS(ω^k_n)−S(ω^k_n)=(ω^k_n)^n−1$ωnk​S(ωnk​)−S(ωnk​)=(ωnk​)n−1

$S(ω^k_n)=\frac{(ω^k_n)^n−1}{ω^k_n−1}$S(ωnk​)=ωnk​−1(ωnk​)n−1​

$S(ω^k_n)=\frac{(ω^n_n)^k−1}{ω^k_n−1}$S(ωnk​)=ωnk​−1(ωnn​)k−1​

$S(ω^k_n)=\frac{1−1}{ω^k_n−1}$S(ωnk​)=ωnk​−11−1​

观察这个式子，不难看出它分母不为$0$0，但是分子为$0$0。
因此，当$K!=0$K!=0时，$S(ω^k_n)=0$S(ωnk​)=0
那当$k=0$k=0时呢？
很显然，$S(ω^0_n)=n$S(ωn0​)=n
$c_k=∑_{j=0}^{n−1}a_j(∑_{i=0}^{n−1}(ω^{j−k}_n)^i)$ck​=j=0∑n−1​aj​(i=0∑n−1​(ωnj−k​)i)

当$j!=k$j!=k时，值为$0$0
当$j=k$j=k时，值为$n$n
$c_k=na_k$ck​=nak​

$a_k=\frac{c_k}{n}$ak​=nck​​
这样我们就得到点值与系数之间的表示了。

---

代码实现

差不多快完结了吧。。。（狂立flag）
从上面一大堆乱七八糟的东西，我们就得到了一个FFT的过程，系数->点值->系数。
引用一下远航之曲大佬的图

递归版

void FFT(int limit, Complex *a, int type) {
    if (limit == 1) return;
    Complex a1[limit >> 1], a2[limit >> 1];
    for (int i = 0; i <= limit; i += 2) a1[i >> 1] = a[i], a2[i >> 1] = a[i + 1];
    FFT(limit >> 1, a1, type); FFT(limit >> 1, a2, type); 
    Complex wn = Complex(cos(2.0 * Pi / limit), type * sin(2.0 * (Pi / limit))), w = Complex(1, 0);
    for (int i = 0; i < (limit >> 1); i ++, w = w * wn) a[i] = a1[i] + w * a2[i], a[i + (limit >> 1)] = a1[i] - w * a2[i];
}

Complex就是前置芝士里面我自己写的结构体。
但是以上这个东西，常数极大，递归本身的常数，还有三角函数所带的常数，到这个这个东西超级慢。
接下来开车了。。。

迭代版

天哪！还有一个这个东西，我讲的稍微快一点。
盗图

是不是有发现了什么东西，分治之后，我们的位置都是二进制翻转后的位置。
那么我们就可以先把这个序列的位置预处理好。然后直接可以对应合并。

 for (int i = 0; i < limit; i ++) r[i] = (r[i >> 1] >> 1)| ((i & 1) << (l - 1));

终极代码

#include <bits/stdc++.h>
#define ms(a, b) memset(a, b, sizeof(a))
#define ll long long
#define ull unsigned long long
#define ms(a, b) memset(a, b, sizeof(a))
#define inf 0x3f3f3f3f
#define db double
#define Pi acos(-1)
#define eps 1e-8
#define N 10000005
using namespace std;
template <typename T> void read(T &x) {
    x = 0; T fl = 1; char ch = 0;
    for (; ch < '0' || ch > '9'; ch = getchar()) if (ch == '-') fl = -1;
    for (; ch >= '0' && ch <= '9'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ 48);
    x *= fl;
}
template <typename T> void write(T x) {
    if (x < 0) x = -x, putchar('-');
    if (x > 9) write(x / 10); putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) { write(x); puts(""); }
struct Complex {
    db x, y;
    Complex (db xx = 0.0, db yy = 0.0) { x = xx, y = yy; }
    Complex operator + (Complex B) { return Complex(x + B.x, y + B.y); }
    Complex operator - (Complex B) { return Complex(x - B.x, y - B.y); }
    Complex operator * (Complex B) { return Complex(x * B.x - y * B.y, x * B.y + y * B.x); }
}a[N], b[N];
int n, m, limit, l;
int r[N];
void FFT(Complex *a, int type) {
    for (int i = 0; i < limit; i ++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 1; mid < limit; mid <<= 1) {
        Complex wn(cos(Pi / mid), type * sin(Pi / mid));
        for (int R = mid << 1, j = 0; j < limit; j += R) {
            Complex w(1, 0);
            for (int k = 0; k < mid; k ++, w = w * wn) {
                Complex x = a[j + k], y = w * a[j + mid + k];
                a[j + k] = x + y; a[j + mid + k] = x - y;
            }
        } 
    }
}
int main() {
    read(n); read(m);
    for (int i = 0; i <= n; i ++) { int x; read(x); a[i].x = 1.0 * x; }
    for (int i = 0; i <= m; i ++) { int x; read(x); b[i].x = 1.0 * x; }
    limit = 1; while (limit <= n + m) limit <<= 1, l ++;
    for (int i = 0; i < limit; i ++) r[i] = (r[i >> 1] >> 1)| ((i & 1) << (l - 1));
    FFT(a, 1); FFT(b, 1);
    for (int i = 0; i <= limit; i ++) a[i] = a[i] * b[i];
    FFT(a, -1); 
    for (int i= 0; i <= n + m; i ++) write((int)(a[i].x / limit + 0.5)), putchar(' ');
    return 0;
}

完结撒花wq11417个字