快速计算多项式相乘系数 FFT快速傅里叶变换

xxcdsg

于 2023-07-07 20:20:46 发布

阅读量389

点赞数

分类专栏：算法，结构体封装计划数论文章标签：算法 c++ 快速傅里叶变换

本文链接：https://blog.csdn.net/xxcdsg/article/details/131603957

版权

算法，结构体封装计划同时被 2 个专栏收录

7 篇文章 0 订阅

订阅专栏

数论

2 篇文章 0 订阅

订阅专栏

快速计算多项式相乘系数 FFT快速傅里叶变换

快速傅里叶变换(FFT)——有史以来最巧妙的算法？

正常求两个多项式乘积

$A(x)=\sum_{i=0}^{n}{A_ix^i},B(x)=\sum_{i=0}^{n}{B_ix^i} \\ C(x)=\sum_{i=0}^{n}\sum_{j=0}^nA_iB_jx^{i+j}$

复杂度为 $O(n^2)$

分别设n+1个在A与B上的点

根据高斯消元法，n+1个不同的点一定能确定一个唯一的n次多项式

$A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2))...(x_n,A(x_n))\}\\ B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2))...(x_n,B(x_n))\}\\$

横坐标相同的点相乘

$C(x)=\{(x_0,A(x_0)B(x_0)),(x_1,A(x_1)B(x_1))...(x_n,A(x_n)B(x_n))\}$

那么就有两个问题：

1.怎么快速得到n+1个点

2.怎么将点值表达法变回系数

如果 $A (x)$ 是一个偶函数，那么 $A (x) = A (- x)$ ，因此带入一个值也就得到了两个点

同理如果 $A (x)$ 是一个奇函数，那么有 $A (x) = - A (- x)$ ，也可以带一个点得到两个点

因此可以将 $A (x)$ 分成奇函数部分与偶函数部分

$P(x)=a_0+a_1x+a_2x^2+...+a_nx^n\\ P(x)=(a_0+a_2x^2+a_4x^4+...)+x(a_1+a_3x^2+a_5x^4+...)\\ P(x)=P_e(x^2)+xP_o(x^2)\\ P(x_i)=P_e(x_i^2)+x_iP_o(x_i^2)\\ P(-x_i)=P_e(x_i^2)-x_iP_o(x_i^2)$

然后如果将 $P_e(x),P_o(x)$ 继续分解

但我们的 $P_e(x),P_o(x)$ 是通过 $P(x)=P_e(x^2)+xP_o(x^2)$ 得来的，因此取相反数也无法让它们一一符号相反

那么我们取复数呢？

$i^2=(-i)^2=-1\\ 1^2=(-1)^2=1\\ 因此有\\ i,-i,1,-1:\ x_1,-x_1,x_2,-x_2\\ -1,1:x_1^2,x_2^2\\ 1:x_1^4 \\可解决n\leq4的问题$

那么更大的呢？

我们令 $w=e^{\frac{2\pi i}{n}}=cos(\frac{2\pi}{n})+isin(\frac{2\pi}{n})$ ，我们可以得到它的性质正好可以解决该问题

请添加图片描述

反过来呢？

将 $w=e^{-\frac{2\pi i}{n}}$ ，然后给每个数都乘上 $\frac{1}{n}$ 即可

//初始化每个位置最终到达的位置 
void init(int k)
{
    int len=1<<k;
	for(int i=0;i<len;i++)
	rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
}
//a表示要操作的系数，n表示序列长度
//若flag为1，则表示FFT，为-1则为IFFT(需要求倒数） 
void fft(cp *a,int n,int flag){ 
    for(int i=0;i<n;i++)
	{
	 //i小于rev[i]时才交换，防止同一个元素交换两次，回到它原来的位置。 
	  if(i<rev[i])swap(a[i],a[rev[i]]);
	}
	for(int h=1;h<n;h*=2)//h是准备合并序列的长度的二分之一
	{
		cp wn=exp(cp(0,flag*pie/h));//求单位根w_n^1 
		for(int j=0;j<n;j+=h*2)//j表示合并到了哪一位
		{
			cp w(1,0);
			for(int k=j;k<j+h;k++)//只扫左半部分，得到右半部分的答案
			{
			    cp x=a[k];
			    cp y=w*a[k+h];
		        a[k]=x+y;  //这两步是蝴蝶变换 
		        a[k+h]=x-y;
		        w*=wn; //求w_n^k 
			}
		 }
	 }
	 //判断是否是FFT还是IFFT 
	 if(flag==-1)
	 for(int i=0;i<n;i++)
     a[i]/=n;
}