代码

// NTT mod 998244353, scalar Montgomery
typedef unsigned long long u64;
typedef unsigned u32;

const u32 MOD = 998244353u;
const u32 NINV = 998244351u;    // -MOD^-1 mod 2^32
const u32 ONE  = 301989884u;    // 2^32 mod MOD
const u32 R2   = 932051910u;    // 2^64 mod MOD
const int MAXL = 1 << 21;

static u32 A[MAXL];
static u32 B[MAXL];
static u32 roots[MAXL];
static u32 roots_inv[MAXL];

static inline u32 mont_mul(u32 x, u32 y) {
    u64 t = (u64)x * y;
    u32 m = (u32)t * NINV;
    u64 u = (t + (u64)m * MOD) >> 32;
    if (u >= MOD) u -= MOD;
    return (u32)u;
}
static inline u32 add_mod(u32 x, u32 y) { u32 s = x + y; return s >= MOD ? s - MOD : s; }
static inline u32 sub_mod(u32 x, u32 y) { return x >= y ? x - y : x + MOD - y; }
static inline u32 to_mont(u32 x) { return mont_mul(x, R2); }

static u32 mont_pow(u32 base, u64 e) {
    u32 r = ONE, b = base;
    while (e) { if (e & 1) r = mont_mul(r, b); b = mont_mul(b, b); e >>= 1; }
    return r;
}

static void ntt(u32 *x, int n, const u32 *rts) {
    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit;
        if (i < j) { u32 t = x[i]; x[i] = x[j]; x[j] = t; }
    }
    for (int len = 2; len <= n; len <<= 1) {
        int half = len >> 1;
        int step = n / len;
        for (int i = 0; i < n; i += len) {
            u32 *y = x + i;
            for (int j = 0; j < half; j++) {
                u32 u = y[j];
                u32 v = mont_mul(y[j + half], rts[j * step]);
                y[j] = add_mod(u, v);
                y[j + half] = sub_mod(u, v);
            }
        }
    }
}

void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
    int L = 1;
    while (L < n + m + 2) L <<= 1;
    for (int i = 0; i <= n; i++) A[i] = to_mont(a[i]);
    for (int i = 0; i <= m; i++) B[i] = to_mont(b[i]);
    for (int i = n + 1; i < L; i++) A[i] = 0;
    for (int i = m + 1; i < L; i++) B[i] = 0;

    u32 w = mont_pow(to_mont(3), (MOD - 1) / L);
    u32 wi = mont_pow(w, MOD - 2);
    u32 cur = ONE;
    for (int i = 0; i < L; i++) { roots[i] = cur; cur = mont_mul(cur, w); }
    cur = ONE;
    for (int i = 0; i < L; i++) { roots_inv[i] = cur; cur = mont_mul(cur, wi); }

    ntt(A, L, roots);
    ntt(B, L, roots);
    for (int i = 0; i < L; i++) A[i] = mont_mul(A[i], B[i]);
    ntt(A, L, roots_inv);

    u32 linv_mont = mont_pow(to_mont((u32)(L % MOD)), MOD - 2);  // L^{-1}·R
    u32 linv_std = mont_mul(linv_mont, 1);                        // L^{-1} standard
    for (int i = 0; i <= n + m; i++) c[i] = mont_mul(A[i], linv_std);
}

评测结果