快速傅里叶变换在算法竞赛中是能高效的计算多项式的算法,可以从 O ( n 2 ) O(n^2) O(n2)的时间复杂度,通过卷积的思想优化到 O ( n l o g n ) O(n logn) O(nlogn),当然在通信方面也有很大的作用。
卷积
对于多项式乘法来说现在给你两个多项式,请你求出他们的和。
例如:
F:
(
a
+
b
x
2
+
c
x
3
…
…
+
d
x
k
)
(a+bx^2+cx^3……+dx^k)
(a+bx2+cx3……+dxk)
P:
(
e
x
+
f
x
2
+
g
x
3
…
…
+
h
x
k
)
(ex+fx^2+gx^3……+hx^k)
(ex+fx2+gx3……+hxk)
那么我们假定a[],b[]分别代表了多项式F和多项式P的系数,数组代表他们的指数
那么我们的程序就这么写
for(int i=1;i<=k;i++){//枚举次数
for(int j=0;j<=i;j++){
ans[i]+=a[j]*b[i-j];
}
}
其实不难发现这个算法的时间复杂度是
O
(
n
2
)
O(n^2)
O(n2)的远远高于我们我们预期的时间要求。
我们当前求的东西其实也叫卷积
a
n
s
[
i
]
=
∑
i
?
j
=
=
k
a
[
i
]
b
[
j
]
ans[i]=\displaystyle\sum_{i?j==k}^{}a[i]b[j]
ans[i]=i?j==k∑a[i]b[j]
其中?为i和j之间的某种运算
两种表示法
对于一个多项式来说我们可以写成
∑
i
=
0
n
−
1
a
[
i
]
x
i
\displaystyle\sum_{i=0}^{n-1}a[i]x^i
i=0∑n−1a[i]xi,它有n个系数
a
[
i
]
a[i]
a[i],这是系数表示法。
同时如果在这个多项式上取
n
+
1
n+1
n+1(只有多取一个点才会确定这一个多项式,否则不是唯一的解),选出来的是:
(
x
0
,
F
(
x
0
)
)
,
(
x
1
,
F
(
x
1
)
)
.
.
.
(
x
n
,
F
(
x
n
)
)
(x_0,F(x_0)),(x_1,F(x_1))...(x_n,F(x_n))
(x0,F(x0)),(x1,F(x1))...(xn,F(xn))
同时如果多项式P的点值表示法为
(
x
0
,
P
(
x
0
)
)
,
(
x
1
,
P
(
x
1
)
)
.
.
.
(
x
n
,
P
(
x
n
)
)
(x_0,P(x_0)),(x_1,P(x_1))...(x_n,P(x_n))
(x0,P(x0)),(x1,P(x1))...(xn,P(xn))
那么我们很容易得到:
P
∗
F
P*F
P∗F的点值表示法是
(
x
0
,
P
(
x
0
)
∗
F
(
x
0
)
)
,
(
x
1
,
P
(
x
1
)
∗
F
(
x
1
)
)
.
.
.
(
x
n
,
P
(
x
n
)
∗
F
(
x
n
)
)
(x_0,P(x_0)*F(x_0)),(x_1,P(x_1)*F(x_1))...(x_n,P(x_n)*F(x_n))
(x0,P(x0)∗F(x0)),(x1,P(x1)∗F(x1))...(xn,P(xn)∗F(xn))
那么
P
(
x
)
∗
F
(
x
)
=
P
∗
F
(
x
)
P(x)*F(x)=P*F(x)
P(x)∗F(x)=P∗F(x)
但是题目一般会要求我们用系数表示法转换成为点值表示法,这里需要用到DFT,如果是点值表示法转换成为系数表示法,那么就是IDFT。
而 FFT 就是通过取某些特殊的x的点值来加速 DFT 和 IDFT 的过程,那么x的点值其实就是我们接下来要介绍的单位复根。
单位复根
单位复根: x n = 1 x^n=1 xn=1的复数解,在复平面上其实就是单位圆上的全部点。
解决了刚刚上面的问题,那么我们现在的问题又来了。
P ∗ F P*F P∗F的点值表示法是 ( x 0 , P ( x 0 ) ∗ F ( x 0 ) ) , ( x 1 , P ( x 1 ) ∗ F ( x 1 ) ) . . . ( x n , P ( x n ) ∗ F ( x n ) ) (x_0,P(x_0)*F(x_0)),(x_1,P(x_1)*F(x_1))...(x_n,P(x_n)*F(x_n)) (x0,P(x0)∗F(x0)),(x1,P(x1)∗F(x1))...(xn,P(xn)∗F(xn))
但是我们需要点值表示法再转换为系数表示法来完成我们的fft整个操作,高斯消元相对而言十分的麻烦,因此我们为了优化这一点,引入了单位复根
w
w
w,单位复根有着十分优良的性质,已知
w
k
=
1
w^k=1
wk=1,两个单位复根相乘,表示为模长相乘,辐角相加,那么我们就可以类似于理解为切蛋糕的想法,把这个蛋糕切成k等分,以上符合条件的部分就可以叫做复根,用
w
w
w表示。
(图源OIwiki 传送门)
同时,单位复根有着非常优良的性质:
ω n k = ω n k % n \large\omega_{n}^k=\omega_{n}^{k\%n} ωnk=ωnk%n
对任意的n, ω n 0 = 1 \ \ \large{\omega_{n}^{0}=1} ωn0=1
ω n k = ( ω n 1 ) k \ \large{\omega_{n}^k=(\omega_{n}^1)^k} ωnk=(ωn1)k
ω n k = ω n 1 ∗ ω n 1 ∗ . . . ∗ ω n 1 (共k个) \omega_{n}^k=\omega_{n}^1*\omega_{n}^1*...*\omega_{n}^1\text{(共k个)} ωnk=ωn1∗ωn1∗...∗ωn1(共k个)
两个n次单位根 ω n j ∗ ω n k = ω n j + k \ \large{\omega_{n}^j*\omega_{n}^k=\omega_{n}^{j+k}} ωnj∗ωnk=ωnj+k
( ω n k ) − 1 = ( ω n k ) − 1 = ω n − k = ω n n − k (ω^k_n)^{-1}=(ω^k_n)^{-1}=ω_n^{-k}=ω_n^{n-k} (ωnk)−1=(ωnk)−1=ωn−k=ωnn−k
ω n n − k ∗ ω n k = ω n n − k + k = ω n n = ω n 0 = 1 ω_n^{n-k}*ω^k_n=\omega_{n}^{n-k+k}=\omega_{n}^{n}=\omega_{n}^{0}=1 ωnn−k∗ωnk=ωnn−k+k=ωnn=ωn0=1
ω p n p k = ω n k ω^{pk}_{pn}=ω^k_n ωpnpk=ωnk
如果n是偶数, ω n ( k + n / 2 ) = − ω n k \ \large{ω^{(k+n/2)}_{n}=-ω^k_n} ωn(k+n/2)=−ωnk
DFT
DFT的基本思路是分而治之
如下
假设
f
(
x
)
=
1
+
2
x
+
x
2
+
2
x
3
+
x
4
+
2
x
5
f(x)=1+2x+x^2+2x^3+x^4+2x^5
f(x)=1+2x+x2+2x3+x4+2x5
那么
f
(
x
)
=
(
1
+
x
2
+
x
4
)
+
x
∗
(
2
+
2
x
2
+
2
x
4
)
f(x)=(1+x^2+x^4)+x*(2+2x^2+2x^4)
f(x)=(1+x2+x4)+x∗(2+2x2+2x4)
设
f
l
(
x
)
=
(
1
+
x
+
x
2
)
,
f
r
(
x
)
=
(
2
+
2
x
+
2
x
2
)
fl(x)=(1+x+x^2),fr(x)=(2+2x+2x^2)
fl(x)=(1+x+x2),fr(x)=(2+2x+2x2)
整理一下:
f
(
x
)
=
f
l
(
x
2
)
+
x
∗
f
r
(
x
2
)
f(x)=fl(x^2)+x*fr(x^2)
f(x)=fl(x2)+x∗fr(x2)
这个就是分治的一个基本流程
根据单位根的性质我们可以得到
f ( w n k ) = f l ( ( w n k ) 2 ) + w n k ∗ f r ( ( w n k ) 2 ) = f l ( w n / 2 k ) + w n k ∗ f r ( w n / 2 k ) f(w_n^k)=fl((w_n^k)^2)+w_n^k*fr((w_n^k)^2)=fl(w_{n/2}^k)+w_n^k*fr(w_{n/2}^k) f(wnk)=fl((wnk)2)+wnk∗fr((wnk)2)=fl(wn/2k)+wnk∗fr(wn/2k)
f ( w n k + n / 2 ) = f l ( ( w n k + n / 2 ) 2 ) + w n k + n / 2 ∗ f r ( ( w n k + n / 2 ) 2 ) = f l ( w n / 2 k ) − w n k ∗ f r ( w n / 2 k ) f(w_n^{k+n/2})=fl((w_n^{k+n/2})^2)+w_n^{k+n/2}*fr((w_n^{k+n/2})^2)=fl(w_{n/2}^k)-w_n^k*fr(w_{n/2}^k) f(wnk+n/2)=fl((wnk+n/2)2)+wnk+n/2∗fr((wnk+n/2)2)=fl(wn/2k)−wnk∗fr(wn/2k)
这两条式子在分治的情况下的作用十分大。
到这里,如果我们知道两个多项式 f l ( x ) fl(x) fl(x)和 f r ( x ) fr(x) fr(x)分别在 ω n / 2 0 , ω n / 2 1 , ω n / 2 2 , . . . , ω n / 2 n / 2 − 1 ω^0_{n/2},ω^1_{n/2},ω^2_{n/2},...,ω^{n/2-1}_{n/2} ωn/20,ωn/21,ωn/22,...,ωn/2n/2−1的点值表示,用 f ( ω n k ) = f l ( ω n / 2 k ) + ω n k f r ( ω n / 2 k ) f(ω^k_n)=fl(ω^k_{n/2})+ω^k_nfr(ω^k_{n/2}) f(ωnk)=fl(ωn/2k)+ωnkfr(ωn/2k)可以 O ( n ) O(n) O(n)求出 f ( x ) f(x) f(x)在 ω n 0 , ω n 1 , ω n 2 , . . . , ω n n / 2 − 1 ω^0_{n},ω^1_{n},ω^2_{n},...,ω^{n/2-1}_{n} ωn0,ωn1,ωn2,...,ωnn/2−1处的点值表示。
用
f
(
ω
n
k
+
n
/
2
)
=
f
l
(
ω
n
/
2
k
)
−
ω
n
k
f
r
(
ω
n
/
2
k
)
f(ω^{k+n/2}_n)=fl(ω^{k}_{n/2})-ω^k_nfr(ω^{k}_{n/2})
f(ωnk+n/2)=fl(ωn/2k)−ωnkfr(ωn/2k)可以
O
(
n
)
O(n)
O(n)求出
f
(
x
)
f(x)
f(x)在
ω
n
n
/
2
,
ω
n
n
/
2
+
1
,
ω
n
n
/
2
+
2
,
.
.
.
,
ω
n
n
−
1
ω^{n/2}_{n},ω^{n/2+1}_{n},ω^{n/2+2}_{n},...,ω^{n-1}_{n}
ωnn/2,ωnn/2+1,ωnn/2+2,...,ωnn−1处的点值表示。
所以如果我们知道两个多项式
f
l
(
x
)
fl(x)
fl(x)和
f
r
(
x
)
fr(x)
fr(x)分别在
ω
n
/
2
0
,
ω
n
/
2
1
,
ω
n
/
2
2
,
.
.
.
,
ω
n
/
2
n
/
2
−
1
ω^0_{n/2},ω^1_{n/2},ω^2_{n/2},...,ω^{n/2-1}_{n/2}
ωn/20,ωn/21,ωn/22,...,ωn/2n/2−1的点值表示,就可以
O
(
n
)
O(n)
O(n)求出
f
(
x
)
f(x)
f(x)在
ω
n
0
,
ω
n
1
,
ω
n
2
,
.
.
.
,
ω
n
n
−
1
ω^0_{n},ω^1_{n},ω^2_{n},...,ω^{n-1}_{n}
ωn0,ωn1,ωn2,...,ωnn−1处的点值表示。
我们就成功的获得了
f
(
x
)
f(x)
f(x)的点值标示。然后就可以递归分治了。
IDFT
我们知道IDFT其实就是DFT的一个逆变换,其实本质上也就是点值表示法到系数表示法的一个转换。
假设
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+a1x+a2x2+...+an−1xn−1
经过DFT变换以后就变成了点值序列A
将多项式做如下变换:
A
(
k
)
=
∑
i
=
0
n
−
1
f
(
i
)
w
n
i
k
A(k)=\displaystyle\sum_{i=0}^{n-1}f(i)w_n^{ik}
A(k)=i=0∑n−1f(i)wnik
经过推证:
n
∗
f
(
k
)
=
∑
i
=
0
n
−
1
(
w
n
−
k
)
i
A
(
i
)
n*f(k)=\displaystyle\sum_{i=0}^{n-1}(w_n^{-k})^iA(i)
n∗f(k)=i=0∑n−1(wn−k)iA(i)
简而言之:就是把DFT中的
w
n
1
w_n^1
wn1换成
w
n
−
1
w_n^{-1}
wn−1,搞完以后除以n就可以了。
蝴蝶变换
这是一个优化的方式:
因为在分治的过程中,我们一定会有很多的拷贝情况
那么分奇偶分治有如下的一个规律
1 2 3 4 5 6 7 8
1 3 5 7 | 2 4 6 8
1 5 | 3 7 | 2 6 | 4 8
1 | 5 | 3 | 7 | 2 | 6 | 4 | 8
对应的二进制为
1 :000 1:000
2: 001 5:100
3: 010 3:010
4: 011 7:110
5: 100 2:001
6: 101 6:101
7: 110 4:011
8: 111 8:111
正好是二进制位的倒转。
那么我们可以用一个递推的做法,来完成如下的操作。
for(int i=0;i<_max;i++){
a[i]=a[i>>1]>>1|((i&1)?_max>>1:0);
//前面的是除了第一位,根据递推来求其他位的反转
//后边的求我们最高位有没有。
//也就是如果原来的数第一位有数,那么现在就是最高位
}
代码实现
构造复数类
struct cp{
double x,y;
//构造函数
cp(double _x=0.0,double _y=0.0){
x=_x,y=_y;
}
cp operator + (cp const &c)const {
return cp(x+c.x,y+c.y);
}
cp operator - (cp const &c)const {
return cp(x-c.x,y-c.y);
}
cp operator * (cp const &c)const {
return cp(x*c.x-y*c.y,x*c.y+y*c.x);
}
};