代码

#include<math.h>
using namespace std;
#define rep(i,l,r) for(int i=(l);i<=(r);++i)
#define rpe(i,r,l) for(int i=(r);i>=(l);--i)

typedef double db;
typedef long long ll;
typedef unsigned long long ull;
typedef unsigned uint;

const int N=1048576+10;
template <class T>inline void swap(T &a,T &b){T c=a;a=b;b=c;}
namespace FFT{
    const db pi=acos(-1);
    struct cp{
	db re,im;
	cp(db _re=0,db _im=0){re=_re;im=_im;}
	cp operator +(cp b){return cp(re+b.re,im+b.im);}
	cp operator -(cp b){return cp(re-b.re,im-b.im);}
	cp operator *(cp b){return cp(re*b.re-im*b.im,re*b.im+im*b.re);}
    };
    int r[N<<1];cp c[N<<1];
    cp w[N<<1];
    inline void fft(cp *a,int f,int n){
	rep(i,0,n-1) if(r[i]>i) swap(a[r[i]],a[i]);
	for(int i=1;i<n;i<<=1){
	    cp wn(cos(pi/i),f*sin(pi/i));
	    w[0]=cp(1,0);rep(k,1,i-1) w[i]=w[i]*wn;
	    for(int j=0,p=(i<<1);j<n;j+=p){
		for(int k=0;k<i;++k){
		    cp x=a[j+k],y=w[k]*a[j+k+i];
		    a[j+k]=x+y;a[j+k+i]=x-y;
		}
	    }
	}
	if(f==-1){rep(i,0,n-1) a[i].re/=n,a[i].im/=n;}
    }
    inline void mul(int n,int m,uint *a,uint *b,uint *cc){
	n+=m;rep(i,0,n) c[i]=cp(a[i],b[i]);
	int l=0;m=n;for(n=1;n<=m;n<<=1) ++l;
	rep(i,0,n-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	rep(i,m+1,n) c[i]=cp(0,0);
	fft(c,1,n);rep(i,0,n-1) c[i]=c[i]*c[i];
	fft(c,-1,n);
	rep(i,0,m) cc[i]=c[i].im/2;
    }
}
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c){
    FFT::mul(n,m,a,b,c);
}

评测结果