文章目錄
- 前言
- ~~Fast Fast TLE~~
- 一、FFT是什么?
- 二、FFT可以干什么?
- 1.多項式乘法
- 2.大數乘法
- 三、FFT怎么做?
- 1. 系數表示法和點值表示法
- 3.如何巧妙地進行DFT? - FFT的高明之處
- 4. *單位復根* - FFT的高明之處(二)
- 5. Inverse Fast Fourier Transform - 從點值表示法回到系數表示法
- 四、FFT的遞回實作
- 五、大數乘法的實作(FFT版)
- 后話(~~廢話~~)
前言
Fast Fast TLE
原本做大數乘法的時候是想偷懶的, 就百度了下大數 A * B 的代碼, 無意中發現有使用FFT的做法, 于是便開始了學習 (受苦) , 本文用以記錄僅在演算法應用層面上我對FFT的理解和一個手搓的大數乘法模版.
需要一定的復數知識, 但不多.
筆者對FFT的理解較為粗淺, 因此在下文可能大量略過某些部分的詳細證明, 請見諒!
一、FFT是什么?
快速傅里葉變換( Fast Fourier Transform ), 是一種對離散傅里葉變換( Discrete Fourier Transform )的優化, 可令離散傅里葉變換的運算量大大降低, 在取樣點數量較多的情況下可顯著提升效率
二、FFT可以干什么?
1.多項式乘法
實際上, FFT可以大大提升卷積1的計算效率, 而當卷積的兩個函式為多項式的時候, FFT則能在多項式乘法上大顯神威.
2.大數乘法
關于乘法, 樸素的演算法則是模擬我們自小學會的豎式乘法, 不難證明這樣的演算法的復雜度為O(n2), 然而這樣的演算法不僅容易精度溢位, 效率也不夠高, 應用FFT則可以使演算法的復雜度降低至O(nlogn).
觀察一個數字, 如123456, 換種角度思考, 可以將其轉化為以下的式子:
123456
=
1
?
1
0
5
+
2
?
1
0
4
+
3
?
1
0
3
+
4
?
1
0
2
+
5
?
1
0
1
+
6
?
1
0
0
123456 = 1 *10^5+2*10^4+3*10^3+4*10^2+5*10^1+6*10^0
123456=1?105+2?104+3?103+4?102+5?101+6?100再轉化一次, 令x = 10, 則可以得到:
1
?
x
5
+
2
?
x
4
+
3
?
x
3
+
4
?
x
2
+
5
?
x
1
+
6
?
x
0
1 *x^5+2*x^4+3*x^3+4*x^2+5*x^1+6*x^0
1?x5+2?x4+3?x3+4?x2+5?x1+6?x0
從這個角度思考, 大數乘法不正是x = 10的情況下的多項式乘法嗎?
三、FFT怎么做?
PS:強烈建議您通過這個視頻進行學習并在需要的情況下結合下文理解, 個人認為這個視頻講得相當易懂
1. 系數表示法和點值表示法
任取一個n-1階2 多項式, 將其表示為
a
0
+
a
1
x
+
a
2
x
2
+
.
.
.
+
a
n
?
1
x
n
?
1
a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}
a0?+a1?x+a2?x2+...+an?1?xn?1實際上我們已經得到了它的系數表示, 即
[
a
0
a
1
a
2
.
.
.
a
n
?
1
]
\left[ \begin{array}{l} a_0&a_1&a_2&...&a_{n-1} \end{array} \right]
[a0??a1??a2??...?an?1??]
那么什么是它的點值表示法呢?
對一個多項式
f
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
.
.
.
+
a
n
?
1
x
n
?
1
f(x)= a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}
f(x)=a0?+a1?x+a2?x2+...+an?1?xn?1任取不小于n個xi帶入并計算得到f(xi)即為其點值表示法:
{
(
x
0
,
f
(
x
0
)
,
(
x
1
,
f
(
x
1
)
,
…
,
(
x
n
?
1
,
f
(
x
n
?
1
)
}
\left\{ (x_0,f(x_0),(x_1,f(x_1), \dots,(x_{n-1},f(x_{n-1}) \right\}
{(x0?,f(x0?),(x1?,f(x1?),…,(xn?1?,f(xn?1?)}
將其化為矩陣形式, 有
[
f
(
x
1
)
f
(
x
2
)
?
f
(
x
?
1
)
]
=
[
1
x
0
x
0
2
…
x
0
n
?
1
1
x
1
x
1
2
…
x
1
n
?
1
…
…
…
?
…
1
x
n
?
1
x
n
?
1
2
…
x
n
?
1
n
?
1
]
[
a
0
a
1
?
a
n
?
1
]
\left[ \begin{array}{l} f(x_1)\\ f(x_2)\\ \vdots \\ f(x-1) \end{array} \right] = \left[ \begin{array}{l} 1&x_0&x_0^2&\dots&x_{0}^{n-1}\\ 1&x_1&x_1^2&\dots&x_{1}^{n-1}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&x_{n-1}&x_{n-1}^2&\dots&x_{n-1}^{n-1} \end{array} \right] \left[ \begin{array}{l} a_0\\ a_1\\ \vdots \\ a_{n-1} \end{array} \right]
??????f(x1?)f(x2?)?f(x?1)???????=??????11…1?x0?x1?…xn?1??x02?x12?…xn?12??……?…?x0n?1?x1n?1?…xn?1n?1??????????????a0?a1??an?1????????
PS:從 系數表示法 得到 點值表示法 的這一步即為DFT(離散傅里葉變換)
由范德蒙德行列式易證得當x0 ~ xn-1互不相同時矩陣可逆, 即當所取的xi不小于矩陣的階加一時, 點值表示法和系數表示法一一對應.
令
h
(
x
)
=
f
(
x
)
g
(
x
)
h(x)=f(x)g(x)
h(x)=f(x)g(x)
在點值表示法的基礎上, 要計算h(x)的運算式, 只需要取相同的xi得出f(xi)和g(xi)并相乘, 再從(xi ,f(xi)g(xi))的點值表示法轉換成系數表示法即可得到h(x)的值.
那么如何從點值表示法轉換成系數表示法呢? 仔細觀察上面的矩陣運算式不難得出:
{ ( x 0 , f ( x 0 ) , ( x 1 , f ( x 1 ) , … , ( x n ? 1 , f ( x n ? 1 ) } \left\{ (x_0,f(x_0),(x_1,f(x_1), \dots,(x_{n-1},f(x_{n-1}) \right\} {(x0?,f(x0?),(x1?,f(x1?),…,(xn?1?,f(xn?1?)}
將其化為矩陣形式, 有
[
a
0
a
1
?
a
n
?
1
]
=
[
1
x
0
x
0
2
…
x
0
n
?
1
1
x
1
x
1
2
…
x
1
n
?
1
…
…
…
?
…
1
x
n
?
1
x
n
?
1
2
…
x
n
?
1
n
?
1
]
?
1
[
f
(
x
1
)
f
(
x
2
)
?
f
(
x
?
1
)
]
\left[ \begin{array}{l} a_0\\ a_1\\ \vdots \\ a_{n-1} \end{array} \right] = \left[ \begin{array}{l} 1&x_0&x_0^2&\dots&x_{0}^{n-1}\\ 1&x_1&x_1^2&\dots&x_{1}^{n-1}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&x_{n-1}&x_{n-1}^2&\dots&x_{n-1}^{n-1} \end{array} \right]^{-1} \left[ \begin{array}{l} f(x_1)\\ f(x_2)\\ \vdots \\ f(x-1) \end{array} \right]
??????a0?a1??an?1????????=??????11…1?x0?x1?…xn?1??x02?x12?…xn?12??……?…?x0n?1?x1n?1?…xn?1n?1?????????1??????f(x1?)f(x2?)?f(x?1)???????
PS:從 點值表示法 得到 系數表示法 的這一步即為Inverse DFT(逆離散傅里葉變換)
3.如何巧妙地進行DFT? - FFT的高明之處
對上述的DFT程序, 一個通常的想法是在x軸上隨便挑幾個點代入計算, 比如說x0= 0, x1 = 1, …然而這樣做的話又重新掉入了O(n2)的窠臼之中, 我們依舊需要對n個甚至以上的點計算n次, 而要探討FFT的做法, 讓我們先著眼于兩個簡單的多項式
P
1
=
x
2
P_1=x^2
P1?=x2
因為這是一個偶函式, 當我們確定f(xi)時通過偶函式的性質可以立刻確定f(-xi) = f(xi), 這樣我們任取一個點, 立刻就可以確定其相反數的函式值了.
而對于
P
2
=
x
3
P_2= x^3
P2?=x3
和P1不同, 這是個奇函式, 因此f(-xi) = -f(xi).
據此, 我們不妨將一個一般的多項式(設n為偶數)分解為奇函式的部分和偶函式的部分, 再從奇函式部分提取一個x:
P
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
.
.
.
+
a
n
?
1
x
n
?
1
令
P
e
(
x
)
=
a
0
+
a
2
x
2
+
.
.
.
+
a
n
?
2
x
n
?
2
,
P
o
(
x
)
=
a
1
+
a
3
x
2
+
.
.
.
+
a
n
?
1
x
n
?
2
則
P
(
x
i
)
=
P
e
(
x
i
)
+
x
i
P
o
(
x
i
)
,
P
(
?
x
i
)
=
P
e
(
x
i
)
?
x
i
P
o
(
x
i
)
P(x) = a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\\ 令P_e(x) = a_0+ a_2x^2+...+a_{n-2}x^{n-2} , P_o(x) = a_1+a_3x^2+...+a_{n-1}x^{n-2}\\ 則P(x_i)=P_e(x_i)+x_iP_o(x_i), \\P(-x_i)=P_e(x_i)-x_iP_o(x_i)
P(x)=a0?+a1?x+a2?x2+...+an?1?xn?1令Pe?(x)=a0?+a2?x2+...+an?2?xn?2,Po?(x)=a1?+a3?x2+...+an?1?xn?2則P(xi?)=Pe?(xi?)+xi?Po?(xi?),P(?xi?)=Pe?(xi?)?xi?Po?(xi?)
若令 p = x2, 則Pe(x)和Po(x)則簡化為以下式子:
P
e
(
x
)
=
a
0
+
a
2
x
2
+
.
.
.
+
a
n
?
2
x
n
?
2
P
o
(
x
)
=
a
1
+
a
3
x
2
+
.
.
.
+
a
n
?
1
x
n
?
2
?
令
p
=
x
2
P
e
(
p
)
=
a
0
+
a
2
p
+
.
.
.
+
a
n
?
2
p
n
?
2
2
P
o
(
p
)
=
a
1
+
a
3
p
+
.
.
.
+
a
n
?
1
p
n
?
2
2
P_e(x) = a_0+ a_2x^2+...+a_{n-2}x^{n-2}\\ P_o(x) = a_1+a_3x^2+...+a_{n-1}x^{n-2}\\ \Updownarrow 令p=x^2\\ P_e(p) = a_0 + a_2p+...+a_{n-2}p^{\frac{n-2}{2}}\\ P_o(p) = a_1 + a_3p+...+a_{n-1}p^{\frac{n-2}{2}}\\
Pe?(x)=a0?+a2?x2+...+an?2?xn?2Po?(x)=a1?+a3?x2+...+an?1?xn?2?令p=x2Pe?(p)=a0?+a2?p+...+an?2?p2n?2?Po?(p)=a1?+a3?p+...+an?1?p2n?2?
如果此時將Pe和Po視為新的式子, 重復上述操作則可將Pe和Po繼續降階直到Pe和Po的式子的階降至0, 那么此時Pe和Po的值便為常數a0和a1.
以多項式P(x) = a0 + a1x + a2x2 + a3x2為例, DFT程序需要構造4個根x1, x2, x3, x4并帶入求值, 以上述方法降階則有:
第一步:
P
0
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
a
3
x
3
?
令
t
=
x
2
P
e
0
(
t
)
=
a
0
+
a
2
t
P
o
0
(
t
)
=
a
1
+
a
3
t
P
0
(
x
i
)
=
P
e
0
(
x
i
2
)
+
x
i
P
o
0
(
x
i
2
)
,
P
0
(
?
x
i
)
=
P
e
0
(
x
i
2
)
?
x
i
P
o
0
(
x
i
2
)
P_0(x) = a_0+a_1x+a_2x^2+a_3x^3\\ \Updownarrow令t = x^2\\ P_{e0}(t) = a_0+a_2t\\ P_{o0}(t)= a_1 + a_3t\\\\ P_0(x_i)=P_{e0}(x_i^2)+x_iP_{o0}(x_i^2), \\P_0(-x_i)=P_{e0}(x_i^2)-x_iP_{o0}(x_i^2)
P0?(x)=a0?+a1?x+a2?x2+a3?x3?令t=x2Pe0?(t)=a0?+a2?tPo0?(t)=a1?+a3?tP0?(xi?)=Pe0?(xi2?)+xi?Po0?(xi2?),P0?(?xi?)=Pe0?(xi2?)?xi?Po0?(xi2?)
第二步:
令
P
1
=
P
e
0
(
t
)
=
a
0
+
a
2
t
,
令
P
2
=
P
o
0
(
t
)
=
a
1
+
a
3
t
,
則
P
e
1
(
t
)
=
a
0
,
則
P
e
2
(
t
)
=
a
1
,
P
o
1
(
t
)
=
a
2
P
e
2
(
t
)
=
a
3
P
1
(
x
i
2
)
=
P
e
1
(
x
i
2
)
+
x
i
2
P
o
1
(
x
i
2
)
=
a
0
+
x
i
2
?
a
2
,
P
2
(
x
i
2
)
=
P
e
2
(
x
i
2
)
+
x
i
2
P
o
2
(
x
i
2
)
=
a
1
+
x
i
2
?
a
3
,
P
1
(
?
x
i
2
)
=
P
e
1
(
?
x
i
2
)
?
x
i
2
P
o
1
(
x
i
2
)
=
a
0
?
x
i
2
?
a
2
P
2
(
?
x
i
2
)
=
P
e
2
(
?
x
i
2
)
?
x
i
2
P
o
2
(
x
i
2
)
=
a
1
?
x
i
2
?
a
3
\begin{array}{c|c} 令P_1=P_{e0}(t) = a_0+a_2t, & 令P_2 = P_{o0}(t)= a_1 + a_3t,\\ 則P_{e1}(t) = a_0, & 則P_{e2}(t) = a_1,\\ P_{o1}(t) = a_2 & P_{e2}(t) = a_3\\ P_1(x_i^2) = P_{e1}(x_i^2)+x_i^2P_{o1}(x_i^2) = a_0 + x_i^2\cdot a_2, & P_2(x_i^2) = P_{e2}(x_i^2)+x_i^2P_{o2}(x_i^2) = a_1 + x_i^2\cdot a_3,\\ P_1(-x_i^2) = P_{e1}(-x_i^2)-x_i^2P_{o1}(x_i^2) = a_0 - x_i^2\cdot a_2 & P_2(-x_i^2) = P_{e2}(-x_i^2)-x_i^2P_{o2}(x_i^2)=a_1 - x_i^2\cdot a_3\\ \end{array}
令P1?=Pe0?(t)=a0?+a2?t,則Pe1?(t)=a0?,Po1?(t)=a2?P1?(xi2?)=Pe1?(xi2?)+xi2?Po1?(xi2?)=a0?+xi2??a2?,P1?(?xi2?)=Pe1?(?xi2?)?xi2?Po1?(xi2?)=a0??xi2??a2??令P2?=Po0?(t)=a1?+a3?t,則Pe2?(t)=a1?,Pe2?(t)=a3?P2?(xi2?)=Pe2?(xi2?)+xi2?Po2?(xi2?)=a1?+xi2??a3?,P2?(?xi2?)=Pe2?(?xi2?)?xi2?Po2?(xi2?)=a1??xi2??a3??
如此, 對P0(x)的快速傅里葉變換程序已經呼之欲出了, 然而這里卻有一個嚴重的問題亟待解決……
4. 單位復根 - FFT的高明之處(二)
也許你已經注意到了其中的問題, 讓我們回到上述第二步, 該如何取xi的值才能使得x02與x12互為相反數呢, 顯然在實數域中這個問題是無法被解決的.
這就是FFT的第二個高明之處, 使用n次單位復根來解決問題.
讓我們繼續使用上面那個多項式P0(x), 我們需要4個數作為x代入多項式并求值, 從而得出P0(x)的點值運算式, 那么我們不妨使用4次單位復根(1, -1, i, -i)作為 xi :
求
P
0
(
x
i
)
的
值
,
其
中
P
0
(
x
)
=
a
0
+
a
1
x
+
a
2
x
2
+
a
3
x
3
,
x
i
=
[
1
,
?
1
,
i
,
?
i
]
P
0
(
1
)
=
P
e
0
(
1
)
+
1
?
P
o
0
(
1
)
,
P
0
(
?
1
)
=
P
e
0
(
1
)
?
1
?
P
o
0
(
1
)
P
0
(
i
)
=
P
e
0
(
?
1
)
+
i
?
P
o
0
(
?
1
)
,
P
0
(
?
1
)
=
P
e
0
(
?
1
)
?
i
?
P
o
0
(
?
1
)
求P_0(x_i)的值, 其中P_0(x) = a_0+a_1x+a_2x^2+a_3x^3 ,x_i = \left[ 1,-1,i,-i \right]\\ P_0(1)=P_{e0}(1)+1 \cdot P_{o0}(1), \\P_0(-1)=P_{e0}(1)-1 \cdot P_{o0}(1)\\ P_0(i)=P_{e0}(-1)+i \cdot P_{o0}(-1), \\P_0(-1)=P_{e0}(-1)-i \cdot P_{o0}(-1)\\
求P0?(xi?)的值,其中P0?(x)=a0?+a1?x+a2?x2+a3?x3,xi?=[1,?1,i,?i]P0?(1)=Pe0?(1)+1?Po0?(1),P0?(?1)=Pe0?(1)?1?Po0?(1)P0?(i)=Pe0?(?1)+i?Po0?(?1),P0?(?1)=Pe0?(?1)?i?Po0?(?1)
P
e
0
(
1
)
=
P
1
(
1
)
=
P
e
1
(
1
)
+
1
?
P
o
1
(
1
)
=
a
0
+
1
?
a
2
,
P
e
0
(
?
1
)
=
P
1
(
?
1
)
=
P
e
1
(
?
1
)
?
1
?
P
o
1
(
1
)
=
a
0
?
1
?
a
2
P
o
0
(
1
)
=
P
1
(
1
)
=
P
e
2
(
1
)
+
1
?
P
o
2
(
1
)
=
a
1
+
1
?
a
3
,
P
o
0
(
?
1
)
=
P
1
(
?
1
)
=
P
e
2
(
?
1
)
?
1
?
P
o
2
(
1
)
=
a
1
?
1
?
a
3
P_{e0}(1) = P_1(1) = P_{e1}(1)+1\cdot P_{o1}(1) = a_0 + 1\cdot a_2, \\ P_{e0}(-1) = P_1(-1) = P_{e1}(-1)-1\cdot P_{o1}(1) = a_0 - 1\cdot a_2\\ P_{o0}(1) = P_1(1) = P_{e2}(1)+1\cdot P_{o2}(1) = a_1 + 1\cdot a_3, \\ P_{o0}(-1) = P_1(-1) = P_{e2}(-1)-1\cdot P_{o2}(1) = a_1 - 1\cdot a_3
Pe0?(1)=P1?(1)=Pe1?(1)+1?Po1?(1)=a0?+1?a2?,Pe0?(?1)=P1?(?1)=Pe1?(?1)?1?Po1?(1)=a0??1?a2?Po0?(1)=P1?(1)=Pe2?(1)+1?Po2?(1)=a1?+1?a3?,Po0?(?1)=P1?(?1)=Pe2?(?1)?1?Po2?(1)=a1??1?a3?
再將值代入上式, 整個FFT程序便結束了.
當然單位復根的其他(我并不熟知)的性質(如可約引理, 等分引理, 求和引理等), 造就了其十分適合用于快速傅里葉變換中.
簡單來說, 對于n次單位復根(其中n是2的正整數次冪), 他們之間兩兩正負配對(第k個根與第k + n/2個根)且平方后同樣兩兩正負配對, 這令其十分適合用于遞回當中.
當然, 這些性質在n是2的正整數次冪時才生效, 當n不是2的正整數次冪時, 可自行將多項式的次數補足, 此時該項系數為0即可.
當然由于我十分的菜, 具體的證明就此略過了, 各位可以看看我推薦的那個視頻或者自行尋找資料進行證明(
5. Inverse Fast Fourier Transform - 從點值表示法回到系數表示法
這一步驟之簡單你絕對無法想象, 觀察上文的IDFT的實作:
[
a
0
a
1
?
a
n
?
1
]
=
[
1
x
0
x
0
2
…
x
0
n
?
1
1
x
1
x
1
2
…
x
1
n
?
1
…
…
…
?
…
1
x
n
?
1
x
n
?
1
2
…
x
n
?
1
n
?
1
]
?
1
[
f
(
x
1
)
f
(
x
2
)
?
f
(
x
?
1
)
]
\left[ \begin{array}{l} a_0\\ a_1\\ \vdots \\ a_{n-1} \end{array} \right] = \left[ \begin{array}{l} 1&x_0&x_0^2&\dots&x_{0}^{n-1}\\ 1&x_1&x_1^2&\dots&x_{1}^{n-1}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&x_{n-1}&x_{n-1}^2&\dots&x_{n-1}^{n-1} \end{array} \right]^{-1} \left[ \begin{array}{l} f(x_1)\\ f(x_2)\\ \vdots \\ f(x-1) \end{array} \right]
??????a0?a1??an?1????????=??????11…1?x0?x1?…xn?1??x02?x12?…xn?12??……?…?x0n?1?x1n?1?…xn?1n?1?????????1??????f(x1?)f(x2?)?f(x?1)???????
讓我們對第二個矩陣進行求逆:
[
1
x
0
x
0
2
…
x
0
n
?
1
1
x
1
x
1
2
…
x
1
n
?
1
…
…
…
?
…
1
x
n
?
1
x
n
?
1
2
…
x
n
?
1
n
?
1
]
?
1
=
[
1
1
1
…
1
1
ω
ω
2
…
ω
n
?
1
…
…
…
?
…
1
ω
n
?
1
ω
2
(
n
?
1
)
…
ω
(
n
?
1
)
(
n
?
1
)
]
?
1
=
1
n
[
1
1
1
…
1
1
ω
?
1
ω
?
2
…
ω
?
(
n
?
1
)
…
…
…
?
…
1
ω
?
(
n
?
1
)
ω
?
2
(
n
?
1
)
…
ω
?
(
n
?
1
)
(
n
?
1
)
]
\left[ \begin{array}{l} 1&x_0&x_0^2&\dots&x_{0}^{n-1}\\ 1&x_1&x_1^2&\dots&x_{1}^{n-1}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&x_{n-1}&x_{n-1}^2&\dots&x_{n-1}^{n-1} \end{array} \right]^{-1} = \left[ \begin{array}{l} 1&1 &1&\dots&1\\ 1&\omega &\omega ^2 &\dots&\omega ^{n-1}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&\omega ^{n-1}&\omega ^{2(n-1)}&\dots&\omega^{(n-1)(n-1)} \end{array} \right]^{-1}=\frac{1}{n}\left[ \begin{array}{l} 1&1 &1&\dots&1\\ 1&\omega^{-1} &\omega ^{-2} &\dots&\omega ^{-(n-1)}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&\omega ^{-(n-1)}&\omega ^{-2(n-1)}&\dots&\omega^{-(n-1)(n-1)} \end{array} \right]
??????11…1?x0?x1?…xn?1??x02?x12?…xn?12??……?…?x0n?1?x1n?1?…xn?1n?1?????????1=??????11…1?1ω…ωn?1?1ω2…ω2(n?1)?……?…?1ωn?1…ω(n?1)(n?1)????????1=n1???????11…1?1ω?1…ω?(n?1)?1ω?2…ω?2(n?1)?……?…?1ω?(n?1)…ω?(n?1)(n?1)???????其中
x
k
=
ω
k
,
ω
=
e
2
π
i
n
x_k = \omega_k, \omega = e^{\frac{2\pi i}{n}}
xk?=ωk?,ω=en2πi?
也就是說:
[
p
0
p
1
?
p
n
?
1
]
=
1
n
[
1
1
1
…
1
1
ω
ω
2
…
ω
n
?
1
…
…
…
?
…
1
ω
n
?
1
ω
2
(
n
?
1
)
…
ω
(
n
?
1
)
(
n
?
1
)
]
[
P
(
x
1
)
P
(
x
2
)
?
P
(
x
?
1
)
]
\left[ \begin{array}{l} p_0\\ p_1\\ \vdots \\ p_{n-1} \end{array} \right] = \frac{1}{n}\left[ \begin{array}{l} 1&1 &1&\dots&1\\ 1&\omega &\omega ^2 &\dots&\omega ^{n-1}\\ \dots&\dots&\dots&\vdots&\dots\\ 1&\omega ^{n-1}&\omega ^{2(n-1)}&\dots&\omega^{(n-1)(n-1)} \end{array} \right] \left[ \begin{array}{l} P(x_1)\\ P(x_2)\\ \vdots \\ P(x-1) \end{array} \right]
??????p0?p1??pn?1????????=n1???????11…1?1ω…ωn?1?1ω2…ω2(n?1)?……?…?1ωn?1…ω(n?1)(n?1)?????????????P(x1?)P(x2?)?P(x?1)???????
不難發現, 要實作IFFT, 只需在FFT的基礎上將ω取倒數, 并在最終結果乘以n分之一即可.
PS:我覺得自己寫的FFT解釋很爛, 在此建議您如果 (肯定) 有看不懂的部分結合剛才推薦的視頻學習
四、FFT的遞回實作
據此, 我們終于可以寫出FFT的代碼了, 我手寫的遞回實作如下:
#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
const double Pi = acos(-1.0); //將常數Pi的值設為arccos(-1)
vector<complex<double> > FFT(vector<complex<double> > &a, int flag) {
//第一個引數為一個多項式的系數, 以次數從小到大的順序, 向量中每一項的實部為該項系數
//flag == 1->FFT, flag == -1->IFFT
int n = a.size();
if(n == 1) //如果當前多項式僅有常數項時直接回傳多項式的值
return a;
vector<complex<double> > Pe, Po; //即原文中的Pe與Po的系數表示法
complex<double> omega(cos(2.0 * Pi / n), sin(flag * 2.0 * Pi / n)), cur(1, 0); //omega為第一個n次復根, cur為第零個n次復根即1
for (unsigned int i = 0; i < n; i++) {
if(i & 1)
Po.push_back(a[i]);
else
Pe.push_back(a[i]);
}
vector<complex<double> > ye = FFT(Pe, flag), yo = FFT(Po, flag), y(n); //遞回求得ye=Pe(xi), yo=Po(xi)
for (int i = 0; i < n / 2; i++) {
y[i] = ye[i] + cur * yo[i]; //求出P(xi)
y[i + n / 2] = ye[i] - cur * yo[i]; //由單位復根的性質可知第k個根與第k + n/2個根互為相反數
cur *= omega; //cur * omega得到下一個復根
}
return y; //回傳最終的系數
}
當然需要注意的是, 如果當前進行的是IFFT時, 需要將最終結果除以n, 其中n為不小于原多項式階+1的最小2的整數次冪.
因此你需要在上述代碼的基礎上增加一個代碼, 向原系數向量補0直到向量中元素的個數為二的整數冪為止.
五、大數乘法的實作(FFT版)
了解了FFT的實作方法, 大數乘法的實作可以說是再簡單不過了(畢竟前面的鬼東西都看懂了), 這里便直接貼上源代碼了,當然有一個非常需要注意的地方, 最侄訓取的多項式系數是有可能大于等于10的, 此時需要做進位處理!
#include <bits/stdc++.h>
#define sync ios::sync_with_stdio(false)
using namespace std;
typedef unsigned long long ull;
const double Pi = acos(-1.0); //將常數Pi的值設為arccos(-1)
vector<complex<double> > FFT(vector<complex<double> > &a, int flag) {
//第一個引數為一個多項式的系數, 以次數從小到大的順序, 向量中每一項的實部為該項系數
//flag == 1->FFT, flag == -1->IFFT
int n = a.size();
if(n == 1) //如果當前多項式僅有常數項時直接回傳多項式的值
return a;
vector<complex<double> > Pe, Po; //即原文中的Pe與Po的系數表示法
complex<double> omega(cos(2.0 * Pi / n), sin(flag * 2.0 * Pi / n)), cur(1, 0); //omega為第一個n次復根, cur為第零個n次復根即1
for (unsigned int i = 0; i < n; i++) {
if(i & 1)
Po.push_back(a[i]);
else
Pe.push_back(a[i]);
}
vector<complex<double> > ye = FFT(Pe, flag), yo = FFT(Po, flag), y(n); //遞回求得ye=Pe(xi), yo=Po(xi)
for (int i = 0; i < n / 2; i++) {
y[i] = ye[i] + cur * yo[i]; //求出P(xi)
y[i + n / 2] = ye[i] - cur * yo[i]; //由單位復根的性質可知第k個根與第k + n/2個根互為相反數
cur *= omega; //cur * omega得到下一個復根
}
return y; //回傳最終的系數
}
vector<complex<double> > read() { //獲取系數, 注意需要從右向左獲取(從0次冪的系數開始)
string num;
cin >> num;
vector<complex<double> > result;
for (int i = num.size() - 1; i >= 0; i--) {
complex<double> tmp(num[i] - '0', 0);
result.push_back(tmp);
}
return result;
}
void solve(vector<complex<double> > &a, vector<complex<double> > &b) { //多項式系數運算式的修飾
complex<double> tmp(0, 0);
int sum = a.size() + b.size();
while (a.size() < sum)
a.push_back(tmp);
while (b.size() < sum)
b.push_back(tmp);
//如果兩式的階不同, 先補齊
int temp = 1;
while (temp < a.size())
temp <<= 1;
//獲取不小于n的最小2的整數次冪
while (a.size() < temp) {
a.push_back(tmp);
b.push_back(tmp);
}
//補齊
}
int main() {
sync;
vector<complex<double> > num1 = read(), num2 = read(), tmp1, tmp2, mid, ans;
solve(num1, num2);
tmp1 = FFT(num1, 1), tmp2 = FFT(num2, 1);
num1.clear();
num2.clear();
for (int i = 0; i < tmp1.size(); i++)
mid.push_back(tmp1[i] * tmp2[i]);
tmp1.clear();
tmp2.clear();
ans = FFT(mid, -1);
bool Ans = false;
int add = 0;
string final;
//進位處理
for (int i = 0; i < ans.size(); i++) {
int op = round(round(ans[i].real()) / ans.size()) + add;
add = 0;
if(op >= 10)
add = op / 10;
final += op % 10 + '0';
}
if(add > 0)
final += add % 10 + '0';
for (int i = final.size() - 1; i >= 0; i--) {
if(final[i] != '0')
Ans = true;
else if(!Ans)
continue;
cout << final[i];
}
cout << '\n';
}
后話(廢話)
- 歷時一個寒假! 終于把這篇萬字長文(x)寫完了, 多謝某csome同學的鼓勵~~(催更)~~ !
- 理解FFT足足花了2天時間, 好難……
- 如果博客有錯漏, 敬請指出并不吝賜教, 感謝!
卷積的定義請參見維基對卷積的定義. ??
多項式的次數最大的項的次數即為多項式的階. ??
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/264155.html
標籤：其他