FFT(NTT)学习笔记
#include<bits/stdc++.h>
using namespace std;
int n1,n2,P=1,n,jw,tot,res[1000010];
char ip1[100010],ip2[100010];
double Pi=3.14159265;
struct comp
{
comp (double xx=0,double yy=0)
{
x=xx,y=yy;
}
double x,y;
};
comp operator + (comp a,comp b)
{
return comp(a.x+b.x,a.y+b.y);
}
comp operator - (comp a,comp b)
{
return comp(a.x-b.x,a.y-b.y);
}
comp operator * (comp a,comp b)
{
return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
comp a[500010],b[500010],tmp[500010];
int FFT(comp *s,int now,int len,int pd)
{
if(!len)
return 0;
for(int i=now;i<now+len*2;i++)
tmp[i]=s[i];
for (int i=now;i<now+len;i++)
{
s[i]=tmp[now+((i-now)*2)];
s[i+len]=tmp[now+(((i-now)*2)|1)];
}
FFT(s,now,len/2,pd);
FFT(s,now+len,len/2,pd);
comp dwg,nowdwg;
dwg=comp(cos(Pi/len),sin(Pi/len)*pd);
nowdwg.x=1;
nowdwg.y=0;
for(int i=now;i<now+len;i++)
{
comp tt=nowdwg*s[i+len];
tmp[i]=s[i]+tt;
tmp[i+len]=s[i]-tt;
nowdwg=nowdwg*dwg;
}
for(int i=now;i<now+len*2;i++)
s[i]=tmp[i];
}
int TS(double x)
{
if(x-floor(x)>=0.5)
{
return ceil(x);
}
else
{
return floor(x);
}
}
int main()
{
scanf("%d",&n);
n--;
n1=n;n2=n;
cin>>ip1;
for(int i=0;i<=n1;i++)
{
a[i].x=ip1[i]-48;
}
cin>>ip2;
for(int i=0;i<=n2;i++)
{
b[i].x=ip2[i]-48;
}
for(;P<=n1+n2;P*=2)
{
}
FFT(a,0,P/2,1);
FFT(b,0,P/2,1);
for(int i=0;i<P;i++)
a[i]=a[i]*b[i];
FFT(a,0,P/2,-1);
/*for(int i=0;i<=n1+n2;i++)
printf("%.0lf ",fabs(a[i].x)/P);*/
for(int i=n1+n2;i>=0;i--)
{
res[tot++]=TS(fabs(a[i].x)/P);
}
for(int i=0;i<=P;i++)
{
res[i+1]+=res[i]/10;
res[i]%=10;
}
while(!res[P])
{
P--;
}
for(int i=P;i>=0;i--)
{
printf("%d",res[i]);
}
}
FFT板子-2
AC代码
#include<bits/stdc++.h>
using namespace std;
int n1,n2,P=1;
double Pi=3.14159265;
struct comp
{
comp (double xx=0,double yy=0)
{
x=xx,y=yy;
}
double x,y;
};
comp operator + (comp a,comp b)
{
return comp(a.x+b.x,a.y+b.y);
}
comp operator - (comp a,comp b)
{
return comp(a.x-b.x,a.y-b.y);
}
comp operator * (comp a,comp b)
{
return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
comp a[500010],b[500010],tmp[500010];
int FFT(comp *s,int now,int len,int pd)
{
if(!len)
return 0;
for(int i=now;i<now+len*2;i++)
tmp[i]=s[i];
for (int i=now;i<now+len;i++)
{
s[i]=tmp[now+((i-now)*2)];
s[i+len]=tmp[now+(((i-now)*2)|1)];
}
FFT(s,now,len/2,pd);
FFT(s,now+len,len/2,pd);
comp dwg,nowdwg;
dwg=comp(cos(Pi/len),sin(Pi/len)*pd);
nowdwg.x=1;
nowdwg.y=0;
for(int i=now;i<now+len;i++)
{
comp tt=nowdwg*s[i+len];
tmp[i]=s[i]+tt;
tmp[i+len]=s[i]-tt;
nowdwg=nowdwg*dwg;
}
for(int i=now;i<now+len*2;i++)
s[i]=tmp[i];
}
int main()
{
scanf("%d%d",&n1,&n2);
for(int i=0;i<=n1;i++)
{
scanf("%lf",&a[i].x);
}
for(int i=0;i<=n2;i++)
{
scanf("%lf",&b[i].x);
}
for(;P<=n1+n2;P*=2)
{
}
FFT(a,0,P/2,1);
FFT(b,0,P/2,1);
for(int i=0;i<P;i++)
a[i]=a[i]*b[i];
FFT(a,0,P/2,-1);
for(int i=0;i<=n1+n2;i++)
printf("%.0lf ",fabs(a[i].x)/P);
}
然后是蝴蝶变换版(更好打?)
#include<bits/stdc++.h>
using namespace std;
int n1,n2,P=1,r[2000010];
double Pi=3.14159265;
struct comp
{
comp (double xx=0,double yy=0)
{
x=xx,y=yy;
}
double x,y;
};
comp operator + (comp a,comp b)
{
return comp(a.x+b.x,a.y+b.y);
}
comp operator - (comp a,comp b)
{
return comp(a.x-b.x,a.y-b.y);
}
comp operator * (comp a,comp b)
{
return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
comp a[500010],b[500010],tmp[500010];
int FFT(comp *s,int pd)
{
for(int i=0;i<P;i++)
{
if(i<r[i])
{
comp tmpp=s[i];
s[i]=s[r[i]];
s[r[i]]=tmpp;
}
}
for(int len=2;len<=P;len*=2)
{
int Len=len/2;
comp dwg(cos(Pi/Len),pd*sin(Pi/Len));
for(int k=0;k<P;k+=len)
{
comp qwq(1,0);
for(int l=k;l<k+Len;l++)
{
comp tt=qwq*s[Len+l];
s[Len+l]=s[l]-tt;
s[l]=s[l]+tt;
qwq=qwq*dwg;
}
}
}
}
int main()
{
scanf("%d%d",&n1,&n2);
for(int i=0;i<=n1;i++)
{
scanf("%lf",&a[i].x);
}
for(int i=0;i<=n2;i++)
{
scanf("%lf",&b[i].x);
}
for(;P<=n1+n2;P*=2)
{
}
for(int i=0;i<P;i++)
r[i]=(r[i/2]/2)+((i%2)?P>>1:0);
FFT(a,1);
FFT(b,1);
for(int i=0;i<P;i++)
{
a[i]=a[i]*b[i];
}
FFT(a,-1);
for(int i=0;i<=n1+n2;i++)
printf("%.0f ",fabs(a[i].x)/P);
}
NTT？
咕咕咕
