• 微信公众号:美女很有趣。 工作之余,放松一下,关注即送10G+美女照片!

[数学]——[多项式]——-快速傅里叶变换(FFT)

互联网 diligentman 2周前 (04-07) 7次浏览

快速傅里叶变换简介

F

a

s

t

 

F

o

u

r

i

e

r

 

T

r

a

n

s

f

o

r

m

Fast Fourier Transform

Fast Fourier Transform ,简称

F

F

T

FFT

FFT。这东西似乎在数字通信领域有大用处:(EDA表示很感动)

计算量小的显著的优点,使得FFT在信号处理技术领域获得了广泛应用,结合高速硬件就能实现对信号的实时处理。例如,对语音信号的分析和合成,对通信系统中实现全数字化的时分制与频分制(TDM/FDM)的复用转换,在频域对信号滤波以及相关分析,通过对雷达、声纳、振动信号的频谱分析以提高对目标的搜索和跟踪的分辨率等等,都要用到FFT。可以说FFT的出现,对数字信号处理学科的发展起了重要的作用。——百度百科(快速傅里叶变换词条)

但是本博客只介绍其在算法竞赛中的应用——快速计算两个多项式的乘积。

假设现在有两个多项式:

A

(

x

)

=

x

2

+

3

x

+

2

A(x)=x^2+3x+2

A(x)=x2+3x+2

B

(

x

)

=

2

x

2

+

1

B(x)=2x^2+1

B(x)=2x2+1,要求

C

(

x

)

=

A

(

x

)

×

B

(

x

)

C(x)=A(x) times B(x)

C(x)=A(x)×B(x)。朴素做法是按照乘法分配律,将其中一个多项式的每一项分别与另一个多项式的每一项相乘,如果这两个多项式最高次数为

N

N

N,这个算法的时间复杂度将会是

O

(

n

2

)

O(n^2)

O(n2)。然而我们可以用

F

F

T

FFT

FFT

O

(

n

l

o

g

n

)

O(nlogn)

O(nlogn)的时间内计算出

C

(

x

)

C(x)

C(x)

F

F

T

FFT

FFT算法的核心,是四个十分精妙的想法,发明这个算法的人应该是天选者。

下面将介绍这四个巧妙的想法,以及

F

F

T

FFT

FFT的代码和实现细节。

1. 多项式的另一种表示

一般怎么用计算机存储多项式?

开一个数组

A

[

i

]

A[i]

A[i],表示多项式

x

i

x^i

xi 项的系数,例如

F

(

x

)

=

4

x

4

+

6

x

3

+

5

x

2

+

3

x

+

2

F(x)=4x^4+6x^3+5x^2+3x+2

F(x)=4x4+6x3+5x2+3x+2,那么

A

[

]

=

{

4

6

5

3

2

}

A[]={4,6,5,3,2}

A[]={46532}。这完全没有问题,形如

F

(

x

)

=

a

i

x

i

F(x)=sum a_ix^i

F(x)=aixi 的多项式表示方法,叫做“系数表示法”。

然而现在为了更快计算多项式乘积,我们使用多项式的另一种表示法,要理解这种表示方法,可以从最简单的多项式开始。设

A

(

x

)

=

a

0

+

a

1

x

A(x)=a_0+a_1x

A(x)=a0+a1x,一次函数就是一种简单的多项式,它的最高幂次为1。

我们知道

A

(

x

)

=

a

0

+

a

1

x

A(x)=a_0+a_1x

A(x)=a0+a1x 不但是一个多项式,还是一条直线。由于两点能确定一条直线,我们在平面中随便取两点

(

x

1

y

1

)

(x_1,y_1)

(x1y1)

(

x

2

y

2

)

(x_2,y_2)

(x2y2),那么系数

a

0

a_0

a0

a

1

a_1

a1 也就随之确定了。而我们现在在讲的“多项式的另一种表示方法”,就是用平面上的点来确定多项式。

高次多项式的图像是一条曲线。显然,无论你在平面上取多少点,也无法唯一确定一条曲线,但却可以确定一个多项式。因为多项式“并不是普通的(或者说任意的)曲线”。这里有一个结论:平面上的

N

+

1

N+1

N+1 个点,能唯一确定一个

N

N

N 次多项式。例如:你在平面上随机取四个点,那么有且只有一个三次函数同时穿过这四个点。这个结论的正确性需要用线性代数的知识来证明,详细写可能需要一篇新的博客,所以我只能简单说明一下正确性。其实是不会

F

(

x

)

=

i

=

0

i

=

n

a

i

x

i

F(x)=sum_{i=0}^{i=n}a_ix^i

F(x)=i=0i=naixi,我们选取的点为

{

x

0

x

1

x

2

.

.

.

x

n

}

{x_0,x_1,x_2,…,x_n}

{x0x1x2...xn}。那么:

[

F

(

x

0

)

F

(

x

1

)

F

(

x

n

)

]

=

[

1

x

0

x

0

2

x

0

n

1

x

1

x

1

2

x

1

n

1

x

n

x

n

2

x

n

n

]

[

a

0

a

1

a

n

]

left[ begin{matrix} F(x_0) \ F(x_1) \ vdots \F(x_n) end{matrix} right]= left[ begin{matrix} 1 & x_0 &x^2_0 & cdots & x_0^n \ 1 & x_1 &x^2_1 & cdots & x_1^n \ vdots & vdots & vdots & ddots & vdots \ 1 & x_n &x^2_n & cdots & x_n^n \ end{matrix} right] left[ begin{matrix} a_0 \ a_1 \ vdots \ a_n end{matrix} right]

F(x0)F(x1)F(xn)=111x0x1xnx02x12xn2x0nx1nxnna0a1an
易证中间矩阵的行列式不为0,所以方程组只有一组解。好了证毕。

这种表示多项式的方法叫做“点值表示法”。假设对于两个多项式

f

(

x

)

f(x)

f(x)

g

(

x

)

g(x)

g(x),我们选取相同的

x

x

x 序列

{

x

0

,

x

1

.

.

.

.

,

x

n

}

{ x_0, x_1….,x_n}

{x0,x1....,xn},那么两个多项式可以表示为:

f

(

x

)

=

(

 

(

x

0

,

f

(

x

0

)

)

,

(

x

1

,

f

(

x

1

)

)

,

.

.

.

,

(

x

n

,

f

(

x

n

)

)

 

)

g

(

x

)

=

(

 

(

x

0

,

g

(

x

0

)

)

,

(

x

1

,

g

(

x

1

)

)

,

.

.

.

,

(

x

n

,

g

(

x

n

)

)

 

)

f(x)=( (x_0,f(x_0)),(x_1,f(x_1)),…,(x_n,f(x_n)) ) \ g(x)=( (x_0,g(x_0)),(x_1,g(x_1)),…,(x_n,g(x_n)) )

f(x)=( (x0,f(x0)),(x1,f(x1)),...,(xn,f(xn)) )g(x)=( (x0,g(x0)),(x1,g(x1)),...,(xn,g(xn)) )

F

(

x

)

=

f

(

x

)

×

g

(

x

)

F(x)=f(x) times g(x)

F(x)=f(x)×g(x),则

F

(

x

)

=

{

(

x

0

,

f

(

x

0

)

g

(

x

0

)

)

,

(

x

1

,

f

(

x

1

)

g

(

x

1

)

)

,

.

.

.

,

(

x

n

,

f

(

x

n

)

g

(

x

n

)

)

}

F(x)={(x_0,f(x_0)g(x_0)),(x_1,f(x_1)g(x_1)),…,(x_n,f(x_n)g(x_n))}

F(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),...,(xn,f(xn)g(xn))}

容易发现,如果已知两个多项式的点值表示,求两个多项式乘积多项式的点值表示,可以在线性时间内处理完。如果能找到一种方法,能快速将一个多项式在系数表示点值表示之间转化,就能实现快速多项式乘法。

2.系数到点值的转换(DFT)

既然任意

N

+

1

N+1

N+1 个点就可以确定一个

N

N

N 阶多项式,那么如果我直接暴力随机选

N

+

1

N+1

N+1 个横坐标,带入多项式求值会怎样呢?容易想到这样做的时间复杂度是

O

(

n

2

)

O(n^2)

O(n2)(找

N

+

1

N+1

N+1个点

O

(

N

)

O(N)

O(N),每个点求值

O

(

N

)

O(N)

O(N)),这样就又回到了原点,没有达到加速的目的,所以需要找另外的方法。

依旧考虑最简单的情况:现在,我想在二次函数

f

(

x

)

=

x

2

f(x)=x^2

f(x)=x2上,任意取 8 个点,并求出这8个点的函数值,有没有什么比较快的方法?容易想到,其实我只需要任意取四个

x

x

x 就行了。因为如果我知道了

x

=

x

0

x=x_0

x=x0 的函数值,我立刻就知道了

x

=

x

0

x=-x_0

x=x0 的函数值,因为

f

(

x

)

=

x

2

f(x)=x^2

f(x)=x2 是一个偶函数。奇函数同理,只不过我需要在函数值上再加一个负号。总之:我们可以利用奇偶性来减少选取点的数目。

下面的问题就是,如何将一个一般的多项式,转换成可以利用奇偶性进行优化的多项式。我们可以将原多项式,按照每一项

x

x

x 的幂次分组,偶数次幂一组,奇数次幂一组,然后再把奇数次幂组中的每一项提出一个

x

x

x 来,就像下面这样:

F

(

x

)

=

3

x

5

+

2

x

4

+

x

3

+

7

x

2

+

5

x

+

1

F

(

x

)

=

(

2

x

4

+

7

x

2

+

1

)

+

x

(

3

x

4

+

x

2

+

5

)

F(x)=3x^5+2x^4+x^3+7x^2+5x+1 Rightarrow \ F(x)=(2x^4+7x^2+1)+x(3x^4+x^2+5)

F(x)=3x5+2x4+x3+7x2+5x+1F(x)=(2x4+7x2+1)+x(3x4+x2+5)
我们可以把左右两部分视为新的,关于

x

2

x^2

x2 的多项式。

F

1

(

x

2

)

=

2

x

4

+

7

x

2

+

1

F

2

(

x

2

)

=

3

x

4

+

x

2

+

5

F_1(x^2) = 2x^4+7x^2+1 \ F_2(x^2) = 3x^4+x^2+5

F1(x2)=2x4+7x2+1F2(x2)=3x4+x2+5
我们只需要解决这两个小多项式的求值,再利用

F

(

x

)

=

F

1

(

x

2

)

+

x

F

2

(

x

2

)

F(x)=F_1(x^2)+xF_2(x^2)

F(x)=F1(x2)+xF2(x2)

F

(

x

)

=

F

1

(

x

2

)

x

F

2

(

x

2

)

F(-x)=F_1(x^2)-xF_2(x^2)

F(x)=F1(x2)xF2(x2) 这两个公式代回即可。最终取的点的形式是:

{

F

(

x

1

)

,

F

(

x

1

)

,

F

(

x

2

)

,

.

.

.

.

,

F

(

x

n

/

2

)

,

F

(

x

n

/

2

)

}

{F(x_1), F(-x_1), F(x_2),….,F(x_{n/2}),F(-x_{n/2})}

{F(x1),F(x1),F(x2),....,F(xn/2),F(xn/2)}。而两个小多项式都只有原多项式一半的规模,所以这个方法的复杂度是

O

(

n

log

n

)

O(nlog{n})

O(nlogn)。但是这个方法还存在一个问题:我们想利用奇偶性和每次选取

x

=

x

0

x=x_0

x=x0

x

=

x

0

x=-x_0

x=x0 两个点,来减少选取点的次数,也就是说我们每次取的是成相反数的两个点。然而当问题变为关于

F

1

(

x

2

)

F_1(x^2)

F1(x2)

f

2

(

x

2

)

f_2(x^2)

f2(x2) 时,可以发现:

x

2

x^2

x2 永远为正,我们不再能取一对相反数了。而解决这个问题的方法,被人们视为

F

F

T

FFT

FFT的精髓。

3.单位复根与快速傅里叶变换

可以看出,问题出在多项式的定义域上。一旦问题的规模缩小,原多项式的定义域会失效。那么如果我们把多项式的定义域扩大到复数域呢?这次不从最简单的情况入手了,这次需要一个逆向思维。下面结合图来说明。

假设递归已经进行到了终点,也就是现在多项式只剩下一次项了。也就是说我们需要选取一个点来求值,那么这个点如何选取?显然我可以选

x

=

1

x=1

x=1。然后考虑上一层递归,我要选两个值,我希望他们俩是相反数,还希望他们俩的平方等于 1。显然我可以取

x

=

1

x=1

x=1

x

=

1

x=-1

x=1。再往上一层,左边是 1,和第一层一样,关键是右边的 -1。我同样希望上一层的两个数,他们互为相反数,而他们的平方等于 -1。这两个数显然是

i

i

i

i

-i

i
[数学]------[多项式]-------快速傅里叶变换(FFT)
到目前为止,我们得到了四个数:1,-1,i,-i 。如果我的原多项式是一个三阶多项式,我只需要取

x

x

x 等于这四个数字就可以。但如果我的多项式阶数更高,这四个数还不够。我需要至少

N

N

N 个才行。这

N

N

N 个数需要满足

x

N

=

1

x^N = 1

xN=1,而之前学过的代数知识告诉我们,

x

N

=

1

x^N=1

xN=1 在复数域内恰好有

N

N

N 个解。也就是说我们前面得到的四个数字,其实是

x

4

=

1

x^4=1

x4=1 的四个解。

x

n

=

1

x^n=1

xn=1 在复数域上的

n

n

n 个解,均匀地分布在复平面的单位圆上。例如:

x

8

=

1

x^8=1

x8=1 的八个解,就是单位圆的八个八等分点。读者可以自行验证。这样的点叫做单位复根。将(1,0)点编号为0,其余点依次编号。我们定义

ω

n

k

omega_{n}^k

ωnk

x

n

=

1

x^n=1

xn=1 的第

k

k

k 号单位复根。
(本来想做个单位圆的图,然而没有比较好的几何工具,就咕咕了)

简单推导可得:

ω

n

k

=

cos

(

2

π

k

n

)

+

i

×

sin

(

2

π

k

n

)

omega_{n}^k=cos(frac{2pi k}{n})+i times sin(frac{2pi k}{n})

ωnk=cos(n2πk)+i×sin(n2πk)

单位复根有三个重要性质,通过简单推导就可以得到:

  • ω

    n

    n

    =

    1

    omega_{n}^n=1

    ωnn=1

  • ω

    n

    k

    =

    ω

    2

    n

    2

    k

    omega_{n}^k=omega_{2n}^{2k}

    ωnk=ω2n2k

  • ω

    2

    n

    k

    +

    n

    =

    ω

    2

    n

    k

    omega_{2n}^{k+n}=-omega_{2n}^k

    ω2nk+n=ω2nk

那么原来的公式

F

(

x

)

=

F

1

(

x

2

)

+

x

F

2

(

x

2

)

F(x)=F_1(x^2)+xF_2(x^2)

F(x)=F1(x2)+xF2(x2),利用单位复根的性质可得:
[数学]------[多项式]-------快速傅里叶变换(FFT)
同理可得:
[数学]------[多项式]-------快速傅里叶变换(FFT)
因此我们求出了

D

F

T

(

G

(

ω

n

/

2

k

)

)

DFT(G(omega^k_{n/2}))

DFT(G(ωn/2k))

D

F

T

(

H

(

ω

n

/

2

k

)

)

DFT(H(omega^k_{n/2}))

DFT(H(ωn/2k)) 后,就可以同时求出

D

F

T

(

f

(

ω

n

k

)

)

DFT(f(omega^k_n))

DFT(f(ωnk))

D

F

T

(

f

(

ω

n

k

+

n

/

2

)

)

DFT(f(omega^{k+n/2}_n))

DFT(f(ωnk+n/2)) 。于是对

G

G

G

H

H

H 分别递归 DFT 即可。考虑到分治 DFT 能处理的多项式长度只能是二的整数次幂,否则在分治的时候左右不一样长,右边就取不到系数了。所以要在第一次 DFT 之前就把序列向上补成长度为

2

m

2^m

2m(高次系数补

0

0

0)、最高项次数为

2

m

1

2^m-1

2m1 的多项式。

所以到现在我们可以发现,最初的原多项式代入的

n

n

n 个点,实际就是

{

ω

n

0

,

ω

n

1

,

ω

n

2

,

.

.

.

,

ω

n

n

1

}

{omega_n^0,omega_n^1,omega_n^2,…,omega_n^{n-1}}

{ωn0,ωn1,ωn2,...,ωnn1},其中

n

n

n 是二的整数次幂。

得到两个多项式乘积后,还需要把这个结果多项式从点值表示变回系数表示。

4.点值到系数的转换(IDFT)

这里我真写不下去了,因为 DFT 和 IDFT 怎么结合到一起的我真啥也没看懂,就知道一个结论:做 DFT 和 IDFT 的区别就在于一个正负号。

1

ω

k

=

ω

k

1

=

e

2

π

i

k

=

cos

(

2

π

k

)

+

i

×

sin

(

2

π

k

)

frac{1}{omega_k}=omega_k^{-1}=e^{-frac{2pi i}{k}}=cos(frac{2pi}{k})+i times sin(-frac{2pi}{k})

ωk1=ωk1=ek2πi=cos(k2π)+i×sin(k2π)
因此实际在做

F

F

T

FFT

FFT时,通常会传入一个参数,来决定是做

D

F

T

DFT

DFT 还是

I

D

F

T

IDFT

IDFT

之后要是弄懂了可能会回来写。

板子

P3803 【模板】多项式乘法(FFT).

P1919 【模板】A*B Problem升级版.


程序员灯塔
转载请注明原文链接:[数学]——[多项式]——-快速傅里叶变换(FFT)
喜欢 (0)