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(n\log_2 n)O(nlog
2
n)的時間解決
[Math Processing Error] FFTFFT名字逼格高，也難懂，其他教程寫得讓人看不太懂，於是自己隨便寫一下
建議對復數、三角函數相關知識有所耳聞 （不會也無所謂）
下面難懂的點我會從網上盜
多項式的系數表示法和點值表示法
[Math Processing Error] FFTFFT其實是一個用[Math Processing Error] O(n\log_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,b\in 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},\theta)(
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,\theta_1)*(a_2,\theta_2)=(a_1a_2,\theta_1+\theta_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}(k\in[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{k\over n}2π+i\sin{k\over n} 2π
ω
n
k
=cos
n
k
2π+isin
n
k
2π
或者說也可以這樣解釋這條式子
注意[Math Processing Error] sin^2\theta+cos^2\theta=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{k\over n}2π+isin{k\over n} 2π=cos{2k\over 2n}2π+isin{2k\over 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^{n\over 2}_n=cos{{n\over 2}\over n}2\pi+isin{{n\over 2}\over n}2\pi=cos\pi+isin\pi=-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^{i\pi}+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^{{n\over 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<{n\over 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_{n\over2})+\omega^k_nA_2(\omega^k_{n\over 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+{n\over2}}_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+{n\over 2}}_n)=A_1(\omega^{2k+n}_n)+\omega^{k+{n\over 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}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^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_{n\over2})-\omega^k_nA_2(\omega^k_{n\over2})
=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+{n\over2}}_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+{n\over2}}_n)A(ω
n
k+
2
n
)的值
現在我們就可以遞歸分治來搞[Math Processing Error] FFTFFT了
每一次回溯時只掃當前前面一半的序列，即可得出后面一半序列的答案
[Math Processing Error] n==1n==1時只有一個常數項，直接[Math Processing Error] returnreturn
時間復雜度[Math Processing Error] O(n\log_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(n\log_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
版權聲明：本文為博主原創文章，轉載請附上博文鏈接！