代码

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>

using std::swap;
using std::reverse;

const double eps(0.4), pi(acos(-1.0));

inline int toInt(double x) {return int(x+(x>0? eps: -eps));}

struct C
{
    double r, i;
    C() {} C(double r): r(r), i(0) {}
    C(double r, double i): r(r), i(i) {}

    inline C operator +(const C &c) const {return C(r+c.r, i+c.i);}
    inline C operator -(const C &c) const {return C(r-c.r, i-c.i);}
    inline C operator *(const C &c) const {return C(r*c.r-i*c.i, r*c.i+i*c.r);}
    inline C conj() {return C(r, -i);}
};

inline void dFT(C* a, int s)
{
    static C b[1<<21|233];
    C *f(a), *g(b);
    for(int k(s>>1); k; k>>=1, swap(f, g))
        for(int i(0); i<s>>1; i+=k)
        {
            C *g0(g+i), *g1(g+(s>>1|i)), *f0(f+(i<<1)), *f1(f+(i<<1|k)), ome(cos(2*pi*i/s), sin(2*pi*i/s));
            for(int j(0); j<k; ++j)
            {
                C h(f1[j]*ome);
                g0[j] = f0[j] + h;
                g1[j] = f0[j] - h;
            }
        }
    for(int i(0); i<s; ++i)
        a[i] = f[i];
}

inline void iDFT(C* a, int s)
{
    dFT(a, s);
    reverse(a+1, a+s);
    double k(1.0/s);
    for(int i(0); i<s; ++i)
        a[i].r*=k, a[i].i*=k;
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *r)
{
    int s(1);
    for(; s<n+m+1; s<<=1);

    static C c[1<<21|233];
    for(int i(0); i<s; ++i)
        c[i] = C(i<=n? a[i]: 0, i<=m? b[i]: 0);
    dFT(c, s);
    c[0] = c[0].r * c[0].i;
    c[s>>1] = c[s>>1].r * c[s>>1].i;
    for(int i(1); i<s>>1; ++i)
    {
        C x(c[i]), y(c[s-i]);
        c[i] = (x+y.conj())*(x-y.conj())*C(0, -0.25);
        c[s-i] = (y+x.conj())*(y-x.conj())*C(0, -0.25);
    }
    iDFT(c, s);
    for(int i(0); i<=n+m; ++i)
        r[i] = toInt(c[i].r);
}

评测结果