小学生都能看懂的FFT!!
作者(小学生)其实学这个学了两个月，但我相信你只要努力，就能成功
好的，废话不多说，正片开始
FFT章节1:了解FFT是干嘛的
oiwiki：FFT支持在 O ( n log n ) O(n\log n) O(nlogn) 的时间内计算两个 n n n次多项式的乘法，比朴素的 O ( n 2 ) O(n^2) O(n2) 算法更高效。由于两个整数的乘法也可以被当作多项式乘法，因此这个算法也可以用来加速大整数的乘法计算。
我自己的理解就是可以加速多项式的运算，比如多项式乘法，多项式求逆，多项式开根…，大部分都可以用FFT来计算
FFT章节2:复数
有些人可能会有点害怕，但是不要怕，这东西虽说是高中~大学的，但是理解起来还是很简单的
首先我们先要弄明白实数
实数是么子东西?实数是由有理数与无理数组成的，有理数又是什么?
有理数是只我们常见的整数，分数，有限小数小数，无限循环小数组成
有人会问:wxy,那无限不循环小数是什么呢?无限不循环小数是无理数，如： π e Φ . . . \pi\ e \ Φ... π e Φ...还有 2 3 . . . \sqrt{2}\ \sqrt{3}... 2 3 ...
如果你看不懂 π e Φ . . . \pi\ e \ Φ... π e Φ...还有 2 3 . . . \sqrt{2}\ \sqrt{3}... 2 3 ...没关系，本文中就只要 π \pi π 就行
π \pi π 是什么?是一个圆的 周长 直径 \frac{周长}{直径} 直径周长，约等于 3.1415926535.... 3.1415926535.... 3.1415926535....,当一个圆的直径为 1 1 1 时，这个圆的周长为 π \pi π
现在明白实数是什么了吧，那么虚数呢?
虚数，是建立在实数之外的一种数，他打破了不能给负数开根号的限制，根号是什么东西，简单来说对于 a > 0 a>0 a>0,都有 a × a = a \sqrt{a}\times\sqrt{a}=a a ×a =a,就是这么个意思，我们去想，在实数轴上面，是不是找不到一个数的平方等于 − 1 -1 −1,因为正数乘上正数等于正数，负数乘上负数也等于正数，根本找不到，所以说就有了虚数这么个东西
虚数: i = − 1 i=\sqrt{-1} i=−1 ,定义就是这么简单，我们称这个 i i i叫做 虚数 i 虚数i 虚数i
小公式（请自行理解）: − a = i a \sqrt{-a}=i\sqrt{a} −a =ia
实数的运算法则大多数都可运用道虚数上,那么复数又是什么呢
复数可以这么理解，复数有三大类：实数，虚数，实数+虚数
其实这三大类都可以用 a + b i a+bi a+bi这个式子来表示，其中 a , b a,b a,b均为实数，第一大类其实就是当 b = 0 b=0 b=0时的情况，第二大类其实就是 a = 0 a=0 a=0时的情况，第三大类其实就是 a , b a,b a,b均为实数的情况
因此 a + b i a+bi a+bi被我们称为虚数的表示方法，当然，还有别的表示方法，下面会讲
c++中有自带的虚数库，叫做 c o m p l e x complex complex，用法如下
#include<complex>//complex头文件
using namespace std;
...
complex<int>x;//定义一个实数部分和虚数部分均为整数的虚数，其实就是a+bi中的a和b都为int类型
x.real();//a+bi中的a
x.imag();//a+bi中的b
当然，我们也可以自己去定义，由于需要加减乘除的运算，这里我赘述一下
复数的加法
:
(
a
+
b
i
)
+
(
c
+
d
i
)
=
a
+
c
+
b
i
+
d
i
=
(
a
+
c
)
+
(
b
+
d
)
i
复数的减法
:
(
a
+
b
i
)
−
(
c
+
d
i
)
=
a
−
c
+
b
i
−
d
i
=
(
a
−
c
)
+
(
b
−
d
)
i
复数的乘法
:
(
a
+
b
i
)
(
c
+
d
i
)
=
a
c
+
a
d
i
+
c
b
i
+
b
d
i
2
因为
i
2
=
−
1
所以
a
c
+
a
d
i
+
c
b
i
+
b
d
i
2
=
(
a
c
−
b
d
)
+
(
a
d
+
c
b
)
i
复数的除法
a
+
b
i
c
+
d
i
=
.
.
.
.
=
a
c
+
b
d
c
2
+
d
2
+
b
c
−
a
d
c
2
+
d
2
i
复数的加法:\\ (a+bi)+(c+di)=a+c+bi+di=(a+c)+(b+d)i\\ 复数的减法:\\ (a+bi)-(c+di)=a-c+bi-di=(a-c)+(b-d)i\\ 复数的乘法:\\ (a+bi)(c+di)=ac+adi+cbi+bdi^2\\ 因为i^2=-1\\ 所以ac+adi+cbi+bdi^2=(ac-bd)+(ad+cb)i\\ 复数的除法\\ \frac{a+bi}{c+di}=....=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i
复数的加法:(a+bi)+(c+di)=a+c+bi+di=(a+c)+(b+d)i复数的减法:(a+bi)−(c+di)=a−c+bi−di=(a−c)+(b−d)i复数的乘法:(a+bi)(c+di)=ac+adi+cbi+bdi2因为i2=−1所以ac+adi+cbi+bdi2=(ac−bd)+(ad+cb)i复数的除法c+dia+bi=....=c2+d2ac+bd+c2+d2bc−adi
注意:在复数中是没有大小关系的
FFT章节3:多项式
了解多项式，我们先要了解单项式，单项式是未知数，数字的乘积，如: 5 x , π x 2 , 3 x 5x,\pi x^2,3x 5x,πx2,3x这种就是单项式
单项式的一般形式可以写为: a x k ax^k axk,其中 a a a为系数， k k k为次数
多项式就是有多个单项式相加减的来，两个次数相同并且未知数相同的单项式我们将他们称之为同类项
一个 n n n次多项式就是有 n n n个多项式相加减，其中，这些单项式的次数分别为 0 , 1 , 2 , 3 , 4 , 5 0,1,2,3,4,5 0,1,2,3,4,5,例如一个 5 5 5次多项式有： 3 + 5 x + 4 x 2 + 2 x 3 + x 4 + x 5 3+5x+4x^2+2x^3+x^4+x^5 3+5x+4x2+2x3+x4+x5, 3 3 3是常数项
一个n次多项式的一般形式为: ∑ i = 0 n a i x i \sum_{i=0}^{n}a_ix^i ∑i=0naixi
其中 a i a_i ai是 i i i次项的系数。 i i i为 i i i次项的系数
当然，我们通常使用函数来表示多项式，比如说多项式 A ( x ) = ∑ i = 0 n a i x i A(x)=\sum_{i=0}^{n}a_ix^i A(x)=∑i=0naixi，我们称这个函数为多项式 A A A
那么为什么要用函数来表示呢，比如说我要求出多项式 A ( x ) = x 2 + 2 x − 3 A(x)=x^2+2x-3 A(x)=x2+2x−3在 x = 2 x=2 x=2时的值，我们就相当于求 A ( 2 ) A(2) A(2)的值，只要把 A A A的式子展开，然后将式子中的 x x x给替换成 2 2 2就行了
FFT章节4:单位根
刚才我们学习了复数，现在我们来讲单位根是什么
要了解单位根，首先要知道什么是复平面
复平面简单来说就是一个平面直角坐标系，只不过 y y y轴上面的不再是 1 , 2 , 3 , 4 , 5... 1,2,3,4,5... 1,2,3,4,5...了，而是 i , 2 i , 3 i , 4 i , 5 i . . . . i,2i,3i,4i,5i.... i,2i,3i,4i,5i....。
在这个复平面上，我们可以表示出复数，通俗来讲，一条数轴可以表示出实数，两条数轴可以表示出虚数和实数，一个复平面可以表示出复数!!!
图解:
在复平面上以
(
0
,
0
)
(0,0)
(0,0)为圆心，祖宗十八代为半径，
1
1
1为半径，画一个圆，这个圆就被我们称为单位圆
单位根章节分块:弧度制
弧度制也就是在单位圆上表示角度的一种方法，比如说，我们要说 3 0 ∘ 30^\circ 30∘，我会改成 π 6 \frac{\pi}{6} 6π，这是怎么改过来的呢，其实有一个很好记的方法，因为单位圆的周长是 2 π 2\pi 2π,所以我们也可以认为 2 π 就是 36 0 ∘ 2\pi就是360^\circ 2π就是360∘，则 π 也就是 18 0 ∘ \pi也就是180^\circ π也就是180∘，为了简便，我就写成 π = 18 0 ∘ \pi=180^\circ π=180∘
弧度制的意义是什么。我们可以这么理解: π = 18 0 ∘ \pi=180^\circ π=180∘的意思相当于在单位圆上从 ( 1 , 0 ) (1,0) (1,0)这个点逆时针旋转 π \pi π，所得的点和原来点组成的角的夹角就是 18 0 ∘ 180^\circ 180∘
接下来我就可以开始讲第二中复数的表示方法了
著名的欧拉公式就是 e i x = cos x + i sin x e^{ix}=\cos x+i\sin x eix=cosx+isinx，这个式子在复平面上是很好理解的，其中 x x x是弧度制数，在单位圆上旋转 x x x，得到的点的横坐标正好就是 cos x \cos x cosx，纵坐标也正好是 sin x \sin x sinx，由于这个是复平面，所以纵坐标 sin x \sin x sinx要乘上一个 i i i，所以得出结论 e i x e^{ix} eix的意义就是从 ( 1 , 0 ) (1,0) (1,0)这个点沿着单位根旋转 x x x,得到的那一个点(其实就是数)
欧拉公式就是这么来的 e i π = cos π + i sin π = − 1 e^{i\pi}=\cos \pi+i\sin \pi=-1 eiπ=cosπ+isinπ=−1，也就是从 ( 1 , 0 ) (1,0) (1,0)这个点沿着单位根旋转 π \pi π,得到的那一个点,正好落在 ( − 1 , 0 ) (-1,0) (−1,0)上
于是，一个复数就可以被表示为 r e i x re^{ix} reix， e i x e^{ix} eix用来确定角度， r r r来确定从 ( 0 , 0 ) (0,0) (0,0)伸出去的长度
那么现在可以开始讲单位根是什么了吧
单位根就相当于把单位圆给分成了 n n n等分，取其中的 1 1 1份
算了说出来吧， w n w_n wn表示将单位单位圆给分成了 n n n等分，其中的 1 1 1份，他的公式为
w n = e i 2 π n w_n=e^{i\frac{2\pi}{n}} wn=ein2π
，那么取其中的两份并不是 2 w n 2w^n 2wn,而是 w n 2 w_n^2 wn2，也就是说n等分分别是 w n 0 , w n 1 , w n 2 , w n 3 , w n 4 , w n 5 . . . . . w_n^0,w_n^1,w_n^2,w_n^3,w_n^4,w_n^5..... wn0,wn1,wn2,wn3,wn4,wn5.....，如下图，讲单位圆分成八等份，那个右上角的那个红点就是 w 8 w_8 w8,最上面那个点就是 w 8 2 w_8^2 w82，以此类推
公式:
w r n r k = e i 2 π k r n r = e i 2 π k n = w n k w_{rn}^{rk}=e^{i\frac{2\pi kr}{nr}}=e^{i\frac{2\pi k}{n}}=w_n^k wrnrk=einr2πkr=ein2πk=wnk
w n n = w n 0 = 1 w_n^{n}=w_n^0=1 wnn=wn0=1(很显然，看图都可以看出来)
w n k + n 2 = − w n k w_n^{k+\frac{n}{2}}=-w_n^k wnk+2n=−wnk(这个看图也可以看出来)
你以为现在已经很难了吗，不！后面更难
(正片开始)FFT章节5:FFT正变换的原理及代码实现
首先我们要先学习一下DFT的原理
我们要将 A A A和 B B B相乘，普通的步骤为，直接相乘，次数相加对吧，写一个伪代码
定义两个整数n,m，表示多项式A和B的项数（不是次数）
定义三个数组A,B,C,分别表示多项式A的系数，多项式B的系数，多项式A和多项式B的乘积的系数
for 0~n-1 cin>>A[i]
for 0~m-1 cin>>B[i]
for i to 1~n-1
for j to 1~m-1
C[i+j]+=A[i]*B[j]//两个单项式相乘为他们的次数相加，系数相乘
这样的时间复杂度是 O ( n 2 ) O(n^2) O(n2)的
我们知道，可以用 n 个点来表示一个 n − 1 次多项式，这个我在我之前的帖子真 . 浅析拉格朗日插值法里面提到过 n个点来表示一个n-1次多项式，这个我在我之前的帖子真.浅析拉格朗日插值法里面提到过 n个点来表示一个n−1次多项式，这个我在我之前的帖子真.浅析拉格朗日插值法里面提到过
FFT的步骤是这样的:
多项式
A
和
B
−
>
将多项式转化为多个点
−
>
将每个点的因变量给乘起来，得到多项式
A
×
B
的多个点
−
>
通过高斯消元得到多项式
A
×
B
多项式A和B->将多项式转化为多个点->将每个点的因变量给乘起来，得到多项式A\times B的多个点->通过高斯消元得到多项式A\times B
多项式A和B−>将多项式转化为多个点−>将每个点的因变量给乘起来，得到多项式A×B的多个点−>通过高斯消元得到多项式A×B
如果单纯使用高斯消元和多点求值，时间复杂度还是
O
(
n
3
)
O(n^3)
O(n3)，虽说有
O
(
n
log
2
n
)
O(n\log^2 n)
O(nlog2n)的多点求值，但是高斯消元依然是
O
(
n
3
)
O(n^3)
O(n3)，我就不建议大家用这种方法
FFT是怎么将多项式转化为多个点的呢?
比如说DFT要将多项式
A
(
x
)
=
x
2
+
2
x
−
3
A(x)=x^2+2x-3
A(x)=x2+2x−3，我们要将他转换为
4
4
4个点，也就是我们自己要选择四个数，然后把这四个数分别带入到
A
(
x
)
A(x)
A(x)中，求出来值，暴力的想法是随机四个值，然后求出这些值在
A
(
x
)
A(x)
A(x)中的值，时间复杂度
O
(
n
2
)
O(n^2)
O(n2)，我们不妨换一种想法，取两个多项式分别是
A
0
(
x
)
和
A
1
(
x
)
A_0(x)和A_1(x)
A0(x)和A1(x),
A
0
(
x
)
A_0(x)
A0(x)中存储的是
A
(
x
)
A(x)
A(x)的偶数次项，
A
1
(
x
)
A_1(x)
A1(x)中存储的是
A
(
x
)
A(x)
A(x)的技术次项
A
(
x
)
=
x
2
+
2
x
−
3
A
0
(
x
2
)
=
x
2
−
3
A
1
(
x
2
)
=
2
x
A
0
(
x
)
=
x
−
3
A
1
(
x
)
=
2
A
(
x
)
=
A
0
(
x
2
)
+
x
A
1
(
x
2
)
A(x)=x^2+2x-3\\ A_0(x^2)=x^2-3\\ A_1(x^2)=2x\\ A_0(x)=x-3\\ A_1(x)=2\\ A(x)=A_0(x^2)+xA_1(x^2)
A(x)=x2+2x−3A0(x2)=x2−3A1(x2)=2xA0(x)=x−3A1(x)=2A(x)=A0(x2)+xA1(x2)
为什么要这么做呢?因为这样一来，
A
0
(
x
2
)
A_0(x^2)
A0(x2)就是一个偶函数，
A
1
(
x
)
A_1(x)
A1(x)就是一个奇函数，我们选的点值只要是一正一负的就行了
意思就是:
A
(
x
)
=
A
0
(
x
2
)
+
x
A
1
(
x
2
)
A
(
−
x
)
=
A
0
(
x
2
)
−
x
A
1
(
x
2
)
A(x)=A_0(x^2)+xA_1(x^2)\\ A(-x)=A_0(x^2)-xA_1(x^2)
A(x)=A0(x2)+xA1(x2)A(−x)=A0(x2)−xA1(x2)
则这的规模只有原来的一半，并且求解
A
0
(
x
)
和
A
1
(
x
)
A_0(x)和A_1(x)
A0(x)和A1(x)也可以使用这个方法，不过在进行之前，我们要将项数
n
,
m
n,m
n,m给变成一个比他们两个的和都要大并且还是
2
2
2的正整数次幂的数，放心，因为它既时项数变大了，但是他的那些高次项的系数依然是0，毫无差别，只不过你的数组要开打
3
3
3倍
定义两个整数n,m，表示多项式A和B的项数（不是次数）
int l;//为多项式的2的正整数次幂的项数
while(l<=n+m)l<<=1;
你以为这就行了嘛?不，这种方法有漏洞，你们尽然没看出来?
我们取一正一负，但是
A
0
,
A
1
A_0,A_1
A0,A1取得是
x
2
x^2
x2，在实数范围内，
x
2
x^2
x2一直是非负整数，所以就取不了，那么什么数可以在平方后还能保证一正一负呢，答案是单位根，我们使用单位根代替
x
x
x来看看
A
(
w
n
k
)
=
A
0
(
w
n
2
k
)
+
w
n
k
A
1
(
x
n
2
k
)
=
A
0
(
w
n
2
k
)
+
w
n
k
A
1
(
x
n
2
k
)
A
(
w
n
k
+
n
2
)
=
A
0
(
w
n
2
k
+
n
)
+
w
n
k
+
n
2
A
1
(
x
n
2
k
+
n
)
=
A
0
(
w
n
2
k
)
−
w
n
k
A
1
(
x
n
2
k
)
A(w_n^k)=A_0(w_n^{2k})+w_n^kA_1(x_n^{2k})\\ =A_0(w_\frac{n}{2}^{k})+w_n^kA_1(x_\frac{n}{2}^{k})\\ A(w_n^{k+\frac{n}{2}})=A_0(w_n^{2k+n})+w_n^{k+\frac{n}{2}}A_1(x_n^{2k+n})\\ =A_0(w_\frac{n}{2}^{k})-w_n^kA_1(x_\frac{n}{2}^{k})
A(wnk)=A0(wn2k)+wnkA1(xn2k)=A0(w2nk)+wnkA1(x2nk)A(wnk+2n)=A0(wn2k+n)+wnk+2nA1(xn2k+n)=A0(w2nk)−wnkA1(x2nk)
并且
A
0
,
A
1
A_0,A_1
A0,A1的长度正好是
n
2
\frac{n}{2}
2n，所以说我们可以将
A
0
,
A
1
A_0,A_1
A0,A1也这样处理就没问题了
注意边界:当项数为 1 1 1时，直接return，因为只剩下常数项了，无论多项式里是什么数答案都是那个常数项
时间复杂度: O ( n log n ) O(n\log n) O(nlogn)
#define PI acos(1.0)
//使用反三角函数求PI
void FFT(complex<double>a[],int n){//多项式和他的项数
if(n==1)return;
complex<double>a0[n/2+1],a1[n/2+1];//两个多项式A0和A1
for(int i=0;i<n;i+=2){//注意，这里必须是<n,因为n是项数
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}//将系数装填进去，理解不了自己在纸上画一画
FFT(a0,n/2),FFT(a1,n/2);//递归求解
complex<double>w(cos(2*PI/n),sin(2*PI/n)),W(1,0);//根据单位根的定义创建单位根
for(int i=0;i<n/2;i++){
a[i]=a0[i]+W*a1[i];
a[i+n/2]=a0[i]-W*a1[i];
W*=w;//W为w_n^i
}
}
FFT章节6:FFT逆变换，IFFT的原理及实现
我们已经得到了 A × B A\times B A×B的点值表示，我们现在要通过他的点来得到他本身
简化一下，意思就是你已经知道了多项式 C C C的点值表示，我们要求出他的系数表示
一个简单的想法，我们知道 C C C的系数，我们只要把他的点值表示带入高斯消元就可以了，具体怎么干去看我的浅析.拉格朗日插值，但是这样的时间复杂度为 O ( n 3 ) O(n^3) O(n3)，是在是不好用
那么傅里叶是怎么干的呢?(注意:以下用了许多
LaTeX
\LaTeX
LATEX，请读者务必反复阅读，如果真读不懂那就直接看结论算了)
我们将多项式
C
的点值表示也看作一个
n
−
1
次多项式，记作
D
将
D
也进行一次
F
F
T
，得到多项式
F
，只不过这次的
F
F
T
的单位根不再是
w
n
k
了，是
w
n
−
k
则
:
F
k
=
∑
i
=
0
n
−
1
D
i
(
w
n
−
k
)
i
=
∑
i
=
0
n
−
1
(
w
n
−
k
)
i
∑
j
=
0
n
−
1
C
j
(
w
n
i
)
j
=
∑
j
=
0
n
−
1
C
j
∑
i
=
0
n
−
1
(
w
n
−
k
)
i
(
w
n
i
)
j
=
∑
j
=
0
n
−
1
C
j
∑
i
=
0
n
−
1
(
w
n
i
)
−
k
(
w
n
i
)
j
=
∑
j
=
0
n
−
1
C
j
∑
i
=
0
n
−
1
(
w
n
i
)
j
−
k
那么这个东西我们怎么对其进行简化呢
?
我们设函数
G
(
n
,
j
,
k
)
=
∑
i
=
0
n
−
1
(
w
n
i
)
j
−
k
那么当
j
=
k
时，
G
(
n
,
j
,
k
)
明显等于
n
当
j
≠
k
时
,
G
(
n
,
j
,
k
)
经过等比数列求和公式得
:
G
(
n
,
j
,
k
)
=
(
w
n
j
−
k
)
n
−
1
w
n
j
−
k
−
1
=
(
w
n
n
)
j
−
k
−
1
w
n
j
−
k
−
1
=
1
−
1
w
n
j
−
k
−
1
=
0
由此可得
:
G
(
n
,
j
,
k
)
=
n
[
j
=
k
]
意思为当
j
=
k
时这个函数等于
n
，在其他情况下等于
0
所以
:
F
k
=
∑
j
=
0
n
−
1
C
j
n
[
j
=
k
]
F
k
=
n
C
k
我们就可以得到，
C
k
=
F
k
n
我们将多项式C的点值表示也看作一个n-1次多项式，记作D\\ 将D也进行一次FFT，得到多项式F，只不过这次的FFT的单位根不再是w_n^k了，是w_n^{-k}\\ 则:\\ F_k=\sum_{i=0}^{n-1}D_i(w_n^{-k})^i\\ =\sum_{i=0}^{n-1}(w_n^{-k})^i\sum_{j=0}^{n-1}C_j(w_n^i)^j\\ =\sum_{j=0}^{n-1}C_j\sum_{i=0}^{n-1}(w_n^{-k})^i(w_n^i)^j\\ =\sum_{j=0}^{n-1}C_j\sum_{i=0}^{n-1}(w_n^i)^{-k}(w_n^i)^j\\ =\sum_{j=0}^{n-1}C_j\sum_{i=0}^{n-1}(w_n^i)^{j-k}\\ 那么这个东西我们怎么对其进行简化呢?\\ 我们设函数G(n,j,k)=\sum_{i=0}^{n-1}(w_n^i)^{j-k}\\ 那么当j=k时，G(n,j,k)明显等于n\\ 当j≠k时,G(n,j,k)经过等比数列求和公式得:\\ G(n,j,k)=\frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}\\ =\frac{(w_n^n)^{j-k}-1}{w_n^{j-k}-1}\\ =\frac{1-1}{w_n^{j-k}-1}\\ =0\\ 由此可得:G(n,j,k)=n[j=k]意思为当j=k时这个函数等于n，在其他情况下等于0\\ 所以:\\ F_k=\sum_{j=0}^{n-1}C_jn[j=k]\\ F_k=nC_k\\ 我们就可以得到，C_k=\frac{F_k}{n}
我们将多项式C的点值表示也看作一个n−1次多项式，记作D将D也进行一次FFT，得到多项式F，只不过这次的FFT的单位根不再是wnk了，是wn−k则:Fk=i=0∑n−1Di(wn−k)i=i=0∑n−1(wn−k)ij=0∑n−1Cj(wni)j=j=0∑n−1Cji=0∑n−1(wn−k)i(wni)j=j=0∑n−1Cji=0∑n−1(wni)−k(wni)j=j=0∑n−1Cji=0∑n−1(wni)j−k那么这个东西我们怎么对其进行简化呢?我们设函数G(n,j,k)=i=0∑n−1(wni)j−k那么当j=k时，G(n,j,k)明显等于n当j=k时,G(n,j,k)经过等比数列求和公式得:G(n,j,k)=wnj−k−1(wnj−k)n−1=wnj−k−1(wnn)j−k−1=wnj−k−11−1=0由此可得:G(n,j,k)=n[j=k]意思为当j=k时这个函数等于n，在其他情况下等于0所以:Fk=j=0∑n−1Cjn[j=k]Fk=nCk我们就可以得到，Ck=nFk
所以，我们只需要得到
F
F
F我们就只需要将它所有的系数全都除以一个
n
n
n我们就得到了
C
C
C
F F F怎么求就不用我多说了吧，在原本的FFT代码上改就行了
#define PI acos(1.0)
//使用反三角函数求PI
void FFT(complex<double>a[],int n,int inv){//多项式和他的项数
if(n==1)return;
complex<double>a0[n/2+1],a1[n/2+1];//两个多项式A0和A1
for(int i=0;i<n;i+=2){//注意，这里必须是<n,因为n是项数
a0[i/2]=a[i];
a1[i/2]=a[i+1];
}//将系数装填进去，理解不了自己在纸上画一画
FFT(a0,n/2),FFT(a1,n/2);//递归求解
complex<double>w(cos(inv*2*PI/n),sin(inv*2*PI/n)),W(1,0);//根据单位根的定义创建单位根，当inv=1时为正变换，当inv=-1时为逆变换
//complex<double>w(cos(2*PI/n),inv*sin(2*PI/n)),W(1,0);可以这么写，因为cos(-a)=cos(a),sin(-a)=-sin(a)
for(int i=0;i<n/2;i++){
a[i]=a0[i]+W*a1[i];
a[i+n/2]=a0[i]-W*a1[i];
W*=w;//W为w_n^i
}
}
好了，现在就是FFT的代码了，时间复杂度为 O ( n log n ) O(n\log n) O(nlogn)，空间复杂度为 O ( n log n ) O(n\log n) O(nlogn)
什么，你问我怎么用，那你肯定是没听前面的，赶紧去听!!!
对了最后输出时要四舍五入才是正确的哦
FFT章节7:（FFT优化）非递归版FFT
刚才我们讲了递归版的FFT，这种FFT可能会然我们递归的次数变多，从而爆栈
并且上面的递归般的做法的空间复杂度也是 O ( n log n ) O(n\log n) O(nlogn)的，也不好用
于是我们便可以使用非递归版FFT
我们先来解决爆栈的问题
归并排序原本是一种分治递归版的做法，但是他可以进行改进
归并排序的非递归版就是一层一层的往上合并，那么我们FFT可以这么做吗?
答案是:可以，只不过难一点
我们FFT要分治先要将这个多项式按奇偶项分开，然后合并
我们可以找一下规律
长度:8
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
什么?你没看出有什么规律?没关系，我来告诉你规律吧
首先，FFT先将长度给定成了二的整数次幂
所以我们观察，8的二进制下有3个0，而这些序列中的数的二进制不超过3位
他们改完过后的位置都在自己二进制反转后的位置上
什么?尼稳我为顺么?沃野布吉岛
画个图更理解:
(000)(001)(010)(011)(100)(101)(110)(111)
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
(000)(100)(010)(110)(001)(101)(011)(111)
上面的二进制和下面的二进制是反的
于是我们便可以先将他们变成这样，然后再合并
那么怎么变成这样呢?
朴素 O ( log n ) O(\log n) O(logn):
int zhuan(int l,int x){//l是0的位数
int sum=0;
for(int i=0;i<=l;i++){
if(x>>i&1){
sum|=(1<<(l-i));
}
}
return sum;
}
那么总共有 O ( n ) O(n) O(n)个数，时间复杂度是 O ( n log n ) O(n\log n) O(nlogn)的
那么怎么优化呢?我们可以用到递推
递推 预处理 O ( n ) ，使用 O ( 1 ) 预处理O(n)，使用O(1) 预处理O(n)，使用O(1)
rev[0]=0;
for(int i=1;i<l;i++){//注意，这里的l是序列长度
//rev[i]代表这是i二进制转换后的数值
rev[i]=rev[i>>1]>>1;
//先将前面l-1位转换了
if(i&1){//再将后面那位转换
rev[i]|=(l>>1);
}
//写简便一点
//rev[i]=(rev[i>>1]>>1)|((i&1)*(l>>1));
}
这样我们就将rev处理完了
这种优化我们成为蝴蝶优化
我们现在处理一下第二个优化
空间优化
我们发现，由于FFT每次修改 a i a_i ai使都要使用到之前的 a i a_i ai，于是我们将之前的 a i a_i ai不用数组保存，只需要用两个变量来保存就可以了
最后，我们将FFT依次合并上去就好了:D
我们处理完之后的FFT就编程了这样:
void FFT(cp* a,int* rev,int n,int inv){
for(int i=0;i<n;i++){
if(rev[i]<i)swap(a[i],a[rev[i]]);//如果没被反转过就反转
}
for(int g=1;g<n;g<<=1){//g是递归的长度/2
cp Omiga=cp(cos(PI/g),inv*sin(PI/g));//
for(int i=0;i<n;i+=(g<<1)){
cp omiga=cp(1,0);
for(int j=0;j<g;j++,omiga*=Omiga){
cp oo=omiga*a[i+j+g];
cp aa=a[i+j];
a[i+j]=aa+oo;
a[i+j+g]=aa-oo;
}
}
}
}
注意:接下来主要讲优化
FFT章节8：(FFT优化)三次变两次
我们FFT不是需要先将 A , B A,B A,B两个多项式先都进行一次 F F T FFT FFT，然后在对 A × B A\times B A×B进行一次 I F F T IFFT IFFT
我们怎么进行优化呢?
我们可以将多项式 A A A的虚数部分给换成多项式 B B B，然后只对他们进行一次FFT，然后求出他们的平方，随后在返回来就好了
那么此时答案在哪里呢?
其实复数的运算对多项式也有效
(
A
+
B
i
)
2
=
(
A
2
−
B
2
)
+
i
(
2
A
B
)
其实复数的运算对多项式也有效\\ (A+Bi)^2\\ =(A^2-B^2)+i(2AB)
其实复数的运算对多项式也有效(A+Bi)2=(A2−B2)+i(2AB)
所以，答案的两倍就是这个进行FFT后的多项式的虚数部分
FFT到这里就完结了，喜欢的小伙伴点个赞，其他优化FFT的方法效果都不高，除非你要大力卡常，下一期，我们将要讲FFT的亲兄弟–NTT的实现方法，不会丢精度哟
对了，今天给大家留两道作业题: