代码

#include <queue>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define MP make_pair
#define ll long long
#define fi first
#define se second
using namespace std;

int read(void) {
    int x = 0; bool f = 0;
    char c = getchar();
    for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
    for (;isdigit(c);c=getchar()) x=x*10+(c^48);
    if (f) x=-x; return x;
}

inline void write(int x, char ed = '\n')
{
	static short st[30];short tp=0;
	if(x<0) putchar('-'),x=-x;
	do st[++tp]=x%10,x/=10; while(x);
	while(tp) putchar('0'|st[tp--]);
	putchar(ed);
}

#define op com operator 
#define con const
typedef double db;
const int N = 3000005;
const double Pi = acos(-1.0);
struct com {
	db x, y;
	com(db a = 0, db b = 0) : x(a) , y(b) {}
	op + (con com &w) con { return com(x + w.x, y + w.y); }
	op - (con com &w) con { return com(x - w.x, y - w.y); }
	op * (con com &w) con { return com(x * w.x - y * w.y, x * w.y + y * w.x); }
	op - (void) con { return com(-x, -y); }
	com mi() con { return com(-y, x); }
}A[N], E[N];

int r[N], lim = 1, L; 
void dft(com *A) {
	for (int i = 1;i < lim; i++)
		if (r[i] > i) swap(A[i], A[r[i]]);
	if (lim >= 2) 
	for (int j = 0;j < lim; j += 2) {
		com x = A[j], y = A[j+1];
		A[j] = x + y, A[j+1] = x - y;
	}
	if (lim >= 4)
	for (int j = 0;j < lim; j += 4) {
		com x = A[j], y = A[j+2];
		A[j] = x + y, A[j+2] = x - y;
		x = A[j+1], y = A[j+3].mi();
		A[j+1] = x + y, A[j+3] = x - y;
	}
	if (lim >= 8)
	for (int j = 0;j < lim; j += 8) {
		com x = A[j], y = A[j+4];
		A[j] = x + y, A[j+4] = x - y;
		x = A[j+1], y = A[j+5] * E[5];
		A[j+1] = x + y, A[j+5] = x - y;
		x = A[j+2], y = A[j+6].mi();
		A[j+2] = x + y, A[j+6] = x - y;
		x = A[j+3], y = A[j+7] * E[7];
		A[j+3] = x + y, A[j+7] = x - y;
	}
	for (int i = 8;i < lim; i <<= 1) {
		for (int j = 0;j < lim; j += (i << 1)) {
			com *f = A + j, *g = f + i, *e = E + i;
			for (int k = 0;k < i; k++) {
				com x = f[k], y = g[k] * e[k];
				f[k] = x + y, g[k] = x - y;
				k++;
				x = f[k], y = g[k] * e[k];
				f[k] = x + y, g[k] = x - y;
			}
		}
	}
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
	for (int i = 0;i <= n; i++) A[i].x = a[i];
	for (int i = 0;i <= m; i++) A[i].y = b[i];
	while (lim <= (n + m)) lim <<= 1, L++;
	int len = lim >> 1;
	for (int i = 1;i < lim; i++)
		r[i] = r[i >> 1] >> 1 | ((i & 1) ? len : 0);
	E[1] = com(1, 0);
	for (int i = 2;i < lim; i <<= 1) {
		com *e0 = E + i / 2, *e1 = E + i;
		com 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;
	}
	dft(A);
	for (int i = 0;i < lim; i++) A[i] = A[i] * A[i];
	dft(A); reverse(A + 1, A + lim); lim *= 2;
	for (int i = 0;i <= n + m; i++)
		c[i] = (unsigned int)(A[i].y / lim + 0.5);
}

评测结果