代码

#include<bits/stdc++.h>
#pragma GCC optimize("Ofast")
#define _SIZE_ 100000
char _b[_SIZE_],*_b1,*_b2;
#define getc() (_b1==_b2?fread(_b,1,_SIZE_,stdin),_b2=_b+_SIZE_,*((_b1=_b)++):*(_b1++))
#define read(__x) {register signed _x(0);register char _c=getc(),_f(0);for(;_c<47;_c=getc())_f=(_c==45);for(;_c>47;_c=getc())_x=_x*10+(_c^48);_f&&(_x=-_x),__x=_x;}
char _d[_SIZE_],*_p=_d;
#define putc(__c) (_p-_d==_SIZE_?fwrite(_d,1,_SIZE_,stdout),_p=_d,*_p++=__c:*_p++=__c)
#define write(__x) {register signed _x=__x;(_x<0)&&(putc(45),_x=-_x);static signed _q[11];register char _t=0;do{_q[_t++]=_x%10,_x/=10;}while(_x);while(_t)putc(_q[--_t]+48);}
#define fwflush() (fwrite(_d,1,_p-_d,stdout),_p=_d)
#define pi 3.14159265358979323846
struct C
{
    double x,y;
    C operator+(const C b)const{return {x+b.x,y+b.y};} 
	C operator-(const C b)const{return {x-b.x,y-b.y};}
	C operator*(const C b)const{return {x*b.x-y*b.y,x*b.y+y*b.x};}
};
const int _lg2[]={0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3};
using namespace std;
int tmp1,n,m,N;
int lg2(int i)
{
	int tmp(0);
	if(i>65535)tmp|=16,i>>=16;
	if(i>255)tmp|=8,i>>=8;
	if(i>15)tmp|=4,i>>=4;
	return tmp+_lg2[i];
}
double co[23],si[23];
void fft(C*a,int n)
{
	if(n==1)return;
	int mid=n>>1;
	C x__={1,0},x_={co[lg2(n)],si[lg2(n)]};
	fft(a,mid),fft(a+mid,mid);
	for(int i=0;i<mid;i++,x__=x__*x_)
	{
		C x(a[i+mid]*(C){x__.x,x__.y});
		a[i+mid]=a[i]-x,a[i]=a[i]+x;
	}
}
void ifft(C*a,int n)
{
	if(n==1)return;
	int mid=n>>1;
	C x__={1,0},x_={co[lg2(n)],si[lg2(n)]};
	ifft(a,mid),ifft(a+mid,mid);
	for(int i=0;i<mid;i++,x__=x__*x_)
	{
		C x(a[i+mid]*(C){x__.x,-x__.y});
		a[i+mid]=a[i]-x,a[i]=a[i]+x;
	}
}
C a[1<<21];
int r[1<<21];
void rev(C*a,int n){for(int i=0;i<n;i++)if(i<r[i])swap(a[i],a[r[i]]);}
void poly_multiply(unsigned *_a, int n, unsigned *_b, int m, unsigned *_c)
{
	for(int i=0;i<=22;i++)co[i]=cos(2*pi/(1<<i)),si[i]=sin(2*pi/(1<<i));
	int l_n(0);
	++n,m=n,N=1<<21;
	for(int i=0;i<N;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l_n-1));
	for(int i=0;i<n;i++)a[i].x=_a[i],a[i].y=_b[i];
	rev(a,N);
	fft(a,N);
	for(int i=0;i<N;i++)a[i]=a[i]*a[i];
	rev(a,N);
	ifft(a,N);
	for(int i=0;i<n+m-1;i++){_c[i]=round(a[i].y/N/2);}
	fwflush();
	return 0;
}

评测结果