前言

快速傅立葉變換 (Fast Fourier Transform)

，即利用計算機計算離散傅立葉變換（DFT）的高效、快速計算方法的統稱，簡稱FFT，於1965年由J。W。庫利和T。W。圖基提出。

FFT用於用來

加速多項式乘法

，可以以

\mathrm O(n\log n)

的時間複雜度解決樸素高精度

\mathrm O(n^2)

複雜度能夠解決的多項式乘法。

複數

係數與點值表示法

對於任意的

n

項

n-1

次整係數多項式

f(x)

，顯而易見我們可以將其表示為

a_0x^{n-1}+a_1x^{n-2}+...+a_{n-1}x^0

的形式。

我們同樣可以用每一項的係數表示

f(x)

，即：

f(x)=\left\{ a_0,a_1,a_2,...,a_{n-1} \right\}\quad\\

這就是

係數表示法

。

當我們將多項式

f(x)

看作一個

函式

時，如果我們將

n

個不同的數

x_0,x_1,...,x_{n-1}

帶入，就能得到

n

個不同的點。

這

n

個不同的點同樣可以看作是

n

個不同的一元

n-1

次方程聯立，這

n

個點或者方程也能反過來

唯一確定

f(x)

。

所以

f(x)

也可以表示為

f(x)=\left\{ (x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1})) \right\}\quad\\

這就是

點值表示法

。

我們假設兩個多項式

f(x),g(x)

分別為：

f(x)=\left\{ (x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n-1},f(x_{n-1})) \right\}\quad\\

g(x)=\left\{ (x_0,g(x_0)),(x_1,g(x_1)),...,(x_{n-1},g(x_{n-1})) \right\}\quad \\

其乘積為：

h(x)=\left\{ (x_0,f(x_0)\cdot g(x_0)),(x_1,f(x_1)\cdot g(x_1)),...,(x_{n-1},f(x_{n-1})\cdot g(x_{n-1})) \right\}\quad \\

顯然點值表示法完成多項式乘法的時間複雜度僅有

\mathrm O(n)

，所以剩下的問題就是如何將係數表示法轉為點值表示法了。

單位根

以單位圓點為起點，單位圓的

n

等分點為終點，在單位圓上可以得到

n

個複數，設幅角為正且最小的複數為

\omega_n

，稱為

n

次單位根

。

由於複數乘法

模數相乘

，

輻角相加

的性質，編號為

k

的點

\omega_n^k=(\omega_n)^k

。

由

尤拉公式

可得

\omega_{n}^{k}=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}\quad\\

顯然

\omega_n^0=\omega_n^n

。

單位根同樣具有以下性質：

\begin{align}\omega_{rn}^{rk}&=\cos \frac{2rk\pi}{rn}+i\sin \frac{2rk\pi}{rn}=\omega_n^k \quad r \in \mathbb{Z}^+\\\omega _{n}^{k+\frac{n}{2}}&=\omega_n^k \cdot (\cos\frac{n}{2}*\frac{2\pi}{n}+i\sin\frac{n}{2}*\frac{2\pi}{n} )\\&=\omega_n^k \cdot (\cos \pi +i\sin \pi)\\&=-\omega _{n}^{k}\\ \color{red}{\overline{\omega_n^k}}&\color{red}{=\cos\frac{2k\pi}n-i\sin\frac{2k\pi}n}\\ &\color{red}{=\cos(2\pi-\frac{2k\pi}n)+i\sin(2\pi-\frac{2k\pi}n)}\\ &\color{red}{=\omega_n^{n-k}}\end{align}\\

上述性質從幾何的角度看就很容易理解。

快速傅立葉變換

我們假設多項式

f(x)=a_0+a_1*x+ \dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}

，那麼我們可以按

a

下標的奇偶性將

f(x)

分為兩半。

也就是

\begin{align}f(x)&=(a_0+a_2{x^2}+a_4{x^4}+\dots+a_{n-2}x^{n-2})\\&+(a_1x+a_3{x^3}+a_5{x^5}+ \dots+a_{n-1}x^{n-1})\end{align}\quad\\

我們令

f_1(x)=a_0+a_2{x^1}+a_4{x^2}+\dots+a_{n-2}x^{\frac{n}{2}-1}\quad\\

f_2(x)=x(a_1+a_3{x}+a_5{x^2}+ \dots+a_{n-1}x^{\frac{n}{2}-1})\quad\\

顯然

f(x)=f_1(x^2)+xf_2(x^2)

。

帶入

\omega_n^k(k<\frac{n}{2})

可得

\begin{align}f(\omega_n^k)&=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})\\&=f_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kf_2(\omega_{\frac{n}{2}}^{k}) \end{align}\quad\\

帶入

\omega_n^{k+\frac{n}{2}}(k<\frac{n}{2})

可得

\begin{align}f(\omega_n^{k+\frac{n}{2}})&=f_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}f_2(\omega_n^{2k+n})\\&=f_1(\omega_n^{2k}*\omega_n^n)-\omega_n^kf_2(\omega_n^{2k}*\omega_n^n)\\&=f_1(\omega_n^{2k})-\omega_n^kf_2(\omega_n^{2k}) \end{align}\quad\\

我們發現上面的兩個式子只有一個常數項不同，所以當我們求出第一個式子的值後可以以

\mathrm O(1)

的時間複雜度求出第二個式子的值。

所以我們就成功將原來的問題縮小到了之前的一半，而縮小後的問題還能繼續縮小。

所以顯然FFT的時間複雜度即為

\mathrm O(nlogn)

。

快速傅立葉逆變換

我們把兩個多項式相乘（或者說求

卷積

），我們跑完FFT並相乘後我們就得到了積的多項式的

點值表示

，現在我們需要將點值表示轉為

係數表示

，這個轉換的過程我們稱為

離散傅立葉逆變換

（

IDFT

）。

如果我們用矩陣將DFT的過程封裝，那麼DFT就相當於：

\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1) ^ {n-1}\\1 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1}) ^ {n-1}\end{bmatrix} \begin{bmatrix}a_0 \\ a _ 1 \\ a _ 2 \\ \vdots \\ a _ {n-1}\\ \end{bmatrix} = \begin{bmatrix}y_0 \\ y _ 1 \\ y _ 2 \\ \vdots \\ y _ {n-1}\\ \end{bmatrix}\quad\\

其中

A=\begin{bmatrix}a_0 \\ a _ 1 \\ a _ 2 \\ \vdots \\ a _ {n-1}\\ \end{bmatrix}\quad\\

意義為多項式的係數表示法，

B=\begin{bmatrix}y_0 \\ y _ 1 \\ y _ 2 \\ \vdots \\ y _ {n-1}\\ \end{bmatrix}\quad\\

意義為多項式的點值表示法，

y_i

即為

f(\omega_n^i)

。

注意到

x^n-1=(x-1)(x^{n-1}+x^{n-2}+\cdots+1)\\

於是在

x=1

時

x^{n-1}+x^{n-2}+\cdots+1=n\\

否則

x^{n-1}+x^{n-2}+\cdots+1=0\\

其中

x

為

n

次單位根。

受

n

次單位根這一特性的啟發，注意到

x^{n-1}+x^{n-2}+\cdots+1

在其形式上與矩陣乘法的關係，可以發現

\begin{align}\begin{bmatrix}1&(\omega_n^{n-i})^1&({\omega_n^{n-i}})^2&\cdots&(\omega_{n}^{n-i})^{n-1}\end{bmatrix}\begin{bmatrix}1\\\ (\omega_n^{j})^1\\(\omega_n^{j})^2\\\vdots\\(\omega_n^{j})^{n-1}\end{bmatrix}&=\sum_{k=0}^{n-1}(\omega_n^{n-i})^k(\omega_n^j)^k\\ &=\sum_{k=0}^{n-1}(\omega_n^{n-i+j})^k\\ &=\begin{cases}n&i=j\\0&i\ne j\end{cases}\end{align}\\

這表明

\mathbf C=\frac1n\begin{bmatrix}1&1&1&\cdots&1\\ 1&(\omega_n^{n-1})^1&({\omega_n^{n-1}})^2&\cdots&(\omega_{n}^{n-1})^{n-1}\\ 1&(\omega_n^{n-2})^1&({\omega_n^{n-2}})^2&\cdots&(\omega_{n}^{n-2})^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&(\omega_n^1)^1&{(\omega_n^1)}^2&\cdots&(\omega_{n}^1)^{n-1}\end{bmatrix}\\

即為

\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1) ^ {n-1}\\1 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1}) ^ {n-1}\end{bmatrix}\\

的逆矩陣。

由上文標紅部分

\begin{align} \mathbf C&=\frac1n\begin{bmatrix}1&1&1&\cdots&1\\ 1&(\omega_n^{n-1})^1&({\omega_n^{n-1}})^2&\cdots&(\omega_{n}^{n-1})^{n-1}\\ 1&(\omega_n^{n-2})^1&({\omega_n^{n-2}})^2&\cdots&(\omega_{n}^{n-2})^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&(\omega_n^1)^1&{(\omega_n^1)}^2&\cdots&(\omega_{n}^1)^{n-1}\end{bmatrix}\\ &=\frac1n\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\overline{\omega_n^1})^1 & (\overline{\omega_n^1})^2 & \cdots & (\overline{\omega_n^1}) ^ {n-1}\\1 & (\overline{\omega_n^2})^1 & (\overline{\omega_n^2})^2 & \cdots & (\overline{\omega_n^2}) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & \overline{(\omega_n^{n-1}})^1 & (\overline{\omega_n^{n-1}})^2 & \cdots & (\overline{\omega_n^{n-1}}) ^ {n-1}\end{bmatrix}\\\end{align}\\

所以

將一個多項式在分治的過程中乘上的單位根變為其共軛複數，分治完的每一項除以

n

即可得到原多項式的每一項係數

。

具體實現（遞迴版）

//LuoguP3803

#include

const

int

NR

=

1

<<

22

；

const

double

eps

=

1e-6

，

pi

=

acos

（

-

1。0

）；

using

namespace

std

；

complex

<

double

>

a

［

NR

］，

b

［

NR

］；

//complex為C++自帶虛數

int

n

，

m

；

void

FFT

（

complex

<

double

>

*

a

，

int

n

，

int

inv

）

//inv為虛部符號，inv為1時FFT，inv為-1時IFFT

{

if

（

n

==

1

）

//如果需要轉換的只有一項就直接返回

return

；

int

mid

=

n

/

2

；

complex

<

double

>

A1

［

mid

+

1

］，

A2

［

mid

+

1

］；

for

（

int

i

=

0

；

i

<=

n

；

i

+=

2

）

//拆分多項式

{

A1

［

i

/

2

］

=

a

［

i

］，

A2

［

i

/

2

］

=

a

［

i

|

1

］；

}

FFT

（

A1

，

mid

，

inv

）；

//遞迴分治

FFT

（

A2

，

mid

，

inv

）；

complex

<

double

>

w0

（

1

，

0

），

wn

（

cos

（

2

*

pi

/

n

），

inv

*

sin

（

2

*

pi

/

n

））；

//單位根

for

（

int

i

=

0

；

i

<

mid

；

++

i

，

w0

*=

wn

）

//合併多項式

{

a

［

i

］

=

A1

［

i

］

+

w0

*

A2

［

i

］；

a

［

i

+

n

/

2

］

=

A1

［

i

］

-

w0

*

A2

［

i

］；

}

}

int

main

（）

{

scanf

（

“%d%d”

，

&

n

，

&

m

）；

for

（

int

i

=

0

；

i

<=

n

；

++

i

）

//輸入第一個多項式

{

double

x

；

scanf

（

“%lf”

，

&

x

）；

a

［

i

］。

real

（

x

）；

//complex型別變數。real（x）意味著將實數部賦為x，real（）返回實數部值

}

for

（

int

i

=

0

；

i

<=

m

；

++

i

）

//輸入第二個多項式

{

double

x

；

scanf

（

“%lf”

，

&

x

）；

b

［

i

］。

real

（

x

）；

}

int

len

=

1

<<

max

（（

int

）

ceil

（

log2

（

n

+

m

）），

1

）；

//由於FFT需要項數為2的整數次方倍，所以多項式次數len為第一個大於n+m的二的正整數次方

FFT

（

a

，

len

，

1

）；

//係數表達轉點值表達

FFT

（

b

，

len

，

1

）；

for

（

int

i

=

0

；

i

<=

len

；

++

i

）

a

［

i

］

=

a

［

i

］

*

b

［

i

］；

//O（n）乘法

FFT

（

a

，

len

，

-

1

）；

//點值表達轉系數表達

for

（

int

i

=

0

；

i

<=

n

+

m

；

++

i

）

//輸出

printf

（

“%。0f ”

，

a

［

i

］。

real

（）

/

len

+

eps

）；

//記得要除n，eps為解決掉精度問題

return

0

；

}

FFT的最佳化

快速數論變換（NTT）

關於某沒有精度損失還可以取模的比FFT快的小常數多項式乘法。jpg

例題