算法学习--快速傅里叶变换笔记

CHSHU
69 阅读2分钟

用途: 计算多项式卷积。

A(x)=a0+a1x+a2x2+⋅⋅⋅+an−1xn−1A(x)=a_{0}+a_{1}x+a_{2}x^{2}+···+a_{n-1}x^{n-1}

B(x)=b0+b1x+b2x2+⋅⋅⋅+bn−1xn−1B(x)=b_{0}+b_{1}x+b_{2}x^{2}+···+b_{n-1}x^{n-1}

求 : A(x)B(x)A(x)B(x)

前置芝士

多项式性质: 用任意 nn 个多项式函数图像上的点，可唯一确定一个 n−1n-1 次多项式。

复数,不会的话右转数学选修2-2.

做法

将 A(x)B(x)A(x)B(x) 的系数表示法转化为点表示法，

即 (x1,A(x1))   (x2,A(x2)) ⋅⋅⋅ (xn,A(xn))(x_{1},A(x_{1}))\ \ \ (x_{2},A(x_{2}))\ ···\ (x_{n},A(x_{n})) 和 (x1,B(x1))   (x2,B(x2)) ⋅⋅⋅ (xn,B(xn))(x_{1},B(x_{1}))\ \ \ (x_{2},B(x_{2}))\ ···\ (x_{n},B(x_{n}))

那么 A(x)B(x)A(x)B(x) 就可以表示为

(x1,A(x1)B(x1))   (x2,A(x2)B(x2)) ⋅⋅⋅ (xn,A(xn)B(xn))(x_{1},A(x_{1})B(x_{1}))\ \ \ (x_{2},A(x_{2})B(x_{2}))\ ···\ (x_{n},A(x_{n})B(x_{n}))

再将其转化为系数表示法。

下面我们考虑如何转化为系数表示法。

点的横坐标可以选择复数域上的单位根。

将复平面中以原点为圆心的单位圆等分为 nn 份，以 ωnk\omega _{n}^{k} 表示从 xx 正半轴一次逆时针取 nn 份中的 kk 份， ωnk\omega _{n}^{k} 就被称为 nn 次单位根。

nn 次单位根的性质:

  • ωni≠ωnj   ∀i≠j\omega ^{i}_{n} \ne \omega^{j}_{n} \ \ \ \forall i\ne j
  • ωnk=cos⁡2kπn+sin⁡2πni\omega^{k}_{n}=\cos \frac{2k\pi}{n}+\sin \frac{2\pi}{n}i
  • ωn0=ωnn=1\omega^{0}_{n}=\omega^{n}_{n}=1
  • ω2n2k=ωnk\omega^{2k}_{2n}=\omega^{k}_{n}
  • ωnk+n2=−ωnk\omega^{k+\frac{n}{2}}_{n}=-\omega^{k}_{n}

那么点表示法选取的横坐标为: ωn0  ωn1  ωn2⋅⋅⋅ωnn−1\omega^{0}_{n}\ \ \omega^{1}_{n}\ \ \omega^{2}_{n}···\omega^{n-1}_{n}

现在，需要快速计算下面这个式子:

(ωn0,A(ωn0)   (ωn1,A(ωn1)) ⋅⋅⋅ (ωnn−1,A(ωnn−1))(\omega^{0}_{n},A(\omega^{0}_{n})\ \ \ (\omega^{1}_{n},A(\omega^{1}_{n}))\ ···\ (\omega^{n-1}_{n},A(\omega^{n-1}_{n}))

快速傅里叶正变换(系数表示法转化为点表示法)

A(x)=a0+a1x+a2x2+⋅⋅⋅+an−1xn−1   n为2的整次幂=(a0+a2x2+a4x4+⋅⋅⋅+an−2xn−2)+(a1x+a3x3+a5x5+⋅⋅⋅an−1xn−1)   奇偶次项分开\begin{aligned} A(x)&=a_{0}+a_{1}x+a_{2}x^{2}+···+a_{n-1}x^{n-1}\ \ \ n为2的整次幂\\ &=(a_{0}+a_{2}x^{2}+a_{4}x^{4}+···+a_{n-2}x^{n-2})+(a_{1}x+a_{3}x^{3}+a_{5}x^{5}+···a_{n-1}x^{n-1})\ \ \ 奇偶次项分开 \end{aligned}

设:

A1(x)=a0+a2x+a4x2+⋅⋅⋅+an−2xn2−1A_{1}(x)=a_{0}+a_{2}x+a_{4}x^{2}+···+a_{n-2}x^{\frac{n}{2}-1}

A2(x)=a1+a3x+a5x2+⋅⋅⋅+an−1xn2−1A_{2}(x)=a_{1}+a_{3}x+a_{5}x^{2}+···+a_{n-1}x^{\frac{n}{2}-1}

那么可以得到:

A(x)=A1(x2)+xA2(x2)A(x)=A_{1}(x^2)+xA_{2}(x^2)

k∈[0,n−1]k\in [0,n-1] , 将值域分为两段 [0,n2−1] , [n2,n−1][0,\frac{n}{2}-1]\ ,\ [\frac{n}{2},n-1]

如果 k∈[0,n2−1]k\in[0,\frac{n}{2}-1] , 那么有 (ωnk)2=ωn2k(\omega^{k}_{n})^{2}=\omega^{2k}_{n}

A(ωnk)=A1(ωn2k)+ωnkA2(ωn2k)=A1(ωn2k)+ωnkA2(ωn2k)\begin{aligned} A(\omega^{k}_{n})&=A_{1}(\omega^{2k}_{n})+\omega^{k}_{n}A_{2}(\omega^{2k}_{n})\\ &=A_{1}(\omega^{k}_{\frac{n}{2}})+\omega^{k}_{n}A_{2}(\omega^{k}_{\frac{n}{2}})\\ \end{aligned}

如果 [n2,n−1][\frac{n}{2},n-1] , 将这一段全部减去 n2\frac{n}{2} , kk 还是遍历 [0,n2−1][0,\frac{n}{2}-1] 将 kk 加上 n2\frac{n}{2} 就行了。

A(ωnk+n2)=A1(ωn2k+n)+ωnk+n2A2(ωn2k+n)=A1(ωn2k)−ωnkA2(ωn2k)=A1(ωn2k)−ωnkA2(ωn2k)\begin{aligned} A(\omega^{k+\frac{n}{2}}_{n})&=A_{1}(\omega^{2k+n}_{n})+\omega^{k+\frac{n}{2}}_{n}A_{2}(\omega^{2k+n}_{n})\\ &=A_{1}(\omega^{2k}_{n})-\omega^{k}_{n}A_{2}(\omega^{2k}_{n})\\ &=A_{1}(\omega^{k}_{\frac{n}{2}})-\omega^{k}_{n}A_{2}(\omega^{k}_{\frac{n}{2}}) \end{aligned}

发现 A(ωnk) , A(ωnk+n2)A(\omega^{k}_{n})\ ,\ A(\omega^{k+\frac{n}{2}}_{n}) 计算非常相似，只要求出 A1(ωn2k) , A2(ωn2k)A_{1}(\omega^{k}_{\frac{n}{2}})\ ,\ A_{2}(\omega^{k}_{\frac{n}{2}}) 即可。

这个过程可以迭代 ，每次将区间分为两段 ，最多会有 log nlog\ n 层 ，总复杂度 O(nlog n)O(nlog\ n) .

迭代划分(假设有 232^3 项)

初始序列                         [x0,x1,x2,x3,x4,x5,x6,x7]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_{0},x_{1},x_{2},x_{3},x_{4},x_{5},x_{6},x_{7}]

奇偶划分一次                  [x0,x2,x4,x6],[x1,x3,x5,x7]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_{0},x_{2},x_{4},x_{6}],[x_{1},x_{3},x_{5},x_{7}]

将这两段继续奇偶划分     [x0,x4],[x2,x6],[x1,x5],[x3,x7]\ \ \ \ [x_{0},x_{4}],[x_{2},x_{6}],[x_{1},x_{5}],[x_{3},x_{7}]

接着分                           [x0],[x4],[x2],[x6],[x1],[x5],[x3],[x7]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_{0}],[x_{4}],[x_{2}],[x_{6}],[x_{1}],[x_{5}],[x_{3}],[x_{7}]

发现最后一层每一位下标是第一层下标二进制反过来(比如第二项 (1)10=(001)2 , (4)10=(100)2(1)_{10}=(001)_2\ ,\ (4)_{10}=(100)_2)

我们可以递推计算有 2bit2^{bit} 项的多项式最后一层的第 ii 位下标 chgichg_{i}。

chgi=(chgi22)∣((i mod 2)∗2bit−1)chg_{i}=(\frac{chg_{\frac{i}{2}}}{2})|((i\ mod\ 2)*2^{bit-1})

也就是

chg[i]=(chg[i>>1]>>1)|((i&1)<<(bit-1));

理解一下就是: 对于一个数ii，找到没有 ii 的第 00 位的数的反转之后的结果，将这个结果右移一位去掉最高位，再把 ii 的第 00 位放到最高位，就实现了二进制反转.

for  examplefor\ \ example : 计算到 66 时，将 (6)10=(110)2(6)_{10}=(110)_{2} 的二进制表示去掉第 00 位，也就是 (3)10=(011)2(3)_{10}=(011)_{2} ，33 的结果之前处理过了，是 (110)2=(6)10(110)_{2}=(6)_{10}，(这个例子好失败啊)，将反转后的结果右移 11 位 变成了 (011)2(011)_{2} ,然后把 66 的第 00 位(还是 00) ，放到这个数的最高位，就变成了 (011)2=(3)10(011)_2=(3)_{10}。

快速傅里叶逆变换(点表示法转化为系数表示法)

有点 (ωn0,A(ωn0))  (ωn1,A(ωn1))⋅⋅⋅(ωnn−1,A(ωnn−1))(\omega^{0}_{n},A(\omega^{0}_{n}))\ \ (\omega^{1}_{n},A(\omega^{1}_{n}))···(\omega^{n-1}_{n},A(\omega^{n-1}_{n})).

求 A(x)=a0+a1x+a2x2+⋅⋅⋅+an−1xn−1A(x)=a_{0}+a_{1}x+a_{2}x^2+···+a_{n-1}x^{n-1}

设 yk=ωnky_{k}=\omega^{k}_{n} 则有

ck=∑i=0n−1yi(ωn−k)i=nakc_{k}=\sum_{i=0}^{n-1}y_{i}(\omega^{-k}_{n})^{i}=na_{k}

证明:

ck=∑i=0n−1yi(ωn−k)i=∑i=0n−1(∑j=0n−1aj(ωni)j))(ωn−k)i=∑i=0n−1(∑j=0n−1aj(ωnj)i(ωn−k)i=∑i=0n−1(∑j=0n−1aj(ωnj−k)i=∑j=0n−1aj(∑i=0n−1(ωnj−k)i)\begin{aligned} c_{k}&=\sum_{i=0}^{n-1}y_{i}(\omega^{-k}_{n})^{i}\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(\omega^{i}_{n})^j))(\omega^{-k}_{n})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(\omega^{j}_{n})^i(\omega^{-k}_{n})^i\\ &=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_{j}(\omega^{j-k}_{n})^i\\ &=\sum_{j=0}^{n-1}a_{j}(\sum_{i=0}^{n-1}(\omega^{j-k}_{n})^i) \end{aligned}

设 S(x)=1+x+x2+⋅⋅⋅+xn−1S(x)=1+x+x^2+···+x^{n-1} .

k≠0k\ne 0 时 , 有：

S(ωnk)=ωn0+ωnk+ωn2k+⋅⋅⋅+ωn(n−1)kS(\omega^{k}_{n})=\omega^{0}_{n}+\omega^{k}_{n}+\omega^{2k}_{n}+···+\omega^{(n-1)k}_{n}

ωnkS(ωnk)=ωnk+ωn2k+ωn3k+⋅⋅⋅+ωn0\omega^{k}_{n}S(\omega^{k}_{n})=\omega^{k}_{n}+\omega^{2k}_{n}+\omega^{3k}_{n}+···+\omega^{0}_{n}

两柿相减,有：

(1−ωnk)S(ωnk)=0(1-\omega^{k}_{n})S(\omega^{k}_{n})=0

故 S(ωnk)=0S(\omega^{k}_{n})=0

k=0k=0 时 , S(ωnk)=nS(\omega^{k}_{n})=n ,代入显然

综上 : 把 SS 代入柿子可知：

ci=naic_{i}=na_{i}

证毕.

设 S(x)=1+x+x2+⋅⋅⋅+xn−1S(x)=1+x+x^2+···+x^{n-1} .

设 B(x)=y0+y1x1+y2x2+⋅⋅⋅+yn−1xn−1B(x)=y_{0}+y_{1}x^1+y_{2}x^2+···+y_{n-1}x^{n-1}.

则有 ci=B(ωn−i)c_{i}=B(\omega^{-i}_{n}) ,然后用上面说的正变换求一下就行。

FFT板子

#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define N 3000005
#define LL long long 
#define LB long double
using namespace std;

int n,m,bit,tot,chg[N];
const LB pi=3.1415926535897932384626433832795028841;//圆周率也能写 acos(-1)
struct complex
{
	LB x,y;
	complex operator+ (const complex &tmp) const//复数加
	{
		return (complex){x+tmp.x,y+tmp.y};
	}
	complex operator- (const complex &tmp) const//复数减
	{
		return (complex){x-tmp.x,y-tmp.y};
	}
	complex operator* (const complex &tmp) const//复数乘
	{
		return (complex){x*tmp.x-y*tmp.y,x*tmp.y+y*tmp.x};
	}
}a[N],b[N];

inline int qr()
{
	char ch=0;int w=1,x=0;
	while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
	return x*w;
}

inline void FFT(complex a[],int opt)
{
	for(register int i=0;i<tot;i++)
		if(i<chg[i])
			swap(a[i],a[chg[i]]);//找到a[i]奇偶位置转化后的值
	for(register int mid=1;mid<tot;mid<<=1)//mid表示现在区间长度的一半，区间长度len=mid<<1
	{
		complex w1=(complex){cos((LB)pi/mid),(LB)opt*sin((LB)pi/mid)};//计算w1;
		for(register int i=0;i<tot;i+=(mid<<1))//计算这一层所有长度为len的区间
		{
			complex wk=(complex){1,0};//wk初始值为w0=(1,0),w[k]=w[k-1]*w1,每次乘上w1就是要求的wk.
			for(register int j=0;j<mid;j++,wk=wk*w1)
			{
				complex x=a[i+j];//所有偶数项
				complex y=wk*a[i+mid+j];//所有奇数项乘wk
				a[i+j]=x+y;//迭代
				a[i+j+mid]=x-y;
			}
		}
	}
}

int main()
{
	n=qr();
	m=qr();
	for(register int i=0;i<=n;i++)
		a[i].x=qr();
	for(register int i=0;i<=m;i++)
		b[i].x=qr();
	while((1<<bit)<m+n+1)//计算一下 log n+m 向上取整
		bit++;
	tot=(1<<bit);
	for(register int i=0;i<tot;i++)//递推第bit层奇偶项转换位置
		chg[i]=(chg[i>>1]>>1)|((i&1)<<(bit-1));
	FFT(a,1);//系数表示法转化为点表示法
	FFT(b,1);
	for(register int i=0;i<tot;i++)
		a[i]=a[i]*b[i];//计算点表示法纵坐标
	FFT(a,-1);//点表示法转化为系数表示法
	for(register int i=0;i<m+n+1;i++)
		printf("%d ",(int)(a[i].x/tot+0.5));//为了防止精度问题加上0.5再向下取整
	printf("\n");
	return 0;
}
avatar
文章
阅读
粉丝