代码

// 1002i: polynomial multiply.  AVX2 NTT mod 998244353 (DIF forward / DIT inverse,
// Shoup twiddles, batched 8-point tail kernel; inverse scaled inside its last stage).
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <immintrin.h>
#include <sys/auxv.h>
typedef unsigned long long U64;
struct DI { U64 abi; const char *in; U64 insz; char *out; U64 outlim; U64 outsz; char *err;
            U64 errlim; U64 errsz; const char *IB; U64 IBlim; char *OB; U64 OBlim; U64 tscfreq; }
            __attribute__((packed));

// Fast AVX2 NTT mod 998244353.
//   fwd(a,n): natural order -> bit-reversed order  (DIF, vectorized over the twiddle index)
//   inv(a,n): bit-reversed  -> natural order       (DIT, same but conjugated)
// The last three DIF stages (len = 8,4,2) are fused into one branch-free 8-point
// kernel per contiguous group of 8, which removes the scalar tail that otherwise
// dominates small sizes.  Pointwise multiplication of two spectra is valid because
// both are in the same bit-reversed order and the inverse consumes that order.

typedef unsigned u32;
typedef unsigned long long u64;
static const u32 MOD = 998244353;

#define AVX2FN __attribute__((target("avx2")))

static inline u32 fpow(u32 a, u32 e) {
    u64 r = 1, b = a;
    while (e) { if (e & 1) r = r * b % MOD; b = b * b % MOD; e >>= 1; }
    return (u32)r;
}
AVX2FN static inline __m256i mulhi8(__m256i a, __m256i b) {
    __m256i lo = _mm256_mul_epu32(a, b);
    __m256i a1 = _mm256_srli_epi64(a, 32), b1 = _mm256_srli_epi64(b, 32);
    __m256i hi = _mm256_mul_epu32(a1, b1);
    return _mm256_blend_epi32(_mm256_srli_epi64(lo, 32), hi, 0xAA);
}
AVX2FN static inline __m256i shoup8(__m256i a, __m256i w, __m256i ws) {
    __m256i q = mulhi8(a, ws);
    __m256i r = _mm256_sub_epi32(_mm256_mullo_epi32(a, w), _mm256_mullo_epi32(q, _mm256_set1_epi32((int)MOD)));
    __m256i ge = _mm256_cmpgt_epi32(r, _mm256_set1_epi32((int)MOD - 1));
    return _mm256_sub_epi32(r, _mm256_and_si256(ge, _mm256_set1_epi32((int)MOD)));
}
AVX2FN static inline __m256i addmod8(__m256i a, __m256i b) {
    __m256i s = _mm256_add_epi32(a, b);
    return _mm256_sub_epi32(s, _mm256_and_si256(_mm256_cmpgt_epi32(s, _mm256_set1_epi32((int)MOD - 1)), _mm256_set1_epi32((int)MOD)));
}
AVX2FN static inline __m256i submod8(__m256i a, __m256i b) {
    __m256i lt = _mm256_cmpgt_epi32(b, a);
    return _mm256_add_epi32(_mm256_sub_epi32(a, b), _mm256_and_si256(lt, _mm256_set1_epi32((int)MOD)));
}
static inline u32 addm(u32 a, u32 b) { u32 s = a + b; return s >= MOD ? s - MOD : s; }
static inline u32 subm(u32 a, u32 b) { return a >= b ? a - b : a + MOD - b; }
static inline u32 mulm(u32 a, u32 b) { return (u32)((u64)a * b % MOD); }

struct Tables2 {
    int n = 0;
    static u32 fw[1 << 19], fws[1 << 19], iw[1 << 19], iws[1 << 19];   // static: judge links -static
    int foff[40] = {0}, ioff[40] = {0};
    // build the twiddle sequence for stage len: values w_len^j, j = 0..half-1
    static void build(u32 *dst, int half, u32 wlen) {
        u32 wls = (u32)(((__uint128_t)wlen << 32) / MOD);        // Shoup constant for wlen
        u64 cur = 1;
        for (int j = 0; j < half; j++) {
            dst[j] = (u32)cur;
            u32 q = (u32)(((u64)cur * wls) >> 32);
            u32 r = (u32)((u64)cur * wlen - (u64)q * MOD);
            cur = r >= MOD ? r - MOD : r;
        }
    }
    // Build the forward and inverse twiddle chains together.  Each chain is a serial
    // modular-multiply dependency chain; interleaving two of them lets the out-of-order
    // core overlap the latencies, which is where most of the 16M-tick build cost goes.
    // Fused build: both twiddle chains run in lockstep (independent serial
    // modular-multiply chains -> ILP) and each entry's Shoup constant
    // ws = floor(w * 2^32 / MOD) is produced in the same iteration, so the extra
    // O(N) pass over both tables (and its memory traffic) disappears entirely.
    static void build2(u32 *d0, u32 *s0out, u32 *d1, u32 *s1out, int half, u32 w0, u32 w1) {
        u32 s0 = (u32)(((__uint128_t)w0 << 32) / MOD);
        u32 s1 = (u32)(((__uint128_t)w1 << 32) / MOD);
        const double inv = 4294967296.0 / (double)MOD;
        u64 c0 = 1, c1 = 1;
        for (int j = 0; j < half; j++) {
            {   // ws for the forward entry
                u32 x = (u32)c0;
                u32 q = (u32)((double)x * inv);
                u64 lhs = (u64)q * MOD, rhs = (u64)x << 32;
                if (lhs > rhs) q--; else if (lhs + MOD <= rhs) q++;
                s0out[j] = q;
            }
            d0[j] = (u32)c0;
            {   // ws for the inverse entry (independent -> overlaps the chain above)
                u32 x = (u32)c1;
                u32 q = (u32)((double)x * inv);
                u64 lhs = (u64)q * MOD, rhs = (u64)x << 32;
                if (lhs > rhs) q--; else if (lhs + MOD <= rhs) q++;
                s1out[j] = q;
            }
            d1[j] = (u32)c1;
            u32 qa = (u32)((c0 * s0) >> 32);
            u32 ra = (u32)(c0 * w0 - (u64)qa * MOD);
            c0 = ra >= MOD ? ra - MOD : ra;
            u32 qb = (u32)((c1 * s1) >> 32);
            u32 rb = (u32)(c1 * w1 - (u64)qb * MOD);
            c1 = rb >= MOD ? rb - MOD : rb;
        }
    }
    // ws = floor(w * 2^32 / MOD) using a double reciprocal plus exact correction
    static void buildws(u32 *dst, const u32 *w, int half) {
        const double inv = 4294967296.0 / (double)MOD;
        for (int j = 0; j < half; j++) {
            u32 x = w[j];
            u32 q = (u32)((double)x * inv);
            u64 lhs = (u64)q * MOD, rhs = (u64)x << 32;
            if (lhs > rhs) q--;
            else if (lhs + MOD <= rhs) q++;
            dst[j] = q;
        }
    }
    void init(int n_) {
        n = n_;
        if (n > (1 << 19)) n = 1 << 19;
        int idx = 0, idx2 = 0;
        for (int len = n; len >= 16; len >>= 1) {
            int b = __builtin_ctz(len);
            foff[b] = idx; ioff[b] = idx2;
            int half = len >> 1;
            build2(fw + idx, fws + idx, iw + idx2, iws + idx2, half,
                   fpow(3, (MOD - 1) / len),
                   fpow(3, MOD - 1 - (MOD - 1) / len));
            idx += half; idx2 += half;
        }
    }
};
u32 Tables2::fw[1 << 19], Tables2::fws[1 << 19], Tables2::iw[1 << 19], Tables2::iws[1 << 19];


AVX2FN static inline void tr8i(__m256i *r) {
    __m256i t0=_mm256_unpacklo_epi32(r[0],r[1]), t1=_mm256_unpackhi_epi32(r[0],r[1]);
    __m256i t2=_mm256_unpacklo_epi32(r[2],r[3]), t3=_mm256_unpackhi_epi32(r[2],r[3]);
    __m256i t4=_mm256_unpacklo_epi32(r[4],r[5]), t5=_mm256_unpackhi_epi32(r[4],r[5]);
    __m256i t6=_mm256_unpacklo_epi32(r[6],r[7]), t7=_mm256_unpackhi_epi32(r[6],r[7]);
    __m256i u0=_mm256_unpacklo_epi64(t0,t2), u1=_mm256_unpackhi_epi64(t0,t2);
    __m256i u2=_mm256_unpacklo_epi64(t1,t3), u3=_mm256_unpackhi_epi64(t1,t3);
    __m256i u4=_mm256_unpacklo_epi64(t4,t6), u5=_mm256_unpackhi_epi64(t4,t6);
    __m256i u6=_mm256_unpacklo_epi64(t5,t7), u7=_mm256_unpackhi_epi64(t5,t7);
    r[0]=_mm256_permute2x128_si256(u0,u4,0x20); r[1]=_mm256_permute2x128_si256(u1,u5,0x20);
    r[2]=_mm256_permute2x128_si256(u2,u6,0x20); r[3]=_mm256_permute2x128_si256(u3,u7,0x20);
    r[4]=_mm256_permute2x128_si256(u0,u4,0x31); r[5]=_mm256_permute2x128_si256(u1,u5,0x31);
    r[6]=_mm256_permute2x128_si256(u2,u6,0x31); r[7]=_mm256_permute2x128_si256(u3,u7,0x31);
}
static inline void k8fwd_scalar(u32 *x, const u32 *w8, const u32 *w4) {
    for (int j = 0; j < 4; j++) { u32 u = x[j], v = x[j + 4]; x[j] = addm(u, v); x[j + 4] = mulm(subm(u, v), w8[j]); }
    for (int b = 0; b < 8; b += 4) for (int j = 0; j < 2; j++) { u32 u = x[b + j], v = x[b + j + 2]; x[b + j] = addm(u, v); x[b + j + 2] = mulm(subm(u, v), w4[j]); }
    for (int b = 0; b < 8; b += 2) { u32 u = x[b], v = x[b + 1]; x[b] = addm(u, v); x[b + 1] = subm(u, v); }
}
static inline void k8inv_scalar(u32 *x, const u32 *w8i, const u32 *w4i) {
    for (int b = 0; b < 8; b += 2) { u32 u = x[b], v = x[b + 1]; x[b] = addm(u, v); x[b + 1] = subm(u, v); }
    for (int b = 0; b < 8; b += 4) for (int j = 0; j < 2; j++) { u32 u = x[b + j], v = mulm(x[b + j + 2], w4i[j]); x[b + j] = addm(u, v); x[b + j + 2] = subm(u, v); }
    for (int j = 0; j < 4; j++) { u32 u = x[j], v = mulm(x[j + 4], w8i[j]); x[j] = addm(u, v); x[j + 4] = subm(u, v); }
}
// generic scalar DIF (natural -> bit-reversed) for tiny n
static void small_fwd(u32 *a, int n) {
    for (int len = n; len >= 2; len >>= 1) {
        int half = len >> 1;
        u32 lw = fpow(3, (MOD - 1) / len);
        for (int i = 0; i < n; i += len) {
            u64 cur = 1;
            for (int j = 0; j < half; j++) {
                u32 u = a[i + j], v = a[i + j + half];
                a[i + j] = addm(u, v); a[i + j + half] = mulm(subm(u, v), (u32)cur);
                cur = cur * lw % MOD;
            }
        }
    }
}
static void small_inv(u32 *a, int n) {
    for (int len = 2; len <= n; len <<= 1) {
        int half = len >> 1;
        u32 lw = fpow(3, MOD - 1 - (MOD - 1) / len);
        for (int i = 0; i < n; i += len) {
            u64 cur = 1;
            for (int j = 0; j < half; j++) {
                u32 u = a[i + j], x = a[i + j + half];
                u32 v = mulm(x, (u32)cur);
                a[i + j] = addm(u, v); a[i + j + half] = subm(u, v);
                cur = cur * lw % MOD;
            }
        }
    }
    u32 ninv = fpow((u32)n, MOD - 2);
    for (int i = 0; i < n; i++) a[i] = mulm(a[i], ninv);
}

AVX2FN static void tail8fwd(u32 *a, int n) {
    u32 g8 = fpow(3, (MOD - 1) / 8), g4 = fpow(3, (MOD - 1) / 4);
    __m256i W8[4], S8[4], W4[2], S4[2];
    for (int j = 0; j < 4; j++) { u32 w = fpow(g8, (u32)j); W8[j] = _mm256_set1_epi32((int)w);
        S8[j] = _mm256_set1_epi32((int)(u32)(((__uint128_t)w << 32) / MOD)); }
    for (int j = 0; j < 2; j++) { u32 w = fpow(g4, (u32)j); W4[j] = _mm256_set1_epi32((int)w);
        S4[j] = _mm256_set1_epi32((int)(u32)(((__uint128_t)w << 32) / MOD)); }
    u32 w8c[4], w4c[2];
    for (int j = 0; j < 4; j++) w8c[j] = fpow(g8, (u32)j);
    for (int j = 0; j < 2; j++) w4c[j] = fpow(g4, (u32)j);
    for (int i = 0; i + 64 <= n; i += 64) {
        __m256i v[8];
        for (int m = 0; m < 8; m++) v[m] = _mm256_loadu_si256((const __m256i*)(const void*)(a + i + 8 * m));
        tr8i(v);
        for (int j = 0; j < 4; j++) { __m256i u = v[j], x = v[j + 4]; v[j] = addmod8(u, x); v[j + 4] = shoup8(submod8(u, x), W8[j], S8[j]); }
        for (int b = 0; b < 8; b += 4) for (int j = 0; j < 2; j++) { __m256i u = v[b + j], x = v[b + j + 2]; v[b + j] = addmod8(u, x); v[b + j + 2] = shoup8(submod8(u, x), W4[j], S4[j]); }
        for (int b = 0; b < 8; b += 2) { __m256i u = v[b], x = v[b + 1]; v[b] = addmod8(u, x); v[b + 1] = submod8(u, x); }
        tr8i(v);
        for (int m = 0; m < 8; m++) _mm256_storeu_si256((__m256i*)(void*)(a + i + 8 * m), v[m]);
    }
    for (int i = (n / 64) * 64; i + 8 <= n; i += 8) k8fwd_scalar(a + i, w8c, w4c);
}
AVX2FN static void tail8inv(u32 *a, int n) {
    u32 g8 = fpow(3, MOD - 1 - (MOD - 1) / 8), g4 = fpow(3, MOD - 1 - (MOD - 1) / 4);
    __m256i W8[4], S8[4], W4[2], S4[2];
    for (int j = 0; j < 4; j++) { u32 w = fpow(g8, (u32)j); W8[j] = _mm256_set1_epi32((int)w);
        S8[j] = _mm256_set1_epi32((int)(u32)(((__uint128_t)w << 32) / MOD)); }
    for (int j = 0; j < 2; j++) { u32 w = fpow(g4, (u32)j); W4[j] = _mm256_set1_epi32((int)w);
        S4[j] = _mm256_set1_epi32((int)(u32)(((__uint128_t)w << 32) / MOD)); }
    for (int i = 0; i + 64 <= n; i += 64) {
        __m256i v[8];
        for (int m = 0; m < 8; m++) v[m] = _mm256_loadu_si256((const __m256i*)(const void*)(a + i + 8 * m));
        tr8i(v);
        for (int b = 0; b < 8; b += 2) { __m256i u = v[b], x = v[b + 1]; v[b] = addmod8(u, x); v[b + 1] = submod8(u, x); }
        for (int b = 0; b < 8; b += 4) for (int j = 0; j < 2; j++) { __m256i u = v[b + j], x = shoup8(v[b + j + 2], W4[j], S4[j]); v[b + j] = addmod8(u, x); v[b + j + 2] = submod8(u, x); }
        for (int j = 0; j < 4; j++) { __m256i u = v[j], x = shoup8(v[j + 4], W8[j], S8[j]); v[j] = addmod8(u, x); v[j + 4] = submod8(u, x); }
        tr8i(v);
        for (int m = 0; m < 8; m++) _mm256_storeu_si256((__m256i*)(void*)(a + i + 8 * m), v[m]);
    }
    u32 w8i[4], w4i[2];
    for (int j = 0; j < 4; j++) w8i[j] = fpow(g8, (u32)j);
    for (int j = 0; j < 2; j++) w4i[j] = fpow(g4, (u32)j);
    for (int i = (n / 64) * 64; i + 8 <= n; i += 8) k8inv_scalar(a + i, w8i, w4i);
}

AVX2FN static void ntt2_fwd(u32 *a, int n, const Tables2 &T) {
    if (n < 16) { small_fwd(a, n); return; }
    for (int len = n; len >= 16; len >>= 1) {
        int half = len >> 1;
        const u32 *w = T.fw + T.foff[__builtin_ctz(len)];
        const u32 *ws = T.fws + T.foff[__builtin_ctz(len)];
        for (int i = 0; i < n; i += len) {
            u32 *lo = a + i, *hi = lo + half;
            int j = 0;
            for (; j + 8 <= half; j += 8) {
                __m256i u = _mm256_loadu_si256((const __m256i*)(const void*)(lo + j));
                __m256i v = _mm256_loadu_si256((const __m256i*)(const void*)(hi + j));
                __m256i W = _mm256_loadu_si256((const __m256i*)(const void*)(w + j));
                _mm256_storeu_si256((__m256i*)(void*)(lo + j), addmod8(u, v));
                _mm256_storeu_si256((__m256i*)(void*)(hi + j), shoup8(submod8(u, v), W, _mm256_loadu_si256((const __m256i*)(const void*)(ws + j))));
            }
            for (; j < half; j++) {
                u32 u = lo[j], v = hi[j];
                u32 wsj = (u32)(((__uint128_t)w[j] << 32) / MOD);
                u32 q = (u32)(((u64)subm(u, v) * wsj) >> 32);
                u32 r = (u32)((u64)subm(u, v) * w[j] - (u64)q * MOD);
                if (r >= MOD) r -= MOD;
                lo[j] = addm(u, v); hi[j] = r;
            }
        }
    }
    tail8fwd(a, n);
}
AVX2FN static void ntt2_inv(u32 *a, int n, const Tables2 &T) {
    if (n < 16) { small_inv(a, n); return; }
    u32 ninv = fpow((u32)n, MOD - 2);
    u32 ninvs = (u32)(((__uint128_t)ninv << 32) / MOD);
    __m256i NINV = _mm256_set1_epi32((int)ninv), NINVS = _mm256_set1_epi32((int)ninvs);
    tail8inv(a, n);
    for (int len = 16; len <= n; len <<= 1) {
        int half = len >> 1;
        const u32 *w = T.iw + T.ioff[__builtin_ctz(len)];
        const u32 *ws = T.iws + T.ioff[__builtin_ctz(len)];
        for (int i = 0; i < n; i += len) {
            u32 *lo = a + i, *hi = lo + half;
            int j = 0;
            for (; j + 8 <= half; j += 8) {
                __m256i u = _mm256_loadu_si256((const __m256i*)(const void*)(lo + j));
                __m256i x = _mm256_loadu_si256((const __m256i*)(const void*)(hi + j));
                __m256i W = _mm256_loadu_si256((const __m256i*)(const void*)(w + j));
                if (len == n) { u = shoup8(u, NINV, NINVS); x = shoup8(x, NINV, NINVS); }
                __m256i v = shoup8(x, W, _mm256_loadu_si256((const __m256i*)(const void*)(ws + j)));
                _mm256_storeu_si256((__m256i*)(void*)(lo + j), addmod8(u, v));
                _mm256_storeu_si256((__m256i*)(void*)(hi + j), submod8(u, v));
            }
            for (; j < half; j++) {
                u32 u = lo[j], x = hi[j];
                if (len == n) { u = mulm(u, ninv); x = mulm(x, ninv); }
                u32 wsj = (u32)(((__uint128_t)w[j] << 32) / MOD);
                u32 q = (u32)(((u64)x * wsj) >> 32);
                u32 v = (u32)((u64)x * w[j] - (u64)q * MOD);
                if (v >= MOD) v -= MOD;
                lo[j] = addm(u, v); hi[j] = subm(u, v);
            }
        }
    }
}


static u32 A[1 << 18], B[1 << 18];
static inline int rdInt(char *&p) {
    while (*p < '0' || *p > '9') p++;
    // SWAR fast path: eight ASCII digits at once
    unsigned long long x;
    memcpy(&x, p, 8);
    unsigned long long d = x - 0x3030303030303030ULL;
    if (!((d | (d + 0x7676767676767676ULL)) & 0x8080808080808080ULL)) {
        unsigned long long t = (d * 10 + (d >> 8)) & 0x00FF00FF00FF00FFULL;
        t = (t * 100 + (t >> 16)) & 0x0000FFFF0000FFFFULL;
        t = (t * 10000 + (t >> 32)) & 0x00000000FFFFFFFFULL;
        p += 8;
        int v = (int)(unsigned)t;
        while (*p >= '0' && *p <= '9') v = v * 10 + (*p++ - '0');
        return v;
    }
    int v = 0;
    while (*p >= '0' && *p <= '9') v = v * 10 + (*p++ - '0');
    return v;
}

static char outbuf[1 << 22];
static char *op = outbuf;
static const char DIG[201] =
    "00010203040506070809" "10111213141516171819" "20212223242526272829"
    "30313233343536373839" "40414243444546474849" "50515253545556575859"
    "60616263646566676869" "70717273747576777879" "80818283848586878889"
    "90919293949596979899";

static inline void writeInt(u32 v) {
    char t[12]; int n = 0;
    while (v >= 100) { u32 r = v % 100; v /= 100; t[n++] = DIG[2 * r + 1]; t[n++] = DIG[2 * r]; }
    if (v >= 10) { t[n++] = DIG[2 * v + 1]; t[n++] = DIG[2 * v]; }
    else t[n++] = char('0' + v);
    while (n) *op++ = t[--n];
}

int main() {
    static char inbuf[1 << 23];
    int total = (int)fread(inbuf, 1, sizeof(inbuf) - 32, stdin);
    for (int z = 0; z < 32; z++) inbuf[total + z] = 0;   /* SWAR reads up to 8 bytes ahead */
    char *p = inbuf;
    int n = rdInt(p), m = rdInt(p);
    int N = 1; while (N < n + m + 1) N <<= 1;
    for (int i = 0; i <= n; i++) A[i] = (u32)rdInt(p);
    for (int i = 0; i <= m; i++) B[i] = (u32)rdInt(p);

    Tables2 T; T.init(N);
    ntt2_fwd(A, N, T);
    ntt2_fwd(B, N, T);
    for (int i = 0; i < N; i++) A[i] = (u32)((u64)A[i] * B[i] % MOD);
    ntt2_inv(A, N, T);
    DI *di = (DI *)getauxval(0x6b637564);
    char *obase = (di && di->out) ? di->out : outbuf;
    op = obase;
    for (int i = 0; i <= n + m; i++) {
        if (i) *op++ = ' ';
        writeInt(A[i]);
    }
    *op++ = '\n';
    if (di && di->out) di->outsz = (U64)(op - obase);
    else fwrite(outbuf, 1, op - outbuf, stdout);
    return 0;
}

评测结果