FFT

最新推荐文章于 2025-02-04 16:17:56 发布

struct_liao

最新推荐文章于 2025-02-04 16:17:56 发布

阅读量232

点赞数

本文链接：https://blog.csdn.net/struct_liao/article/details/104525742

版权

学习笔记专栏收录该内容

1 篇文章

订阅专栏

想法：

FFT用来解决多项式相乘，A(x)=(a1*x^1+a2*x2^2.......+an*x^n)比如两个1000000位的数相乘可以看作两个n=1000000项的x=10的多项式相乘。

我们用点来代表多项式，n个点能解出n项多项式方程的系数，比如面两个1e6项多项式称出来我们就有2e6项，那么我们两边多项式都弄2e6个点然后y相乘就行了。

然鹅这样无论是造点还是解系数方程复杂度都是n^2的。

我们想能不能带个特殊的x进去让不同的点之间找到多项式内部可评估的联系（即每次改变x只有部分多项式改变），这样产生点就不是n^2复杂度。

运算规则：

DFT应运而生，DFT采用复数的概念，记 $\omega _{n}^{1}$ 表示为把单位圆角度分成n份逆时针第一份角度和圆的交点的位置即（cosθ，i*sinθ）。（精髓所在周期性，用单位元复数即没有值的改变，只有方向的改变）

我们要知道复数运算规则，相加相当于向量相加，相乘相当于角度相加模相乘（至于为什么可以分开来乘算一下），相除a/b相当于把a的角度-b的角度再除以b模的平方。即

思路：

我们要找 $A(\omega _{n}^{0})$ 、 $A(\omega _{n}^{1})$ ...... $A(\omega _{n}^{n-1})$ 的值，我们发现由于替换成单位元复数形式由于周期性 $A(\omega _{k}^{n})$ 已经有了联系，即不是所有的多项式幂次方都变化，k变化的时候有些由于周期性没有变化且 $A(\omega _{n}^{k})$ 没有重复，我们想用递推的方式求 $A(\omega _{k}^{n})$ ，把k看成二进制，任意k可以由不同 $\frac{n}{2}$ 、 $\frac{n}{4}$ 、 $\frac{n}{n}$ 构成（假设n为2的幂次方），假设k= $\frac{n}{2}$ + $\frac{n}{4}$ + $\frac{n}{n}$ ， $A(\omega _{n}^{k})$ = $A(\omega _{n}^{\frac{n}{2}+\frac{n}{4}+\frac{n}{n}})$ ，相当于一开始有n条向量在x正半轴的一个圆，你要拨动不同的向量让他指向正确的位置，而这个拨动不能是 $A(\omega _{n}^{\frac{n}{2}+\frac{n}{4}+\frac{n}{n}})=A(\omega _{n}^{\frac{n}{2}})+A(\omega _{n}^{\frac{n}{4}})+A(\omega _{n}^{\frac{n}{n}})$ ，而是要在多项式内部用单位圆复数乘法来拨动。因此很多如何将向量分类、打包然后一起拨动是关键。

那么我们观察内部多项式受影响的项，受 $\frac{n}{2}$ 影响的、受 $\frac{n}{4}$ 影响的........受 $\frac{n}{n}$ 影响的进行分类，然后拨与不拨衍生出不同的 $A(\omega _{n}^{k})$ ，我们发现假设n=8，k=1时要拨动的有8个，间隔为1，k=2时要拨动的有4个（因为有些向量转了一圈又和其他的重叠了，我们把这些打包在一起一起拨动），间隔为4，k=4时要拨动2个包，间隔为2。

第四层每组数都有 $A(\omega _{n}^{0})$ ，第三层每组数都有 $A(\omega _{n}^{0})$ 和 $A(\omega _{n}^{1})$ ，第二层每组数有 $A(\omega _{n}^{0})$ 、 $A(\omega _{n}^{1})$ 、 $A(\omega _{n}^{2})$ 、 $A(\omega _{n}^{3})$ 。

这样我们就能在nlogn的复杂度求出所有的 $A(\omega _{n}^{k})$ 。

那么考虑由点求系数，即IDFT。

即 $G=DFT(F)$

脑补代入的过程,可得 $G[k]=\sum_{i=0}^{n-1}(\omega _{n}^{k})^{i}F[i]$

那么结论就是 $n*F[k]=\sum_{i=0}^{n-1}(\omega _{n}^{-k})G[i]$

证明:

右边 = $\sum_{i=0}^{n-1}(\omega _{n}^{-k})^iG[i]$

代入可得:

$=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}\omega^{ij}_n\omega^{-ik}_nF[j]$

$=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}\omega^{i(j-k)}_nF[j]$

分类讨论:

① j==kj==k
那么贡献就是 $\sum\limits_{i=0}^{n-1}\omega^{0}_nF[k]=nF[k]$
② j≠kj≠k ,设 p=j-kp=j−k
那么贡献就是 $\sum\limits_{i=0}^{n-1}\omega^{ip}_nF[k+p]$
$=\omega^{p}_n(\sum\limits_{i=0}^{n-1}\omega^{i}_n)F[k+p]$

等比数列求和可以得到:

$(\sum\limits_{i=0}^{n-1}\omega^{i}_n)=\dfrac{\omega_n^{np}-1}{\omega_n^{p}-1}=\dfrac{\omega_n^{0}-1}{\omega_n^{p}-1}$

所以这部分贡献为0,我们就成功地证明了上述结论。

引用：https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji

其实证明也没那么复杂，我们相当于 $n*F[k]=\sum_{i=0}^{n-1}(\omega _{n}^{-k})G[i]$ ，这看起来十分自私只考虑把这个k次幂的向量给拨动回去，实际上其他的p次幂的向量p-k!=0，相当于i=0时向量在0，i=1时在 $\omega _{n}^{p-k}$ ，i=2时在 $\omega _{n}^{2*(p-k)}$ ，所以0~n-1走完后这个向量在图的位置是两两关于原点对称的（可以根据模运算证明），相加等于0。而我们要求的那个次方幂因为次次都在X正半轴不会走，所以不会对称消除。

模版题：P1919 【模板】A*B Problem升级版（FFT快速傅里叶）

代码：

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
using namespace std;
double pai=acos(-1.0);
double co[2400010],si[2400010];
int ans[2400010];
int a[2400010],b[2400010];
struct complex{
	double real,imag;
	complex(double x=0,double y=0){real=x;imag=y;}
	complex operator +(complex &b){return complex(real+b.real,imag+b.imag);}
	complex operator -(complex &b){return complex(real-b.real,imag-b.imag);}
	complex operator *(complex &b){return complex(real*b.real-imag*b.imag,real*b.imag+imag*b.real);}
}tmp[2400010],ret[2400010],A[2400010]; 

void fft(int n,int flag){
	int i;
	for(i=1;i<n;){//i代表这一步要执行w(i/n)，也代表每个系数携带几个w答案即w0~w(0+i-1)
		int gap=(n/i)/2;//下一个i时多少步一个循环
		for(int d=0;d<=gap-1;d++){//第几部分系数集合
			int v=d+gap,now=d*i,now1=v*i,st=2*d*i,st1=(2*d+1)*i;
			complex op(co[d*i],flag*si[d*i]);
			for(int p=0;p<i;p++){
				tmp[p+st]=ret[p+now]+ret[p+now1];
				tmp[p+st1]=(ret[p+now]-ret[p+now1])*op;//因为这两系数对应的幂*op后刚好关于原点对称，而且下一部要合并，所以为了减少乘法运算就先减再乘
			}
		}
		i<<=1; 
		for(int d=0;d<=gap-1;d++)for(int p=0;p<i;p++)ret[p+d*i]=tmp[p+d*i];
		
	}
}
int read1(int a[]){
	char p=0;
	while(p<'0' || p>'9')p=getchar();
	for(int i=0;;i++){
		a[i]=p-'0';
		p=getchar();
		if(p<'0' || p>'9')return i+1;
	}
}
void home(int n){
	for(int i=0;i<=n;i++)co[i]=cos(i*2*pai/n),si[i]=sin(i*2*pai/n);
	return;
}
int main(){
	freopen("P1919_10.in","r",stdin);
	freopen("P1919_10.out","w",stdout);
	int i,n,n1,n2,l,r;
	n1=read1(a);
	n2=read1(b);
	n=1;
	while(n<(n1+n2))n<<=1;
	home(n); 
	for(i=0;i<=n-1;i++)ret[i]=a[n1-i-1];
	fft(n,1);
	for(i=0;i<n;i++)A[i]=ret[i];
	
	for(i=0;i<=n-1;i++)ret[i]=b[n2-i-1];
	fft(n,1);
	for(i=0;i<n;i++)A[i]=A[i]*ret[i];
	
	for(i=0;i<=n-1;i++)ret[i]=A[i];
	fft(n,-1);
	for(i=0;i<n;i++){
		ans[n-i]+=ret[i].real/n+0.4;
		if(ans[n-i]>=10){
			ans[n-i-1]+=ans[n-i]/10;
			ans[n-i]%=10;
		}
	}
	i=0;
	while(!ans[i] && i<=n)i++;
	for(;i<=n;i++)printf("%d",ans[i]);
	return 0;
}