代码

#include <algorithm>
#include <math.h>
typedef long long ll;
using namespace std;
const int N = 4e6 + 5;
namespace { // radix_4_ntt
typedef long long valueType;
const valueType P = 1004535809, G = 3;
valueType qpow(valueType x, valueType y) {
    valueType res = 1;
    for (; y; y >>= 1, x = x * x % P) {
        if (y & 1) res = res * x % P;
    }
    return res;
}
valueType inv(valueType x) { return qpow(x, P - 2); }
const valueType IMAG = P - qpow(G, P - 1 >> 2);
int getlog(int n) {
    int lgn = 0;
    while (n >>= 1) ++lgn;
    return lgn;
}
void revbin(int n, valueType x[]) {
    for (int i = 0, j = 0; i < n; ++i) {
        if (i < j) std::swap(x[i], x[j]);
        for (int k = n >> 1; (j ^= k) < k; k >>= 1) {
        }
    }
}
void dif_fft_core(int n, valueType x[]) {
    int lgn = getlog(n);
    for (int i = n; i >= 4 << (lgn & 1); i >>= 2) {
        valueType w = qpow(G, (P - 1) / i);
        for (int j = 0, l = i >> 2; j < n; j += i) {
            valueType o1 = 1, o2 = 1, o3 = 1;
            for (int k = 0; k < l; ++k) {
                valueType A = x[j + k], B = x[j + k + l], C = x[j + k + 2 * l],
                          D = x[j + k + 3 * l], IB = B * IMAG % P, ID = D * IMAG % P;
                x[j + k] = (A + B + C + D) % P;
                x[j + k + l] = o2 * (A - B + C - D) % P;
                x[j + k + 2 * l] = o1 * (A - IB - C + ID) % P;
                x[j + k + 3 * l] = o3 * (A + IB - C - ID) % P;
                o1 = o1 * w % P, o2 = o1 * o1 % P, o3 = o1 * o2 % P;
            }
        }
    }
    if (lgn & 1) {
        for (int j = 0; j < n; j += 2) {
            valueType A = x[j], B = x[j + 1];
            x[j] = (A + B) % P, x[j + 1] = (A - B) % P;
        }
    }
}
void dit_fft_core(int n, valueType x[]) {
    int lgn = getlog(n);
    if (lgn & 1) {
        for (int j = 0; j < n; j += 2) {
            valueType A = x[j], B = x[j + 1];
            x[j] = (A + B) % P, x[j + 1] = (A - B) % P;
        }
    }
    for (int i = 4 << (lgn & 1); i <= n; i <<= 2) {
        valueType w = qpow(G, (P - 1) / i);
        for (int j = 0, l = i >> 2; j < n; j += i) {
            valueType o1 = 1, o2 = 1, o3 = 1;
            for (int k = 0; k < l; ++k) {
                valueType A = x[j + k], C = o2 * x[j + k + l] % P, B = o1 * x[j + k + 2 * l] % P,
                          D = o3 * x[j + k + 3 * l] % P, IB = B * IMAG % P, ID = D * IMAG % P;
                x[j + k] = (A + B + C + D) % P;
                x[j + k + l] = (A - IB - C + ID) % P;
                x[j + k + 2 * l] = (A - B + C - D) % P;
                x[j + k + 3 * l] = (A + IB - C - ID) % P;
                o1 = o1 * w % P, o2 = o1 * o1 % P, o3 = o1 * o2 % P;
            }
        }
    }
}
void dft(int n, valueType x[]) { dif_fft_core(n, x); }
void idft(int n, valueType x[]) {
    dit_fft_core(n, x), std::reverse(x + 1, x + n);
    valueType in = P - (P - 1) / n;
    for (int i = 0; i < n; ++i) x[i] = (x[i] + P) * in % P;
}
}; // namespace
ll ta[N], tb[N];
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
    for (int i = 0; i <= n; ++i) ta[i] = a[i];
    for (int i = 0; i <= m; ++i) tb[i] = b[i];
    int len = 2;
    while (len < n + m + 2) len <<= 1;
    dft(len, ta), dft(len, tb);
    for (int i = 0; i < len; ++i) ta[i] = ta[i] * tb[i] % P;
    idft(len, ta);
    for (int i = 0; i <= n + m; ++i) c[i] = ta[i];
}

评测结果