FFT

可以用 O(

nlog_{2}n

)的時間將一個用係數表示的多項式轉換成它的點值表示的演算法

多項式有兩種表示方法:係數表示法和點值表示法。

對於多項式

f(x)=\Sigma_{i=0}^{n-1}a_{i}x^{i}

係數表示法:用每一項的係數來表示

f(x)

,即

f(x)=\left\{ a_{0},a_{1},a_{2},...,a_{n-1}\right\}

點值表示法

:用平面直角座標系中n個不同的點來表示一個多項式。

證明:

多項式乘法的處理利器-FFT&&NTT

那麼係數矩陣設為B=

多項式乘法的處理利器-FFT&&NTT

矩陣的行列式等於矩陣轉置的行列式。將B轉置即位範德蒙行列式,因此B的行列式|B|為:

\prod_{1\leq i<j\leq n}\left( x_{i} -x_{j} \right)\ne0

所以B滿秩,存在唯一解,所以任意不同的n個點可以唯一確定一個n次多項式。

點表示法乘法:

多項式乘法的處理利器-FFT&&NTT

現在的主要問題是如何快速的將多項式的係數表示法轉為點表示法進行運算,然後快速的轉為係數表示法,DFT可以在

O(nlog_{2}n)

時間內將係數表示法轉為為點表示法,IDFT可以在

O(nlog_{2}n)

時間內將點表示法轉為係數表示法。

DFT

複數:定義

i^{2}

=−1,一個複數 z 可以表示為z=a+bi (a,b∈R ) 其中a稱為實部,b稱為虛部,i稱為虛數單位。

同時可以把複數看成複平面直角座標系上的一個點,如圖所示:

多項式乘法的處理利器-FFT&&NTT

其中x是實數集中的座標軸,y軸是虛數單位i的座標軸。

這個點 (2,3)表示的複數就是 2+3i或者想象它代表的向量為(2,3)。

一個複數z的模定義為它到原點的距離,記為

\left| z \right|=\sqrt{a^{2}+b^{2}}

複數的運算:

多項式乘法的處理利器-FFT&&NTT

多項式乘法的處理利器-FFT&&NTT

進入正題:

為了方便求解n一般取2的整次冪。如果暴力的找n個點帶入多項式求解,複雜度將達到

O(x^{2})

,但如果帶入一些特殊的點將大大運算。

多項式乘法的處理利器-FFT&&NTT

複數域單位圓上所有的點任何次方的結果模長都為1,我們可以不用做全部的次方運算。

考慮將圓八等分:

多項式乘法的處理利器-FFT&&NTT

橙色點即為n=8時要取的點,從(1,0)點開始,逆時針從0號開始標號,標到 7號

記編號為k的點代表的複數值為

\omega_{n}^{k}

,根據模長相乘,極角相加可知

\left(\omega_{n}^{1} \right)=\omega_{n}^{k}

,其中

\omega_{n}^{1}

稱為n次單位根,每個

\omega

都可以求出:

\omega_{n}^{k}=cos\frac{k}{n}2\pi+isin\frac{k}{n}2\pi

用圖表示為:

多項式乘法的處理利器-FFT&&NTT

w_{n}^{0},w_{n}^{1},...,w_{n}^{n-1}

即為要代入的

x_{n}^{0},x_{n}^{1},...,x_{n}^{n-1}

單位根的性質:

1。

w_{n}^{k}=w_{2n}^{2k}

證明:

多項式乘法的處理利器-FFT&&NTT

2。

w_{n}^{k+\frac{n}{2}}=-w_{n}^{k}

證明:

多項式乘法的處理利器-FFT&&NTT

3。

w_{n}^{0}=w_{n}^{n}=1+0i

簡化運算證明:

多項式乘法的處理利器-FFT&&NTT

多項式乘法的處理利器-FFT&&NTT

觀察上述式子可知:如果已知

A_{1}\left(w_{\frac{n}{2}}^{k} \right)

A_{2}\left(w_{\frac{n}{2}}^{k} \right)

可以同時知道

A\left(w_{n}^{k+\frac{n}{2}} \right)

A\left(w_{n}^{k} \right)

。因此可以用遞迴在

O(nlog_{2}n)

時間內完成DFT運算。

程式碼如圖所示:

IDFT

現在快速把點值表示式轉化為係數表示式。

多項式乘法的處理利器-FFT&&NTT

於是我們發現只需要把

w_{n}^{i}

換成

w_{n}^{-i}

,然後DFT,然後將所得的結果除以n,就可以實現IDFT了。

遞迴版FFT程式碼:

void

fft

int

n

complex

<

double

>*

buffer

int

offset

int

step

complex

<

double

>*

epsilon

{

if

n

==

1

return

int

m

=

n

>>

1

fft

m

buffer

offset

step

<<

1

epsilon

);

fft

m

buffer

offset

+

step

step

<<

1

epsilon

);

for

int

k

=

0

k

!=

m

++

k

{

int

pos

=

2

*

step

*

k

temp

k

=

buffer

pos

+

offset

+

epsilon

k

*

step

*

buffer

pos

+

offset

+

step

];

temp

k

+

m

=

buffer

pos

+

offset

-

epsilon

k

*

step

*

buffer

pos

+

offset

+

step

];

}

for

int

i

=

0

i

!=

n

++

i

buffer

i

*

step

+

offset

=

temp

i

];

}

迭代版FFT和蝴蝶操作:

遞迴可以實現FFT,但常數過大,考慮迭代最佳化操作。

多項式乘法的處理利器-FFT&&NTT

接下來考慮如果翻轉二進位制數:

多項式乘法的處理利器-FFT&&NTT

#include

#include

#include

#include

using

namespace

std

const

int

N

=

300010

const

double

PI

=

acos

-

1

);

int

n

m

struct

Complex

{

double

x

y

Complex

operator

+

const

Complex

&

t

const

{

return

{

x

+

t

x

y

+

t

y

};

}

Complex

operator

-

const

Complex

&

t

const

{

return

{

x

-

t

x

y

-

t

y

};

}

Complex

operator

*

const

Complex

&

t

const

{

return

{

x

*

t

x

-

y

*

t

y

x

*

t

y

+

y

*

t

x

};

}

}

a

N

],

b

N

];

int

rev

N

],

bit

tot

void

fft

Complex

a

[],

int

inv

{

for

int

i

=

0

i

<

tot

i

++

if

i

<

rev

i

])

swap

a

i

],

a

rev

i

]]);

for

int

mid

=

1

mid

<

tot

mid

<<=

1

//從下到上迭代

{

auto

w1

=

Complex

({

cos

PI

/

mid

),

inv

*

sin

PI

/

mid

)});

for

int

i

=

0

i

<

tot

i

+=

mid

*

2

{

auto

wk

=

Complex

({

1

0

});

for

int

j

=

0

j

<

mid

j

++

wk

=

wk

*

w1

{

auto

x

=

a

i

+

j

],

y

=

wk

*

a

i

+

j

+

mid

];

a

i

+

j

=

x

+

y

a

i

+

j

+

mid

=

x

-

y

}

}

}

}

int

main

()

{

scanf

“%d%d”

&

n

&

m

);

for

int

i

=

0

i

<=

n

i

++

scanf

“%lf”

&

a

i

]。

x

);

for

int

i

=

0

i

<=

m

i

++

scanf

“%lf”

&

b

i

]。

x

);

while

((

1

<<

bit

<

n

+

m

+

1

bit

++

tot

=

1

<<

bit

for

int

i

=

0

i

<

tot

i

++

rev

i

=

rev

i

>>

1

>>

1

|

((

i

&

1

<<

bit

-

1

));

fft

a

1

),

fft

b

1

);

//DFT

for

int

i

=

0

i

<

tot

i

++

a

i

=

a

i

*

b

i

fft

a

-

1

);

//IDFT

for

int

i

=

0

i

<=

n

+

m

i

++

printf

“%d ”

int

)(

a

i

]。

x

/

tot

+

0。5

));

return

0

}

NTT

多項式乘法的處理利器-FFT&&NTT

多項式乘法的處理利器-FFT&&NTT

多項式乘法的處理利器-FFT&&NTT

NTT程式碼(51nod1028):

#include

#include

#include

#include

#include

using

namespace

std

typedef

long

long

LL

const

LL

mod

=

998244353

const

int

N

=

1

<<

20

LL

Pow

LL

x

LL

y

){

if

y

return

1LL

LL

xx

=

Pow

x

y

/

2

);

xx

=

xx

*

xx

%

mod

if

y

&

1LL

xx

=

xx

*

x

%

mod

return

xx

}

LL

A

B

a

N

],

b

N

],

R

N

],

g

N

],

n

L

char

str

N

];

void

read

(){

scanf

“%s”

str

);

A

=

strlen

str

);

for

int

i

=

0

i

<

A

i

++

a

A

-

i

-

1

=

str

i

-

‘0’

scanf

“%s”

str

);

B

=

strlen

str

);

for

int

i

=

0

i

<

B

i

++

b

B

-

i

-

1

=

str

i

-

‘0’

}

void

NTT

LL

a

N

],

int

n

){

for

int

i

=

0

i

<

n

i

++

if

i

<

R

i

])

swap

a

i

],

a

R

i

]]);

for

int

d

=

1

t

=

n

>>

1

d

<

n

d

<<=

1

t

>>=

1

for

int

i

=

0

i

<

n

i

+=

d

<<

1

))

for

int

j

=

0

j

<

d

j

++

){

LL

tmp

=

g

t

*

j

*

a

i

+

j

+

d

%

mod

a

i

+

j

+

d

=

a

i

+

j

-

tmp

+

mod

%

mod

a

i

+

j

=

a

i

+

j

+

tmp

%

mod

}

}

int

main

(){

read

();

for

n

=

1

L

=

0

n

<=

A

+

B

n

<<=

1

L

++

);

for

int

i

=

0

i

<

n

i

++

R

i

=

R

i

>>

1

>>

1

|

((

i

&

1

<<

L

-

1

));

g

0

=

1

g

1

=

Pow

3

,(

mod

-

1

/

n

);

for

int

i

=

2

i

<

n

i

++

g

i

=

g

i

-

1

*

g

1

%

mod

NTT

a

n

),

NTT

b

n

);

for

int

i

=

0

i

<

n

i

++

a

i

=

a

i

*

b

i

%

mod

g

0

=

1

g

1

=

Pow

g

1

],

mod

-

2

);

for

int

i

=

2

i

<

n

i

++

g

i

=

g

i

-

1

*

g

1

%

mod

NTT

a

n

);

LL

Inv

=

Pow

n

mod

-

2

);

for

int

i

=

0

i

<

n

i

++

a

i

=

a

i

*

Inv

%

mod

for

int

i

=

0

i

<

n

-

1

i

++

a

i

+

1

+=

a

i

/

10

a

i

%=

10

int

d

for

d

=

n

-

1

d

&&!

a

d

];

d

——

);

for

int

i

=

d

i

>=

0

i

——

putchar

a

i

+

‘0’

);

return

0

}

寫在後面:博主本人時間有限(懶),對於一些基本的概念有些直接從其他部落格截圖,如果原作者看到可聯絡我刪除。