FFT前言
快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。
FFT（Fast Fourier Transformation） 是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
——百度百科
FFT（Fast Fourier Transformation），中文名快速傅里叶变换，用来加速多项式乘法
朴素高精度乘法时间[Math Processing Error] O(n^2)O(n
2
)，但[Math Processing Error] FFTFFT能[Math Processing Error] O(nlog_2 n)O(nlog
2
n)的时间解决
[Math Processing Error] FFTFFT名字逼格高，也难懂，其他教程写得让人看不太懂，于是自己随便写一下
建议对复数、三角函数相关知识有所耳闻 （不会也无所谓）
下面难懂的点我会从网上盗
多项式的系数表示法和点值表示法
[Math Processing Error] FFTFFT其实是一个用[Math Processing Error] O(nlog_2n)O(nlog
2
n)的时间将一个用系数表示的多项式转换成它的点值表示的算法
多项式的系数表示和点值表示可以互相转换
系数表示法
一个n-1次n项多项式[Math Processing Error] f(x)f(x)可以表示为[Math Processing Error] f(x)=sum^{n-1}_{i=0}a_ix^if(x)=∑
i=0
n−1
a
i
x
i
也可以用每一项的系数来表示[Math Processing Error] f(x)f(x)，即
[Math Processing Error] f(x)={a_0,a_1,a_2,...,a_{n-1} }
f(x)={a
0
,a
1
,a
2
,...,a
n−1
}
这就是系数表示法，也就是平时数学课上用的方法
点值表示法
把多项式放到平面直角坐标系里面，看成一个函数
把[Math Processing Error] nn个不同的[Math Processing Error] xx代入，会得出[Math Processing Error] nn个不同的[Math Processing Error] yy，在坐标系内就是[Math Processing Error] nn个不同的点
那么这[Math Processing Error] nn个点唯一确定该多项式，也就是有且仅有一个多项式满足[Math Processing Error] ∀k,f(x_k)=y_k∀k,f(x
k
)=y
k
理由很简单，把[Math Processing Error] nn条式子联立起来成为一个有n条方程的n元方程组，每一项的系数都可以解出来
那么[Math Processing Error] f(x)f(x)还可以用[Math Processing Error] f(x)={(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))}f(x)={(x
0
,f(x
0
)),(x
1
,f(x
1
)),(x
2
,f(x
2
)),...,(x
n−1
,f(x
n−1
))}来表示
这就是点值表示法
高精度乘法下两种多项式表示法的区别
对于两个用系数表示的多项式，我们把它们相乘
设两个多项式分别为[Math Processing Error] A(x),B(x)A(x),B(x)
我们要枚举[Math Processing Error] AA每一位的系数与[Math Processing Error] BB每一位的系数相乘
那么系数表示法做多项式乘法时间复杂度[Math Processing Error] O(n^2)O(n
2
)
但两个用点值表示的多项式相乘，只需要[Math Processing Error] O(n)O(n)的时间
什么意思呢？
设两个点值多项式分别为
[Math Processing Error] f(x)={(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))}
f(x)={(x
0
,f(x
0
)),(x
1
,f(x
1
)),(x
2
,f(x
2
)),...,(x
n−1
,f(x
n−1
))}
[Math Processing Error] g(x)={(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),...,(x_{n-1},g(x_{n-1}))}
g(x)={(x
0
,g(x
0
)),(x
1
,g(x
1
)),(x
2
,g(x
2
)),...,(x
n−1
,g(x
n−1
))}
设它们的乘积是[Math Processing Error] h(x)h(x)，那么
[Math Processing Error] h(x)={(x_0,f(x_0)·g(x_0)),(x_1,f(x_1)·g(x_1)),...,(x_{n-1},f(x_{n-1})·g(x_{n-1}))}
h(x)={(x
0
,f(x
0
)⋅g(x
0
)),(x
1
,f(x
1
)⋅g(x
1
)),...,(x
n−1
,f(x
n−1
)⋅g(x
n−1
))}
所以这里的时间复杂度只有一个枚举的[Math Processing Error] O(n)O(n)
突然感觉高精度乘法能[Math Processing Error] O(n)O(n)暴艹一堆题？
但是朴素的系数表示法转点值表示法的算法还是[Math Processing Error] O(n^2)O(n
2
)的，逆操作类似
朴素系数转点值的算法叫DFT（离散傅里叶变换），点值转系数叫IDFT（离散傅里叶逆变换）
难道高精度乘法只能[Math Processing Error] O(n^2)O(n
2
)了吗？
DFT前置知识&技能
复数
毕竟高中有所以不多说
我们把形如a+bi（a,b均为实数）的数称为复数，其中a称为实部，b称为虚部，i称为虚数单位。当虚部等于零时，这个复数可以视为实数；当z的虚部不等于零时，实部等于零时，常称z为纯虚数。复数域是实数域的代数闭包，也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入，经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作，此概念逐渐为数学家所接受。
——百度百科
初中数学老师会告诉你没有[Math Processing Error] sqrt{-1}
−1
，但仅限[Math Processing Error] RR
扩展至复数集[Math Processing Error] CC，定义[Math Processing Error] i^2=-1i
2
=−1，一个复数[Math Processing Error] zz可以表示为[Math Processing Error] z=a+bi(a,bin R)z=a+bi(a,b∈R)
其中[Math Processing Error] aa称为实部，[Math Processing Error] bb称为虚部，[Math Processing Error] ii称为虚数单位
在复数集中就可以用[Math Processing Error] ii表示负数的平方根，如[Math Processing Error] sqrt{-7}=sqrt{7}i
−7
=
7
i
还可以把复数看成复平面直角坐标系上的一个点，比如下面
[Math Processing Error] xx轴就是实数集中的坐标轴，[Math Processing Error] yy轴就是虚数单位[Math Processing Error] ii轴
这个点[Math Processing Error] (2,3)(2,3)表示的复数就是[Math Processing Error] 2+3i2+3i，或者想象它代表的向量为[Math Processing Error] (2,3)(2,3)
其实我们还可以自己想象 （其实没有这种表达方式） 它可以表示为[Math Processing Error] (sqrt{13}, heta)(
13
,θ)
一个复数[Math Processing Error] zz的模定义为它到原点的距离，记为[Math Processing Error] |z|=sqrt{a^2+b^2}∣z∣=
a
2
+b
2
一个复数[Math Processing Error] z=a+biz=a+bi的共轭复数为[Math Processing Error] a-bia−bi（虚部取反），记为[Math Processing Error] overline{z}=a-bi
z
=a−bi
复数的运算
复数不像点或向量，它和实数一样可以进行四则运算
设两个复数分别为[Math Processing Error] z_1=a+bi,z_2=c+diz
1
=a+bi,z
2
=c+di，那么
[Math Processing Error] z_1+z_2=(a+c)+(b+d)i
z
1
+z
2
=(a+c)+(b+d)i
[Math Processing Error] z_1z_2=(ac−bd)+(ad+bc)i
z
1
z
2
=(ac−bd)+(ad+bc)i
复数相加也满足平行四边形法则
这张是从网上盗的
即[Math Processing Error] AB+AD=ACAB+AD=AC
复数相乘还有一个值得注意的小性质
[Math Processing Error] (a_1, heta_1)*(a_2, heta_2)=(a_1a_2, heta_1+ heta_2)
(a
1
,θ
1
)∗(a
2
,θ
2
)=(a
1
a
2
,θ
1
+θ
2
)
即模长相乘，极角相加
DFT（离散傅里叶变换）
一定注意从这里开始所有的[Math Processing Error] nn都默认为[Math Processing Error] 22的整数次幂
对于任意系数多项式转点值，当然可以随便取任意[Math Processing Error] nn个[Math Processing Error] xx值代入计算
但是暴力计算[Math Processing Error] x_k^0,x_k^1,...,x_k^{n-1}(kin[0,n))x
k
0
,x
k
1
,...,x
k
n−1
(k∈[0,n))当然是[Math Processing Error] O(n^2)O(n
2
)的时间
其实可以代入一组神奇的[Math Processing Error] xx，代入以后不用做那么多的次方运算
这些[Math Processing Error] xx当然不是乱取的，而且取这些[Math Processing Error] xx值应该就是 傅里叶 的主意了
考虑一下，如果我们代入一些[Math Processing Error] xx，使每个[Math Processing Error] xx的若干次方等于[Math Processing Error] 11，我们就不用做全部的次方运算了
[Math Processing Error] ±1±1是可以的，考虑虚数的话[Math Processing Error] ±i±i也可以，但只有这四个数远远不够
傅里叶说：这个圆圈上面的点都可以做到
以原点为圆心，画一个半径为[Math Processing Error] 11的单位圆
那么单位圆上所有的点都可以经过若干次次方得到[Math Processing Error] 11
傅里叶说还要把它给[Math Processing Error] nn等分了，比如此时[Math Processing Error] n=8n=8
橙色点即为[Math Processing Error] n=8n=8时要取的点，从[Math Processing Error] (1,0)(1,0)点开始，逆时针从[Math Processing Error] 00号开始标号，标到[Math Processing Error] 77号
记编号为[Math Processing Error] kk的点代表的复数值为[Math Processing Error] omega_n^kω
n
k
，那么由模长相乘，极角相加可知[Math Processing Error] (omega_n^1)^k=omega_n^k(ω
n
1
)
k
=ω
n
k
其中[Math Processing Error] omega_n^1ω
n
1
称为[Math Processing Error] nn次单位根，而且每一个[Math Processing Error] omegaω都可以求出 （我三角函数不好）
[Math Processing Error] omega_n^k=cos{kover n}2π+isin{kover n} 2π
ω
n
k
=cos
n
k
2π+isin
n
k
2π
或者说也可以这样解释这条式子
注意[Math Processing Error] sin^2 heta+cos^2 heta=1sin
2
θ+cos
2
θ=1什么的，就容易理解了
那么[Math Processing Error] omega^0_n,omega^1_n,...,omega^{n-1}_nω
n
0
,ω
n
1
,...,ω
n
n−1
即为我们要代入的[Math Processing Error] x_0,x_1,...,x_{n-1}x
0
,x
1
,...,x
n−1
单位根的一些性质
推[Math Processing Error] FFTFFT的过程中需要用到[Math Processing Error] omegaω的一些性质
[Math Processing Error] omega^k_n=omega^{2k}_{2n}ω
n
k
=ω
2n
2k
它们表示的点（或向量）表示的复数是相同的
证明
[Math Processing Error] omega^k_n=cos{kover n}2π+isin{kover n} 2π=cos{2kover 2n}2π+isin{2kover 2n} 2π=omega^{2k}_{2n}ω
n
k
=cos
n
k
2π+isin
n
k
2π=cos
2n
2k
2π+isin
2n
2k
2π=ω
2n
2k
[Math Processing Error] omega^{k+{n over 2}}_n=-omega_n^kω
n
k+
2
n
=−ω
n
k
它们表示的点关于原点对称，所表示的复数实部相反，所表示的向量等大反向
证明
[Math Processing Error] omega^{nover 2}_n=cos{{nover 2}over n}2pi+isin{{nover 2}over n}2pi=cospi+isinpi=-1ω
n
2
n
=cos
n
2
n
2π+isin
n
2
n
2π=cosπ+isinπ=−1
（这个东西和[Math Processing Error] e^{ix}=cosx+isinxe
ix
=cosx+isinx与[Math Processing Error] e^{ipi}+1=0e
iπ
+1=0有点关系，我不会就不讲了）
[Math Processing Error] omega^0_n=omega^n_nω
n
0
=ω
n
n
都等于[Math Processing Error] 11，或[Math Processing Error] 1+0i1+0i
FFT（快速傅里叶变换）
虽然[Math Processing Error] DFTDFT搞出来一堆很牛逼的[Math Processing Error] omegaω作为代入多项式的[Math Processing Error] xx值
但只是代入这类特殊[Math Processing Error] xx值法的变换叫做[Math Processing Error] DFTDFT而已，还是要代入单位根暴力计算
DFT还是暴力[Math Processing Error] O(n^2)O(n
2
)……
但[Math Processing Error] DFTDFT可以分治来做，于是 FFT（快速傅里叶变换） 就出来了
首先设一个多项式[Math Processing Error] A(x)A(x)
[Math Processing Error] A(x)=sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}
A(x)=
i=0
∑
n−1
a
i
x
i
=a
0
+a
1
x+a
2
x
2
+...+a
n−1
x
n−1
按[Math Processing Error] A(x)A(x)下标的奇偶性把[Math Processing Error] A(x)A(x)分成两半，右边再提一个[Math Processing Error] xx
[Math Processing Error] A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})
A(x)=(a
0
+a
2
x
2
+...+a
n−2
x
n−2
)+(a
1
x+a
3
x
3
+...+a
n−1
x
n−1
)
[Math Processing Error] =(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})
=(a
0
+a
2
x
2
+...+a
n−2
x
n−2
)+x(a
1
+a
3
x
2
+...+a
n−1
x
n−2
)
两边好像非常相似，于是再设两个多项式[Math Processing Error] A_1(x),A_2(x)A
1
(x),A
2
(x)，令
[Math Processing Error] A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{nover 2}-1}
A
1
(x)=a
0
+a
2
x+a
4
x
2
+...+a
n−2
x
2
n
−1
[Math Processing Error] A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{n over 2}-1}
A
2
(x)=a
1
+a
3
x+a
5
x
2
+...+a
n−1
x
2
n
−1
比较明显得出
[Math Processing Error] A(x)=A_1(x^2)+xA_2(x^2)
A(x)=A
1
(x
2
)+xA
2
(x
2
)
再设[Math Processing Error] k<{nover 2}k<
2
n
，把[Math Processing Error] omega^k_nω
n
k
作为[Math Processing Error] xx代入[Math Processing Error] A(x)A(x)（接下来几步变换要多想想单位根的性质）
[Math Processing Error] A(omega^k_n)=A_1((omega^k_n)^2)+omega^k_nA_2((omega^k_n)^2)
A(ω
n
k
)=A
1
((ω
n
k
)
2
)+ω
n
k
A
2
((ω
n
k
)
2
)
[Math Processing Error] =A_1(omega^{2k}_n)+omega^k_nA_2(omega^{2k}_n)=A_1(omega^k_{nover2})+omega^k_nA_2(omega^k_{nover 2})
=A
1
(ω
n
2k
)+ω
n
k
A
2
(ω
n
2k
)=A
1
(ω
2
n
k
)+ω
n
k
A
2
(ω
2
n
k
)
那么对于[Math Processing Error] A(omega^{k+{nover2}}_n)A(ω
n
k+
2
n
)的话，代入[Math Processing Error] omega^{k+{n over 2}}_nω
n
k+
2
n
[Math Processing Error] A(omega^{k+{nover 2}}_n)=A_1(omega^{2k+n}_n)+omega^{k+{nover 2}}_nA_2(omega^{2k+n}_n)
A(ω
n
k+
2
n
)=A
1
(ω
n
2k+n
)+ω
n
k+
2
n
A
2
(ω
n
2k+n
)
[Math Processing Error] =A_1(omega^{2k}_nomega^n_n)-omega^k_nA_2(omega^{2k}_nomega^n_n)
=A
1
(ω
n
2k
ω
n
n
)−ω
n
k
A
2
(ω
n
2k
ω
n
n
)
[Math Processing Error] =A_1(omega^{2k}_n)-omega^k_nA_2(omega^{2k}_n)=A_1(omega^k_{nover2})-omega^k_nA_2(omega^k_{nover2})
=A
1
(ω
n
2k
)−ω
n
k
A
2
(ω
n
2k
)=A
1
(ω
2
n
k
)−ω
n
k
A
2
(ω
2
n
k
)
发现了什么？
[Math Processing Error] A(omega^k_n)A(ω
n
k
)和[Math Processing Error] A(omega^{k+{nover2}}_n)A(ω
n
k+
2
n
)两个多项式后面一坨东西只有符号不同
就是说，如果已知[Math Processing Error] A_1(omega^k_n)A
1
(ω
n
k
)和[Math Processing Error] A_2(omega^k_n)A
2
(ω
n
k
)的值，我们就可以同时知道[Math Processing Error] A(omega^k_n)A(ω
n
k
)和[Math Processing Error] A(omega^{k+{nover2}}_n)A(ω
n
k+
2
n
)的值
现在我们就可以递归分治来搞[Math Processing Error] FFTFFT了
每一次回溯时只扫当前前面一半的序列，即可得出后面一半序列的答案
[Math Processing Error] n==1n==1时只有一个常数项，直接[Math Processing Error] returnreturn
时间复杂度[Math Processing Error] O(nlog_2n)O(nlog
2
n)
IFFT（快速傅里叶逆变换）
想一下，我们不仅要会[Math Processing Error] FFTFFT，还要会IFFT（快速傅里叶逆变换）
我们把两个多项式相乘 （也叫求卷积），做完两遍[Math Processing Error] FFTFFT也知道了积的多项式的点值表示
可我们平时用系数表示的多项式，点值表示没有意义
怎么把点值表示的多项式快速转回系数表示法？
[Math Processing Error] IDFTIDFT暴力[Math Processing Error] O(n^2)O(n
2
)做？其实也可以用[Math Processing Error] FFTFFT用[Math Processing Error] O(nlog_2n)O(nlog
2
n)的时间搞
你有没有想过为什么傅里叶是把[Math Processing Error] omega^k_nω
n
k
作为[Math Processing Error] xx代入而不是别的什么数？
原因是有的但是有我也看不懂
由于我是沙雕所以只用记住一个结论
一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以[Math Processing Error] nn即为原多项式的每一项系数
意思就是说[Math Processing Error] FFTFFT和[Math Processing Error] IFFTIFFT可以一起搞
朴素版FFT板子
#include<complex>
#define cp complex<double>
void fft(cp *a,int n,int inv)//inv是取共轭复数的符号
{
if (n==1)return;
int mid=n/2;
static cp b[MAXN];
fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1];
fo(i,0,n-1)a[i]=b[i];
fft(a,mid,inv),fft(a+mid,mid,inv);//分治
fo(i,0,mid-1)
{
cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv取决是否取共轭复数
b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid];
}
fo(i,0,n-1)a[i]=b[i];
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
两个多项式[Math Processing Error] a,ba,b相乘再转系数多项式[Math Processing Error] cc，通常只用打这么一小段
cp a[MAXN],b[MAXN];
int c[MAXN];
fft(a,n,1),fft(b,n,1);//1系数转点值
fo(i,0,n-1)a[i]*=b[i];
fft(a,n,-1);//-1点值转系数
fo(i,0,n-1)c[i]=(int)(a[i]/n+0.5);//注意精度
1
2
3
4
5
6
很明显，FFT只能处理[Math Processing Error] nn为[Math Processing Error] 22的整数次幂的多项式
所以在[Math Processing Error] FFTFFT前一定要把[Math Processing Error] nn调到[Math Processing Error] 22的次幂
这个板子看着好像很优美，但是
递归常数太大，要考虑优化…
FFTの优化——迭代版FFT
这个图也是盗的
这个很容易发现点什么吧？
每个位置分治后的最终位置为其二进制翻转后得到的位置
这样的话我们可以先把原序列变换好，把每个数放在最终的位置上，再一步一步向上合并
一句话就可以[Math Processing Error] O(n)O(n)预处理第[Math Processing Error] ii位最终的位置[Math Processing Error] rev[i]rev[i]
fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
1
至于蝴蝶变换它死了其实是我不会
真·FFT板子
void fft(cp *a,int inv)
{
int bit=0;
while ((1<<bit)<n)bit++;
fo(i,0,n-1)
{
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次（就是没交换）
}
for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
{
cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根，pi的系数2已经约掉了
for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度，i是合并到了哪一位
{
cp omega(1,0);
for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分，得到右半部分的答案
{
cp x=a[i+j],y=omega*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
}
}
}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
这个板子好像不是那么好背
至少这个板子已经很优美了
FFT后记
本人版权意识薄弱……
本博客部分知识学习于
https://www.cnblogs.com/RabbitHu/p/FFT.html
https://www.cnblogs.com/zwfymqz/p/8244902.html?mType=Group#_label3
https://blog.csdn.net/ggn_2015/article/details/68922404
[Math Processing Error] NTTNTT我来了
---------------------
作者：路人黑的纸巾
来源：CSDN
原文：https://blog.csdn.net/enjoy_pascal/article/details/81478582
版权声明：本文为博主原创文章，转载请附上博文链接！