您的位置:首页 > 其它

洛谷 P3803 【模板】多项式乘法(FFT)

2018-03-24 21:14 337 查看
算法导论上讲得很好,暂时先不写了。补充一点书上没有的:项数不是2的整数幂时,要在前面添若干项补成2的整数幂,补的各项系数为0。

注意细节:数组不能只开到n+m (2×106)n+m (2×106),因为FFT要求项数必须补到2的整数幂,而大于等于n+mn+m的最小2的整数幂约为26万,所以数组大小要开到这个以上。另外,输出的时候注意前面补的那些项就不用输了(感觉自己好智障)。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
typedef double db;
const int maxn=3e6+10;
const db PI=acos(-1.0);
int n,m,N=1;

inline db read(){
db ret=0,sign=1;
char ch=getchar();
while(ch<'0'||ch>'9') {if(ch=='-') sign=-1; ch=getchar();}
while(ch>='0'&&ch<='9') {ret=ret*10+ch-'0'; ch=getchar();}
return ret*sign;
}

struct complex{
db x,y;
}a[maxn],b[maxn];
complex operator + (complex p,complex q){return (complex){p.x+q.x,p.y+q.y};}
complex operator - (complex p,complex q){return (complex){p.x-q.x,p.y-q.y};}
complex operator * (complex p,complex q){return (complex){p.x*q.x-p.y*q.y,p.x*q.y+p.y*q.x};}
complex operator / (complex p,db q){return (complex){p.x/q,p.y/q};}

int rev[maxn],bit=0;
inline void init(){
while(N<n+m+1) N<<=1,bit++;
for(int i=0;i<N;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}

void fft(complex* p,int len,int dft){
for(int i=0;i<len;i++)
if(i<rev[i]) swap(p[i],p[rev[i]]);
for(int i=1;i<len;i<<=1){
complex wn=(complex){cos(PI/i),dft*sin(PI/i)};
for(int j=0;j<len;j+=i<<1){
complex wnk=(complex){1,0};
for(int k=j;k<j+i;k++){
complex x=p[k],y=wnk*p[k+i];
p[k]=x+y,p[k+i]=x-y;
wnk=wnk*wn;
}
}
}
if(dft==-1) for(int i=0;i<len;i++)
p[i]=p[i]/(db)len;
}

int main(){
n=read(),m=read();
for(int i=0;i<=n;i++)
a[i].x=read();
for(int i=0;i<=m;i++)
b[i].x=read();
init();
fft(a,N,1),fft(b,N,1);
for(int i=0;i<N;i++)
a[i]=a[i]*b[i];
fft(a,N,-1);
for(int i=0;i<=n+m;i++){
if(i!=0) putchar(' ');
printf("%d",(int)(a[i].x+0.5));
}
putchar('\n');
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: