代码

#include <math.h>

#define N 1048576
#define pi 3.1415926535897932384626
#define op cplx operator
typedef double fp;
struct cplx{
	fp c,s;
	cplx(fp _c=0,fp _s=0):c(_c),s(_s){}
	op+(const cplx&w)const{return cplx(c+w.c,s+w.s);}
	op-(const cplx&w)const{return cplx(c-w.c,s-w.s);}
	op*(const cplx&w)const{return cplx(c*w.c-s*w.s,s*w.c+c*w.s);}
	op-()const{return cplx(-c,-s);}
	cplx mi()const{return cplx(-s,c);}
}A[N],B[N],C[N],E[N];
typedef cplx*const cc;
int rev[N];
__attribute((opitmize("-o3")))
void dft(cc X){
	for(int j=0;j<N;j+=2){
		cplx x=X[j],y=X[j+1];
		X[j]=x+y;
		X[j+1]=x-y;
	}
	for(int i=2;i<N;i<<=1){
		for(int j=0;j<N;j+=i<<1){
			cc f=X+j,g=f+i,e=E+i;
			for(int k=0;k<i;++k){
				cplx x=f[k],y=g[k]*e[k];
				f[k]=x+y;
				g[k]=x-y;
			}
		}
	}
}
__attribute((opitmize("-o3")))
void cal(unsigned*a,int n,cc A){
	for(int i=0;i<=n;i+=2)A[rev[i>>1]]=cplx(a[i],i==n?0:a[i+1]);
	dft(A);
	A[N]=A[0];
}
__attribute((opitmize("-o3")))
void trans(cplx&a0,cplx&a1,const cplx&v0,const cplx&v1){
	a0=cplx(v0.c+v1.c,v0.s-v1.s);
	a1=cplx(-v0.s-v1.s,v0.c-v1.c);
}
__attribute((opitmize("-o3")))
void poly_multiply(unsigned*a,int n,unsigned*b,int m,unsigned*c){
	for(int i=1;i<N;++i)rev[i]=rev[i>>1]>>1|(i&1)<<19;
	E[1]=cplx(1,0);
	for(int i=2;i<N;i<<=1){
		cc e0=E+i/2,e1=E+i;
		cplx w(cos(pi/i),sin(pi/i));
		for(int j=0;j<i;j+=2)e1[j]=e0[j>>1],e1[j+1]=e1[j]*w;
	}
	cal(a,n,A);
	cal(b,m,B);
	for(int i=0;i<N;++i){
		cplx A0,A1,B0,B1;
		trans(A0,A1,A[i],A[N-i]);
		trans(B0,B1,B[i],B[N-i]);
		cplx e=i<N/2?E[N/2+i]:-E[i];
		C[rev[i]]=A0*B0+A1*B1*e+(A0*B1+A1*B0).mi();
	}
	for(int i=0;i<N;++i)C[i].s*=-1;
	dft(C);
	for(int i=0;i<=n+m;i+=2){
		cplx x=C[i>>1]-C[N-1];
		c[i]=int(x.c/(N*4)+.5);
		if(i==n+m)break;
		c[i+1]=int(x.s/(N*4)+.5);
	}
}

评测结果