代码

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define per(i,a,b) for(int i=(a);i>=(b);i--)
#define REP(i,n) for(int i=(0);i<(n);i++)
#define fi first
#define se second
#define pb push_back
#define mp make_pair
using namespace std;
typedef pair<int,int> pii;
typedef vector<int> vi;
typedef long long ll;
typedef double ld;
 
const ld pi=acos(-1),inf=1e99;
struct C{
	ld r,i;
	C operator+(const C &o)const{return (C){r+o.r,i+o.i};}
	C operator-(const C &o)const{return (C){r-o.r,i-o.i};}
	C operator/(const ld&o)const{return (C){r/o,i/o};}
	C operator*(const C &o)const{return (C){r*o.r-i*o.i,r*o.i+i*o.r};}
	C mulR(const ld&o)const{return (C){r*o,i*o};}
	C mulI(const ld&o)const{return (C){-i*o,r*o};};
	C conj()const{return (C){r,-i};}
};
typedef vector<ld> po;
 
namespace Poly{
	const int N=(1<<20)+5; int mx;
	void dft(C *a,int n){
		static C W[N<<1],*H[30],*las=W;
		for(;mx<n;mx++){
			H[mx]=las; C w=(C){1,0};
			C wn=(C){cos(pi/(1<<mx)),sin(pi/(1<<mx))};
			REP(i,1<<mx)*las++=w,w=w*wn;
		}
		for(int i=0,j=0;i<(1<<n);i++){
			if(j<i)swap(a[i],a[j]);
			for(int k=1<<(n-1);(j^=k)<k;k>>=1);
		}
		int m=1<<n;
		for(int k=0,d=1;k<n;k++,d<<=1)
			for(int i=0;i<m;i+=d<<1){
				C *l=a+i,*r=a+i+d,*w=H[k];
				for(int j=0;j<d;++j,++l,++r){
					const C t=(*w++)*(*r);
					*r=*l-t,*l=*l+t;
				}
			}
	}
	void Ddft(C *a,C *b,int n){
		int m=1<<n;
		REP(i,m)a[i].i=b[i].r;
		dft(a,n);
		REP(i,m)b[i]=a[i?m-i:0].conj();
		REP(i,m){
			C t=(a[i]+b[i]).mulR(0.5);
			b[i]=(a[i]-b[i]).mulI(-0.5);
			a[i]=t;
		}
	}
	void idft(C *a,int n){
		reverse(a+1,a+(1<<n));
		dft(a,n); ld pw=1<<n;
		REP(i,1<<n)a[i]=a[i]/pw;
	}
	void Idft(C *a,int n){
		static C b[N]; int m=1<<(n-1); C w=(C){1,0};
		C wn=(C){cos(pi/(1<<(n-1))),-sin(pi/(1<<(n-1)))};
		REP(i,m){
			b[i]=(a[i]+a[m+i]).mulR(0.5)
				+(a[i]-a[m+i]).mulI(0.5)*w;
			w=w*wn;
		}
		idft(b,n-1);
		REP(i,m){
			a[i<<1]=(C){b[i].r,0};
			a[i<<1|1]=(C){b[i].i,0};
		}
	}
	po mul(po &A,po &B){
		if(min(A.size(),B.size())<=16){
			po res(A.size()+B.size());
			REP(i,A.size())REP(j,B.size())
				res[i+j]+=A[i]*B[j];
			return res;
		}
		static C a[N],b[N];
		int aim=A.size()+B.size()-1,n=1;
		while((1<<n)<aim)n++;
		REP(i,1<<n)a[i]=b[i]=(C){0,0};
		REP(i,A.size())a[i].r=A[i];
		REP(i,B.size())b[i].r=B[i];
		Ddft(a,b,n);
		REP(i,1<<n)a[i]=a[i]*b[i];
		Idft(a,n); po res(aim);
		REP(i,aim)res[i]=a[i].r;
		return res;
	}
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c){
n++,m++;
po A(n,0),B(m,0);
REP(i,n) A[i]=a[i];
REP(i,m) B[i]=b[i];
po C=Poly::mul(A,B);
REP(i,n+m-1) c[i]=C[i];
}

评测结果