Fast Fourier Transportation(FFT)
·多项式的表达
系数表达
对于一个次数界为n的多项式 A ( x ) = ∑ j = 0 n − 1 a j x j A(x)=\sum_{j=0}^{n-1}{a_jx^j} A(x)=∑j=0n−1ajxj而言,其系数表达是由一个系数组成的向量 a = ( a 0 , a 1 , . . . , a n − 1 ) a=(a_0,a_1,...,a_{n-1}) a=(a0,a1,...,an−1)。
点值表达
一个次数界为n的多项式A(x)的点值表达就是一个由n个点值对所组成的集合
(
x
0
,
y
0
)
,
(
x
1
,
y
1
)
,
.
.
.
,
(
x
n
−
1
,
y
n
−
1
)
{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})}
(x0,y0),(x1,y1),...,(xn−1,yn−1)
使得对k=0,1,…,n-1,所有
x
k
x_k
xk各不相同,
y
k
=
A
(
x
k
)
y_k=A(x_k)
yk=A(xk)
简单的求点值运算我们可以随意代入n个不相同的数,然后得出点对,时间复杂度
Θ
(
n
2
)
\Theta(n^2)
Θ(n2)。后面可以看到,如果我们用一点巧妙的取值,可以使时间复杂度优化到
Θ
(
n
lg
2
n
)
\Theta(n\lg_2n)
Θ(nlg2n)。
求值计算的逆(从一个多项式的点值表达确定的系数表达形式)称为插值。
·多项式运算
C j = ∑ k = 0 j A k B j − k C_j=\sum_{k=0}^{j}{A_kB_{j-k}} Cj=k=0∑jAkBj−k
上述是多项式的乘法,我们把C称为A和B的卷积(convolution),表示成 C = A ⨂ B C=A\bigotimes B C=A⨂B。
FFT的主要思路是首先把A和B转成点值表达,然后得到C的点值表达,再逆着做一遍,得到C的系数表达。
·DFT与FFT
上述做法太慢,我们要用 Θ ( n 2 ) \Theta(n^2) Θ(n2)的时间把每个多项式转成点值表达,然后再用 Θ ( n 2 ) \Theta(n^2) Θ(n2)的时间转回来,明显很慢,还不如暴力。
我们想,可不可以使用某些特殊的数,使得每次可以做一次运算就可以得到多个数的呢?答案是有的:单位根复数根
n次单位复数根是满足
ω
n
=
1
\omega^n=1
ωn=1的复数
ω
\omega
ω。n次单位复数根恰好有n个,对于k=0,1,…,n-1,这些根是
e
π
i
k
/
n
e^{\pi ik/n}
eπik/n。为了解释这个表达式,我们用复数的指数形式来定义:
e
i
u
=
c
o
s
(
u
)
+
i
s
i
n
(
u
)
e^{iu}=cos(u)+isin(u)
eiu=cos(u)+isin(u)
也就是给定一个单位圆,上面均匀地分布着n个向量,如图:
·关于n次单位复根
以上图为例我们可以发现,每一个n(这里是8)次单位复根都是一个向量,他们在乘法意义下形成一个群。
引理1:(消去引理)
对
于
任
意
整
数
n
⩾
0
,
k
⩾
0
,
以
及
d
>
0
,
对于任意整数n\geqslant0,k\geqslant0,以及d>0,
对于任意整数n⩾0,k⩾0,以及d>0,
ω d k d k = ω n k \omega^{dk}_{dk}=\omega^{k}_{n} ωdkdk=ωnk
证明
ω
d
k
d
k
=
(
e
2
π
i
/
d
n
)
d
k
=
(
e
2
π
i
/
n
)
k
=
ω
n
k
\omega^{dk}_{dk}=(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k=\omega^{k}_{n}
ωdkdk=(e2πi/dn)dk=(e2πi/n)k=ωnk
引理2:(折半引理)
如
果
n
>
0
为
偶
数
,
那
么
n
个
n
次
单
位
根
的
平
方
集
合
就
是
n
/
2
个
n
/
2
次
单
位
根
的
集
合
如果n>0为偶数,那么n个n次单位根的平方集合就是n/2个n/2次单位根的集合
如果n>0为偶数,那么n个n次单位根的平方集合就是n/2个n/2次单位根的集合
证明
(
ω
n
k
+
n
/
2
)
2
=
ω
n
2
k
+
n
=
ω
n
2
k
ω
n
n
=
(
ω
n
k
)
2
(\omega^{k+n/2}_{n})^2=\omega^{2k+n}_n=\omega^{2k}_n\omega^n_n=(\omega^k_n)^2
(ωnk+n/2)2=ωn2k+n=ωn2kωnn=(ωnk)2
因此
ω
n
k
\omega^k_n
ωnk与
ω
n
k
+
n
/
2
\omega^{k+n/2}_n
ωnk+n/2平方相等。
引理3:(求和引理)
对
于
任
意
整
数
n
⩾
0
和
不
能
被
n
整
除
的
非
负
整
数
k
,
有
对于任意整数n\geqslant0和不能被n整除的非负整数k,有
对于任意整数n⩾0和不能被n整除的非负整数k,有
∑ j = 0 n − 1 ( ω n k ) j = 0 \sum_{j=0}^{n-1}(\omega^k_n)^j=0 j=0∑n−1(ωnk)j=0
证明
∑
j
=
0
n
−
1
(
ω
n
k
)
j
=
(
ω
n
k
)
0
(
1
−
(
ω
n
k
)
n
)
1
−
ω
n
k
=
(
ω
n
n
)
k
−
1
ω
n
k
−
1
=
(
1
)
k
−
1
ω
n
k
−
1
=
0
\sum_{j=0}^{n-1}(\omega^k_n)^j=\frac{(\omega^k_n)^0(1-(\omega^k_n)^n)}{1-\omega^{k}_{n}}=\frac{(\omega^n_n)^k-1}{\omega^{k}_{n}-1}=\frac{(1)^k-1}{\omega^{k}_{n}-1}=0
j=0∑n−1(ωnk)j=1−ωnk(ωnk)0(1−(ωnk)n)=ωnk−1(ωnn)k−1=ωnk−1(1)k−1=0
因为要求k不能被n整除,而且仅当k被n整除时
ω
n
k
=
1
\omega^k_n=1
ωnk=1成立,同时保证分母不为0。
DFT
回顾一下,我们希望计算次数界为n的多项式
A
(
x
)
=
∑
j
=
0
n
−
1
a
j
x
j
A(x)=\sum_{j=0}^{n-1}{a_jx^j}
A(x)=j=0∑n−1ajxj
在
ω
n
0
,
ω
n
1
,
.
.
.
,
ω
n
n
−
1
\omega^0_n,\omega^1_n,...,\omega^{n-1}_n
ωn0,ωn1,...,ωnn−1处的值。假设A以系数形式给出,接下来定义结果
y
k
y_k
yk:
y
k
=
A
(
ω
n
k
)
=
∑
j
=
0
n
−
1
a
j
ω
n
k
j
y_k=A(\omega^k_n)=\sum_{j=0}^{n-1}a_j\omega^{kj}_n
yk=A(ωnk)=j=0∑n−1ajωnkj
向量
y
=
(
y
0
,
y
1
,
.
.
.
,
y
n
−
1
)
y=(y_0,y_1,...,y_{n-1})
y=(y0,y1,...,yn−1)就是系数向量
a
=
(
a
0
,
a
1
,
.
.
.
,
a
n
−
1
)
a=(a_0,a_1,...,a_{n-1})
a=(a0,a1,...,an−1)的离散傅里叶变换(DFT)。我们也记作
y
=
D
F
T
n
(
a
)
y=DFT_n(a)
y=DFTn(a)。
FFT
通过使用一种称为快速傅里叶变换(FFT)的方法,利用复根的特殊性质,我们就可以在 Θ ( n lg n ) \Theta(n\lg n) Θ(nlgn)的时间内计算出 D F T n ( a ) DFT_n(a) DFTn(a)。
注意:通篇的n我们假设是2的整数次幂。
FFT利用分治策略,采用A(x)中偶数下标的系数与奇数下标的系数,分别定义两个新的次数界为n/2的多项式
A [ 0 ] ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n / 2 − 1 A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1} A[0](x)=a0+a2x+a4x2+...+an−2xn/2−1
A [ 1 ] ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n / 2 − 1 A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1} A[1](x)=a1+a3x+a5x2+...+an−1xn/2−1
可以很容易发现
A
(
x
)
=
A
[
0
]
(
x
2
)
+
x
A
[
1
]
(
x
2
)
A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)
A(x)=A[0](x2)+xA[1](x2)
所以原问题转化为求两个次数界为n/2的多项式
A
[
0
]
(
x
)
A^{[0]}(x)
A[0](x)和
A
[
1
]
(
x
)
A^{[1]}(x)
A[1](x)在点
(
ω
n
0
)
2
,
(
ω
n
1
)
2
,
.
.
.
,
(
ω
n
n
−
1
)
2
(\omega^0_n)^2,(\omega^1_n)^2,...,(\omega^{n-1}_n)^2
(ωn0)2,(ωn1)2,...,(ωnn−1)2的取值。
所以我们可以发现在求出
A
[
0
]
(
x
2
)
A^{[0]}(x^2)
A[0](x2)和
A
[
1
]
(
x
2
)
A^{[1]}(x^2)
A[1](x2)以后,可以算出两个复根的结果:
y
k
=
y
k
[
0
]
+
ω
n
k
y
k
[
1
]
=
A
[
0
]
(
ω
n
2
k
)
+
ω
n
k
A
[
1
]
(
ω
n
2
k
)
=
A
(
ω
n
k
)
y_k=y^{[0]}_k+\omega^k_ny^{[1]}_k =A^{[0]}(\omega^{2k}_n)+\omega^k_nA^{[1]}(\omega^{2k}_n) =A(\omega^k_n)
yk=yk[0]+ωnkyk[1]=A[0](ωn2k)+ωnkA[1](ωn2k)=A(ωnk)
还有
y
k
+
(
n
/
2
)
=
y
k
[
0
]
−
ω
n
k
y
k
[
1
]
=
y
k
[
0
]
+
ω
n
k
+
(
n
/
2
)
y
k
[
1
]
=
A
[
0
]
(
ω
n
2
k
)
+
ω
k
+
(
n
/
2
)
A
[
1
]
(
ω
n
2
k
)
y_{k+(n/2)}=y^{[0]}_k-\omega^{k}_ny^{[1]}_k =y^{[0]}_k+\omega^{k+(n/2)}_ny^{[1]}_k =A^{[0]}(\omega^{2k}_n)+\omega^{k+(n/2)}A^{[1]}(\omega^{2k}_n)
yk+(n/2)=yk[0]−ωnkyk[1]=yk[0]+ωnk+(n/2)yk[1]=A[0](ωn2k)+ωk+(n/2)A[1](ωn2k)
=
A
[
0
]
(
ω
n
2
k
+
n
)
+
ω
k
+
(
n
/
2
)
A
[
1
]
(
ω
n
2
k
+
n
)
=
A
(
ω
n
k
+
(
n
/
2
)
)
=A^{[0]}(\omega^{2k+n}_n)+\omega^{k+(n/2)}A^{[1]}(\omega^{2k+n}_n) =A(\omega^{k+(n/2)}_n)
=A[0](ωn2k+n)+ωk+(n/2)A[1](ωn2k+n)=A(ωnk+(n/2))
所以就有代码:
void FFT(comp *a,int n,int inv){
if(n==1) return;
int mid=n/2;
for (int i=0;i<mid;++i) c[i]=a[i*2],c[i+mid]=a[i*2+1];
for (int i=0;i<n;++i) a[i]=c[i];
FFT(a,mid,inv);
FFT(a+mid,mid,inv);
comp wn={cos(2.0*pi/n),inv*sin(2.0*pi/n)},w={1,0};
for (int i=0;i<mid;++i,w=w*wn){
c[i]=a[i]+w*a[i+mid];
c[i+mid]=a[i]-w*a[i+mid];
}
for (int i=0;i<n;++i) a[i]=c[i];
}
·在单位复数根的插值
现在我们展示如何在单位复数根处插值来完成多项式乘法方案,使得我们把一个多项式从点值表达转换回系数表达。
我们可以把DFT写成矩阵乘积
[
y
0
y
1
y
2
⋮
y
n
−
1
]
=
[
1
1
1
1
⋯
1
1
ω
n
ω
n
2
ω
n
3
⋯
ω
n
n
−
1
1
ω
n
2
ω
n
4
ω
n
6
⋯
ω
n
2
(
n
−
1
)
1
ω
n
3
ω
n
6
ω
n
9
⋯
ω
n
3
(
n
−
1
)
⋮
⋮
⋮
⋮
⋱
⋮
1
ω
n
n
−
1
ω
n
2
(
n
−
1
)
ω
n
3
(
n
−
1
)
⋯
ω
n
(
n
−
1
)
(
n
−
1
)
]
[
a
0
a
1
a
2
⋮
a
n
−
1
]
\left[ \begin{matrix} y_0\\ y_1\\ y_2\\ \vdots\\ y_{n-1} \end{matrix} \right]= \left[ \begin{matrix} 1 & 1 & 1 & 1 &\cdots& 1\\ 1 & \omega_n & \omega^2_n & \omega^3_n &\cdots& \omega^{n-1}_n\\ 1 & \omega^2_n & \omega^4_n & \omega^6_n &\cdots& \omega^{2(n-1)}_n\\ 1 & \omega^3_n & \omega^6_n & \omega^9_n &\cdots& \omega^{3(n-1)}_n\\ \vdots & \vdots& \vdots& \vdots& \ddots& \vdots\\ 1 & \omega^{n-1}_n & \omega^{2(n-1)}_n & \omega^{3(n-1)}_n &\cdots& \omega^{(n-1)(n-1)}_n \end{matrix} \right] \left[ \begin{matrix} a_0\\ a_1\\ a_2\\ \vdots\\ a_{n-1} \end{matrix} \right]
⎣⎢⎢⎢⎢⎢⎡y0y1y2⋮yn−1⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡1111⋮11ωnωn2ωn3⋮ωnn−11ωn2ωn4ωn6⋮ωn2(n−1)1ωn3ωn6ωn9⋮ωn3(n−1)⋯⋯⋯⋯⋱⋯1ωnn−1ωn2(n−1)ωn3(n−1)⋮ωn(n−1)(n−1)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡a0a1a2⋮an−1⎦⎥⎥⎥⎥⎥⎤
尴尬的是跑得贼慢:
随便卡卡就爆了…
分治难免递归,递归常数大。
于是,考虑改进。
·蝴蝶变换
盗图一张
可以发现,每个下标的二进制形式反过来就是它们最后在序列中的位置。
于是就有了迭代打法。
void FFT(Moon *a,int inv){
int i,j,len;
for (i=0;i<n;++i)
if(i<R[i])swap(a[i],a[R[i]]);
for (len=2;len<=n;len<<=1){
int half=len/2;
Moon w,wn={cos(Pi/half),inv*sin(Pi/half)};
for (j=0;j<n-i;j+=len,w={1,0}){
for (i=0;i<half;++i,w=w*wn){
Moon q=w*a[j+half+i],Q=a[j+i];
a[j+half+i]=Q-q;
a[j+i]=Q+q;
}
}
}
}
int main(){
for (i=0;i<n;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(p-1));
}