想法:
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采用复数的概念,记表示为把单位圆角度分成n份逆时针第一份角度和圆的交点的位置即(cosθ,i*sinθ)。(精髓所在周期性,用单位元复数即没有值的改变,只有方向的改变)
我们要知道复数运算规则,相加相当于向量相加,相乘相当于角度相加模相乘(至于为什么可以分开来乘算一下),相除a/b相当于把a的角度-b的角度再除以b模的平方。即
思路:
我们要找、
......
的值,我们发现由于替换成单位元复数形式由于周期性
已经有了联系,即不是所有的多项式幂次方都变化,k变化的时候有些由于周期性没有变化且
没有重复,我们想用递推的方式求
,把k看成二进制,任意k可以由不同
、
、
构成(假设n为2的幂次方),假设k=
+
+
,
=
,相当于一开始有n条向量在x正半轴的一个圆,你要拨动不同的向量让他指向正确的位置,而这个拨动不能是
,而是要在多项式内部用单位圆复数乘法来拨动。因此很多如何将向量分类、打包然后一起拨动是关键。
那么我们观察内部多项式受影响的项,受影响的、受
影响的........受
影响的进行分类,然后拨与不拨衍生出不同的
,我们发现假设n=8,k=1时要拨动的有8个,间隔为1,k=2时要拨动的有4个(因为有些向量转了一圈又和其他的重叠了,我们把这些打包在一起一起拨动),间隔为4,k=4时要拨动2个包,间隔为2。
第四层每组数都有,第三层每组数都有
和
,第二层每组数有
、
、
、
。
这样我们就能在nlogn的复杂度求出所有的。
那么考虑由点求系数,即IDFT。
即
脑补代入的过程,可得
那么结论就是
证明:
右边 =
代入可得:
分类讨论:
① j==kj==k
那么贡献就是
② j≠kj≠k ,设 p=j-kp=j−k
那么贡献就是等比数列求和可以得到:
所以这部分贡献为0,我们就成功地证明了上述结论。
引用:https://www.luogu.com.cn/blog/command-block/fft-xue-xi-bi-ji
其实证明也没那么复杂,我们相当于,这看起来十分自私只考虑把这个k次幂的向量给拨动回去,实际上其他的p次幂的向量p-k!=0,相当于i=0时向量在0,i=1时在
,i=2时在
,所以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;
}