代码

#pragma GCC optimize("O3,unroll-loops,omit-frame-pointer")
#pragma GCC target("arch=skylake")

#ifndef DUCK_FASTIO_H
#define DUCK_FASTIO_H

typedef unsigned long duck_u64;
typedef long duck_i64;

typedef struct {
    duck_u64 abi_version;
    const char *stdin_ptr;
    duck_u64 stdin_size;
    char *stdout_ptr;
    duck_u64 stdout_limit;
    duck_u64 stdout_size;
    char *stderr_ptr;
    duck_u64 stderr_limit;
    duck_u64 stderr_size;
    const char *ib_ptr;
    duck_u64 ib_limit;
    char *ob_ptr;
    duck_u64 ob_limit;
    duck_u64 tsc_frequency;
} __attribute__((packed)) DuckInfo;

static __attribute__((always_inline)) inline DuckInfo *duck_info(long argc, char **argv) {
    char **p = argv + argc + 1;
    while (*p) ++p;
    duck_u64 *aux = (duck_u64 *)(p + 1);
    while (aux[0]) {
        if (aux[0] == 0x6b637564UL) return (DuckInfo *)aux[1];
        aux += 2;
    }
    return (DuckInfo *)0;
}

static __attribute__((always_inline)) inline duck_u64 duck_read_u64(const char **cursor) {
    const char *p = *cursor;
    while ((unsigned char)(*p - '0') > 9) ++p;
    duck_u64 value = 0;
    do {
        value = value * 10 + (unsigned char)(*p - '0');
        ++p;
    } while ((unsigned char)(*p - '0') <= 9);
    *cursor = p;
    return value;
}

static __attribute__((always_inline)) inline duck_i64 duck_read_i64(const char **cursor) {
    const char *p = *cursor;
    while (*p != '-' && (unsigned char)(*p - '0') > 9) ++p;
    int negative = *p == '-';
    p += negative;
    duck_u64 value = 0;
    do {
        value = value * 10 + (unsigned char)(*p - '0');
        ++p;
    } while ((unsigned char)(*p - '0') <= 9);
    *cursor = p;
    return negative ? -(duck_i64)value : (duck_i64)value;
}

static __attribute__((always_inline)) inline char *duck_write_u64(char *out, duck_u64 value) {
    char tmp[24];
    unsigned n = 0;
    do {
        tmp[n++] = (char)('0' + value % 10);
        value /= 10;
    } while (value);
    do *out++ = tmp[--n]; while (n);
    return out;
}

static __attribute__((always_inline)) inline char *duck_write_i64(char *out, duck_i64 value) {
    if (value < 0) {
        *out++ = '-';
        return duck_write_u64(out, (duck_u64)(-value));
    }
    return duck_write_u64(out, (duck_u64)value);
}

static __attribute__((always_inline, noreturn)) inline void duck_exit(void) {
    __asm__ volatile("mov $60,%%eax;xor %%edi,%%edi;syscall" ::: "rax", "rdi", "rcx", "r11", "memory");
    __builtin_unreachable();
}

#endif

#include <immintrin.h>

typedef unsigned int u32;
typedef unsigned long u64;

enum { MOD = 998244353u, MAXL = 1 << 21 };
static u32 fa[MAXL], fb[MAXL];
static u32 twiddle[MAXL / 2];
static u32 twiddle_shoup[MAXL / 2];

static __attribute__((always_inline)) inline u32 mont_mul(u32 a, u32 b) {
    u64 t = (u64)a * b;
    u32 m = (u32)t * 0x3b7fffffu;
    return (u32)((t + (u64)m * MOD) >> 32);
}

static u32 mod_pow(u32 a, u32 e) {
    u64 r = 1;
    while (e) {
        if (e & 1) r = r * a % MOD;
        a = (u32)((u64)a * a % MOD);
        e >>= 1;
    }
    return (u32)r;
}

__attribute__((target("avx2"), always_inline))
static inline __m256i mont_mul8(__m256i a, __m256i b) {
    const __m256i qi = _mm256_set1_epi32((int)0x3b7fffffu);
    const __m256i qm = _mm256_set1_epi32((int)MOD);
    __m256i te = _mm256_mul_epu32(a, b);
    __m256i to = _mm256_mul_epu32(_mm256_srli_epi64(a, 32),
                                  _mm256_srli_epi64(b, 32));
    __m256i me = _mm256_mul_epu32(te, qi);
    __m256i mo = _mm256_mul_epu32(to, qi);
    __m256i ue = _mm256_srli_epi64(
        _mm256_add_epi64(te, _mm256_mul_epu32(me, qm)), 32);
    __m256i uo = _mm256_srli_epi64(
        _mm256_add_epi64(to, _mm256_mul_epu32(mo, qm)), 32);
    return _mm256_blend_epi32(ue, _mm256_slli_si256(uo, 4), 0xaa);
}

__attribute__((target("avx2"), always_inline))
static inline __m256i mulhi_u32x8(__m256i a, __m256i b) {
    __m256i pe = _mm256_srli_epi64(_mm256_mul_epu32(a, b), 32);
    __m256i po = _mm256_srli_epi64(
        _mm256_mul_epu32(_mm256_srli_epi64(a, 32),
                         _mm256_srli_epi64(b, 32)), 32);
    return _mm256_blend_epi32(pe, _mm256_slli_si256(po, 4), 0xaa);
}

__attribute__((target("avx2"), always_inline))
static inline __m256i shoup_mul8(__m256i a, __m256i w, __m256i ws) {
    __m256i qhat = mulhi_u32x8(a, ws);
    return _mm256_sub_epi32(_mm256_mullo_epi32(a, w),
                            _mm256_mullo_epi32(qhat,
                                              _mm256_set1_epi32((int)MOD)));
}

__attribute__((target("avx2"), always_inline))
static inline void forward_block(u32 *a, int p, int half, int j) {
    const __m256i m2 = _mm256_set1_epi32((int)(2u * MOD));
    __m256i x = _mm256_loadu_si256((const __m256i *)(a + p + j));
    __m256i y = _mm256_loadu_si256((const __m256i *)(a + p + j + half));
    __m256i sum = _mm256_add_epi32(x, y);
    sum = _mm256_min_epu32(sum, _mm256_sub_epi32(sum, m2));
    __m256i diff = _mm256_sub_epi32(_mm256_add_epi32(x, m2), y);
    diff = _mm256_min_epu32(diff, _mm256_sub_epi32(diff, m2));
    __m256i w = _mm256_loadu_si256((const __m256i *)(twiddle + j));
    _mm256_storeu_si256((__m256i *)(a + p + j), sum);
    _mm256_storeu_si256((__m256i *)(a + p + j + half),
                        mont_mul8(diff, w));
}

__attribute__((target("avx2"), always_inline))
static inline __m256i reduce_sum(__m256i x, __m256i y) {
    __m256i m2 = _mm256_set1_epi32((int)(2u * MOD));
    __m256i z = _mm256_add_epi32(x, y);
    return _mm256_min_epu32(z, _mm256_sub_epi32(z, m2));
}

__attribute__((target("avx2"), always_inline))
static inline __m256i reduce_diff(__m256i x, __m256i y) {
    __m256i m2 = _mm256_set1_epi32((int)(2u * MOD));
    __m256i z = _mm256_sub_epi32(_mm256_add_epi32(x, m2), y);
    return _mm256_min_epu32(z, _mm256_sub_epi32(z, m2));
}

/* Fuse two consecutive DIF layers.  Besides halving the array traffic, this
   lets both transforms share the same pair of packed twiddle streams. */
__attribute__((target("avx2"), always_inline))
static inline void forward_radix4_block(u32 *a, int p, int quarter, int j) {
    __m256i x0 = _mm256_loadu_si256((const __m256i *)(a + p + j));
    __m256i x1 = _mm256_loadu_si256((const __m256i *)(a + p + j + quarter));
    __m256i x2 = _mm256_loadu_si256((const __m256i *)(a + p + j + 2 * quarter));
    __m256i x3 = _mm256_loadu_si256((const __m256i *)(a + p + j + 3 * quarter));
    __m256i s0 = reduce_sum(x0, x2);
    __m256i s1 = reduce_sum(x1, x3);
    __m256i w0 = _mm256_loadu_si256((const __m256i *)(twiddle + j));
    __m256i w1 = _mm256_loadu_si256((const __m256i *)(twiddle + quarter + j));
    __m256i w2 = _mm256_loadu_si256((const __m256i *)(twiddle_shoup + j));
    __m256i d0 = mont_mul8(reduce_diff(x0, x2), w0);
    __m256i d1 = mont_mul8(reduce_diff(x1, x3), w1);
    _mm256_storeu_si256((__m256i *)(a + p + j), reduce_sum(s0, s1));
    _mm256_storeu_si256((__m256i *)(a + p + j + quarter),
                        mont_mul8(reduce_diff(s0, s1), w2));
    _mm256_storeu_si256((__m256i *)(a + p + j + 2 * quarter),
                        reduce_sum(d0, d1));
    _mm256_storeu_si256((__m256i *)(a + p + j + 3 * quarter),
                        mont_mul8(reduce_diff(d0, d1), w2));
}

__attribute__((target("avx2")))
static void forward_fused16(u32 *a, int n) {
    __m256i w1 = _mm256_loadu_si256((const __m256i *)twiddle);
    __m256i w2 = _mm256_setr_epi32(
        (int)twiddle[0], (int)twiddle[2], (int)twiddle[4], (int)twiddle[6],
        (int)twiddle[0], (int)twiddle[2], (int)twiddle[4], (int)twiddle[6]);
    for (int i = 0; i < n; i += 16) {
        __m256i x = _mm256_loadu_si256((const __m256i *)(a + i));
        __m256i y = _mm256_loadu_si256((const __m256i *)(a + i + 8));
        __m256i s = reduce_sum(x, y);
        __m256i d = mont_mul8(reduce_diff(x, y), w1);
        __m256i sl = _mm256_permute2x128_si256(s, s, 0x00);
        __m256i sh = _mm256_permute2x128_si256(s, s, 0x11);
        __m256i dl = _mm256_permute2x128_si256(d, d, 0x00);
        __m256i dh = _mm256_permute2x128_si256(d, d, 0x11);
        __m256i sums = _mm256_permute2x128_si256(
            reduce_sum(sl, sh), reduce_sum(dl, dh), 0x20);
        __m256i diffs = _mm256_permute2x128_si256(
            reduce_diff(sl, sh), reduce_diff(dl, dh), 0x20);
        diffs = mont_mul8(diffs, w2);
        _mm256_storeu_si256((__m256i *)(a + i),
                            _mm256_permute2x128_si256(sums, diffs, 0x20));
        _mm256_storeu_si256((__m256i *)(a + i + 8),
                            _mm256_permute2x128_si256(sums, diffs, 0x31));
    }
}

__attribute__((target("avx2")))
static void forward_fused4(u32 *a, int n) {
    const __m256i ix = _mm256_setr_epi32(0, 1, 0, 1, 4, 5, 4, 5);
    const __m256i iy = _mm256_setr_epi32(2, 3, 2, 3, 6, 7, 6, 7);
    const __m256i swap = _mm256_setr_epi32(1, 0, 3, 2, 5, 4, 7, 6);
    __m256i w = _mm256_setr_epi32(
        (int)twiddle[0], (int)twiddle[1], (int)twiddle[0], (int)twiddle[1],
        (int)twiddle[0], (int)twiddle[1], (int)twiddle[0], (int)twiddle[1]);
    for (int i = 0; i < n; i += 8) {
        __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
        __m256i x = _mm256_permutevar8x32_epi32(v, ix);
        __m256i y = _mm256_permutevar8x32_epi32(v, iy);
        __m256i q = _mm256_blend_epi32(
            reduce_sum(x, y), mont_mul8(reduce_diff(x, y), w), 0xcc);
        __m256i z = _mm256_permutevar8x32_epi32(q, swap);
        _mm256_storeu_si256((__m256i *)(a + i),
            _mm256_blend_epi32(reduce_sum(q, z), reduce_diff(z, q), 0xaa));
    }
}

__attribute__((target("avx2")))
static void forward_small(u32 *a, int n, int len) {
    if (len == 2) {
        const __m256i swap_index = _mm256_setr_epi32(1, 0, 3, 2, 5, 4, 7, 6);
        for (int i = 0; i < n; i += 8) {
            __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
            __m256i z = _mm256_permutevar8x32_epi32(v, swap_index);
            __m256i s = reduce_sum(v, z);
            __m256i d = reduce_diff(z, v);
            _mm256_storeu_si256((__m256i *)(a + i),
                                _mm256_blend_epi32(s, d, 0xaa));
        }
    } else if (len == 4) {
        __m256i w = _mm256_setr_epi32(
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[0], (int)twiddle[1]);
        for (int i = 0; i < n; i += 8) {
            __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
            __m256i x = _mm256_permute4x64_epi64(v, 0xa0);
            __m256i y = _mm256_permute4x64_epi64(v, 0xf5);
            __m256i s = reduce_sum(x, y);
            __m256i d = mont_mul8(reduce_diff(x, y), w);
            _mm256_storeu_si256((__m256i *)(a + i),
                                _mm256_blend_epi32(s, d, 0xcc));
        }
    } else {
        __m256i w = _mm256_setr_epi32(
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[2], (int)twiddle[3],
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[2], (int)twiddle[3]);
        for (int i = 0; i < n; i += 8) {
            __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
            __m256i x = _mm256_permute2x128_si256(v, v, 0x00);
            __m256i y = _mm256_permute2x128_si256(v, v, 0x11);
            __m256i s = reduce_sum(x, y);
            __m256i d = mont_mul8(reduce_diff(x, y), w);
            _mm256_storeu_si256((__m256i *)(a + i),
                _mm256_permute2x128_si256(s, d, 0x20));
        }
    }
}

__attribute__((target("avx2")))
static void forward_both(u32 *a, u32 *b, int n) {
    const u32 mod2 = 2u * MOD;
    int len = n;
    for (; len >= 32; len >>= 2) {
        int half = len >> 1, quarter = len >> 2;
        u32 step = mont_mul(mod_pow(3, (MOD - 1u) / (u32)len),
                            932051910u);
        twiddle[0] = 301989884u;
        for (int j = 1; j < half; ++j)
            twiddle[j] = mont_mul(twiddle[j - 1], step);
        for (int j = 0; j < quarter; ++j)
            twiddle_shoup[j] = twiddle[j << 1];
        for (int i = 0; i < n; i += len) {
            for (int j = 0; j < quarter; j += 8) {
                forward_radix4_block(a, i, quarter, j);
                forward_radix4_block(b, i, quarter, j);
            }
        }
    }
    if (len == 16) {
        u32 step = mont_mul(mod_pow(3, (MOD - 1u) / 16u), 932051910u);
        twiddle[0] = 301989884u;
        for (int j = 1; j < 8; ++j) twiddle[j] = mont_mul(twiddle[j - 1], step);
        forward_fused16(a, n);
        forward_fused16(b, n);
        len = 4;
    }
    if (len == 4) {
        u32 step = mont_mul(mod_pow(3, (MOD - 1u) / 4u), 932051910u);
        twiddle[0] = 301989884u;
        twiddle[1] = mont_mul(twiddle[0], step);
        forward_fused4(a, n);
        forward_fused4(b, n);
        len = 1;
    }
    for (; len >= 2; len >>= 1) {
        int half = len >> 1;
        u32 step = mont_mul(mod_pow(3, (MOD - 1u) / (u32)len),
                            932051910u);
        twiddle[0] = 301989884u;
        for (int j = 1; j < half; ++j)
            twiddle[j] = mont_mul(twiddle[j - 1], step);
        if (len <= 8) {
            forward_small(a, n, len);
            forward_small(b, n, len);
            continue;
        }
        for (int i = 0; i < n; i += len) {
            int j = 0;
            for (; j + 8 <= half; j += 8) {
                forward_block(a, i, half, j);
                forward_block(b, i, half, j);
            }
            for (; j < half; ++j) {
                u32 w = twiddle[j];
                u32 ax = a[i + j], ay = a[i + j + half];
                u32 as = ax + ay;
                if (as >= mod2) as -= mod2;
                u32 ad = ax + mod2 - ay;
                if (ad >= mod2) ad -= mod2;
                a[i + j] = as;
                a[i + j + half] = mont_mul(ad, w);
                u32 bx = b[i + j], by = b[i + j + half];
                u32 bs = bx + by;
                if (bs >= mod2) bs -= mod2;
                u32 bd = bx + mod2 - by;
                if (bd >= mod2) bd -= mod2;
                b[i + j] = bs;
                b[i + j + half] = mont_mul(bd, w);
            }
        }
    }
}

__attribute__((target("avx2"), always_inline))
static inline void inverse_radix4_block(u32 *a, int p, int quarter, int j) {
    __m256i x0 = _mm256_loadu_si256((const __m256i *)(a + p + j));
    __m256i x1 = _mm256_loadu_si256((const __m256i *)(a + p + j + quarter));
    __m256i x2 = _mm256_loadu_si256((const __m256i *)(a + p + j + 2 * quarter));
    __m256i x3 = _mm256_loadu_si256((const __m256i *)(a + p + j + 3 * quarter));
    __m256i w2 = _mm256_loadu_si256((const __m256i *)(twiddle_shoup + j));
    x1 = mont_mul8(x1, w2);
    x3 = mont_mul8(x3, w2);
    __m256i s0 = reduce_sum(x0, x1);
    __m256i d0 = reduce_diff(x0, x1);
    __m256i s1 = reduce_sum(x2, x3);
    __m256i d1 = reduce_diff(x2, x3);
    __m256i w0 = _mm256_loadu_si256((const __m256i *)(twiddle + j));
    __m256i w1 = _mm256_loadu_si256((const __m256i *)(twiddle + quarter + j));
    s1 = mont_mul8(s1, w0);
    d1 = mont_mul8(d1, w1);
    _mm256_storeu_si256((__m256i *)(a + p + j), reduce_sum(s0, s1));
    _mm256_storeu_si256((__m256i *)(a + p + j + quarter), reduce_sum(d0, d1));
    _mm256_storeu_si256((__m256i *)(a + p + j + 2 * quarter), reduce_diff(s0, s1));
    _mm256_storeu_si256((__m256i *)(a + p + j + 3 * quarter), reduce_diff(d0, d1));
}

__attribute__((target("avx2")))
static void inverse_fused4(u32 *a, int n) {
    const __m256i swap = _mm256_setr_epi32(1, 0, 3, 2, 5, 4, 7, 6);
    const __m256i ix = _mm256_setr_epi32(0, 1, 0, 1, 4, 5, 4, 5);
    const __m256i iy = _mm256_setr_epi32(2, 3, 2, 3, 6, 7, 6, 7);
    __m256i w = _mm256_setr_epi32(
        (int)twiddle[0], (int)twiddle[1], (int)twiddle[0], (int)twiddle[1],
        (int)twiddle[0], (int)twiddle[1], (int)twiddle[0], (int)twiddle[1]);
    for (int i = 0; i < n; i += 8) {
        __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
        __m256i z = _mm256_permutevar8x32_epi32(v, swap);
        __m256i q = _mm256_blend_epi32(
            reduce_sum(v, z), reduce_diff(z, v), 0xaa);
        __m256i x = _mm256_permutevar8x32_epi32(q, ix);
        __m256i y = mont_mul8(_mm256_permutevar8x32_epi32(q, iy), w);
        _mm256_storeu_si256((__m256i *)(a + i),
            _mm256_blend_epi32(reduce_sum(x, y), reduce_diff(x, y), 0xcc));
    }
}

__attribute__((target("avx2")))
static void inverse_fused16(u32 *a, int n) {
    __m256i w1 = _mm256_loadu_si256((const __m256i *)twiddle);
    __m256i w2 = _mm256_setr_epi32(
        (int)twiddle[0], (int)twiddle[2], (int)twiddle[4], (int)twiddle[6],
        (int)twiddle[0], (int)twiddle[2], (int)twiddle[4], (int)twiddle[6]);
    for (int i = 0; i < n; i += 16) {
        __m256i x = _mm256_loadu_si256((const __m256i *)(a + i));
        __m256i y = _mm256_loadu_si256((const __m256i *)(a + i + 8));
        __m256i even = _mm256_permute2x128_si256(x, y, 0x20);
        __m256i odd = _mm256_permute2x128_si256(x, y, 0x31);
        odd = mont_mul8(odd, w2);
        __m256i sums = reduce_sum(even, odd);
        __m256i diffs = reduce_diff(even, odd);
        __m256i left = _mm256_permute2x128_si256(sums, diffs, 0x20);
        __m256i right = _mm256_permute2x128_si256(sums, diffs, 0x31);
        right = mont_mul8(right, w1);
        _mm256_storeu_si256((__m256i *)(a + i), reduce_sum(left, right));
        _mm256_storeu_si256((__m256i *)(a + i + 8), reduce_diff(left, right));
    }
}

__attribute__((target("avx2"), always_inline))
static inline void inverse_block(u32 *a, int p, int half, int j) {
    const __m256i m2 = _mm256_set1_epi32((int)(2u * MOD));
    __m256i x = _mm256_loadu_si256((const __m256i *)(a + p + j));
    __m256i y = _mm256_loadu_si256((const __m256i *)(a + p + j + half));
    __m256i w = _mm256_loadu_si256((const __m256i *)(twiddle + j));
    y = mont_mul8(y, w);
    __m256i sum = _mm256_add_epi32(x, y);
    sum = _mm256_min_epu32(sum, _mm256_sub_epi32(sum, m2));
    __m256i diff = _mm256_sub_epi32(_mm256_add_epi32(x, m2), y);
    diff = _mm256_min_epu32(diff, _mm256_sub_epi32(diff, m2));
    _mm256_storeu_si256((__m256i *)(a + p + j), sum);
    _mm256_storeu_si256((__m256i *)(a + p + j + half), diff);
}

__attribute__((target("avx2")))
static void inverse_small(u32 *a, int n, int len) {
    if (len == 2) {
        const __m256i swap_index = _mm256_setr_epi32(1, 0, 3, 2, 5, 4, 7, 6);
        for (int i = 0; i < n; i += 8) {
            __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
            __m256i z = _mm256_permutevar8x32_epi32(v, swap_index);
            __m256i s = reduce_sum(v, z);
            __m256i d = reduce_diff(z, v);
            _mm256_storeu_si256((__m256i *)(a + i),
                                _mm256_blend_epi32(s, d, 0xaa));
        }
    } else if (len == 4) {
        __m256i w = _mm256_setr_epi32(
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[0], (int)twiddle[1]);
        for (int i = 0; i < n; i += 8) {
            __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
            __m256i x = _mm256_permute4x64_epi64(v, 0xa0);
            __m256i y = mont_mul8(_mm256_permute4x64_epi64(v, 0xf5), w);
            __m256i s = reduce_sum(x, y);
            __m256i d = reduce_diff(x, y);
            _mm256_storeu_si256((__m256i *)(a + i),
                                _mm256_blend_epi32(s, d, 0xcc));
        }
    } else {
        __m256i w = _mm256_setr_epi32(
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[2], (int)twiddle[3],
            (int)twiddle[0], (int)twiddle[1],
            (int)twiddle[2], (int)twiddle[3]);
        for (int i = 0; i < n; i += 8) {
            __m256i v = _mm256_loadu_si256((const __m256i *)(a + i));
            __m256i x = _mm256_permute2x128_si256(v, v, 0x00);
            __m256i y = mont_mul8(
                _mm256_permute2x128_si256(v, v, 0x11), w);
            __m256i s = reduce_sum(x, y);
            __m256i d = reduce_diff(x, y);
            _mm256_storeu_si256((__m256i *)(a + i),
                _mm256_permute2x128_si256(s, d, 0x20));
        }
    }
}

__attribute__((target("avx2")))
static void inverse(u32 *a, int n) {
    const u32 mod2 = 2u * MOD;
    int tail = n;
    while (tail >= 32) tail >>= 2;
    int len = 2;
    if (tail == 16) {
        u32 root4 = mod_pow(3, (MOD - 1u) / 4u);
        u32 step4 = mont_mul(mod_pow(root4, MOD - 2u), 932051910u);
        twiddle[0] = 301989884u;
        twiddle[1] = mont_mul(twiddle[0], step4);
        inverse_fused4(a, n);
        u32 root16 = mod_pow(3, (MOD - 1u) / 16u);
        u32 step16 = mont_mul(mod_pow(root16, MOD - 2u), 932051910u);
        twiddle[0] = 301989884u;
        for (int j = 1; j < 8; ++j) twiddle[j] = mont_mul(twiddle[j - 1], step16);
        inverse_fused16(a, n);
        len = 32;
    }
    for (; len <= tail; len <<= 1) {
        int half = len >> 1;
        u32 root = mod_pow(3, (MOD - 1u) / (u32)len);
        u32 step = mont_mul(mod_pow(root, MOD - 2u), 932051910u);
        twiddle[0] = 301989884u;
        for (int j = 1; j < half; ++j)
            twiddle[j] = mont_mul(twiddle[j - 1], step);
        if (len <= 8) {
            inverse_small(a, n, len);
            continue;
        }
        for (int i = 0; i < n; i += len) {
            int j = 0;
            for (; j + 8 <= half; j += 8) inverse_block(a, i, half, j);
            for (; j < half; ++j) {
                u32 w = twiddle[j];
                u32 x = a[i + j];
                u32 y = mont_mul(a[i + j + half], w);
                u32 sum = x + y;
                if (sum >= mod2) sum -= mod2;
                u32 diff = x + mod2 - y;
                if (diff >= mod2) diff -= mod2;
                a[i + j] = sum;
                a[i + j + half] = diff;
            }
        }
    }
    for (int outer = tail << 2; outer <= n; outer <<= 2) {
        int half = outer >> 1, quarter = outer >> 2;
        u32 root = mod_pow(3, (MOD - 1u) / (u32)outer);
        u32 step = mont_mul(mod_pow(root, MOD - 2u), 932051910u);
        twiddle[0] = 301989884u;
        for (int j = 1; j < half; ++j)
            twiddle[j] = mont_mul(twiddle[j - 1], step);
        for (int j = 0; j < quarter; ++j)
            twiddle_shoup[j] = twiddle[j << 1];
        for (int i = 0; i < n; i += outer)
            for (int j = 0; j < quarter; j += 8)
                inverse_radix4_block(a, i, quarter, j);
    }
    u32 inv_n = mod_pow((u32)n, MOD - 2u);
    __m256i vinv = _mm256_set1_epi32((int)inv_n);
    __m256i vmod = _mm256_set1_epi32((int)MOD);
    for (int i = 0; i < n; i += 8) {
        __m256i x = mont_mul8(
            _mm256_loadu_si256((const __m256i *)(a + i)), vinv);
        x = _mm256_min_epu32(x, _mm256_sub_epi32(x, vmod));
        _mm256_storeu_si256((__m256i *)(a + i), x);
    }
}

__attribute__((target("avx2")))
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
    static const u32 digit_mont[10] = {
        0u, 301989884u, 603979768u, 905969652u, 209715183u,
        511705067u, 813694951u, 117440482u, 419430366u, 721420250u
    };
    const int size = 1 << 21;
    const int need = n + m + 1;
    for (int i = 0; i <= n; ++i) fa[i] = digit_mont[a[i]];
    for (int i = n + 1; i < size; ++i) fa[i] = 0;
    for (int i = 0; i <= m; ++i) fb[i] = digit_mont[b[i]];
    for (int i = m + 1; i < size; ++i) fb[i] = 0;
    forward_both(fa, fb, size);
    for (int i = 0; i < size; i += 8) {
        __m256i x = _mm256_loadu_si256((const __m256i *)(fa + i));
        __m256i y = _mm256_loadu_si256((const __m256i *)(fb + i));
        _mm256_storeu_si256((__m256i *)(fa + i), mont_mul8(x, y));
    }
    inverse(fa, size);
    for (int i = 0; i < need; ++i) c[i] = fa[i];
}


enum { NN=1<<20, DECIMAL_DIGITS=1000000, INPUT_LIMBS=500000 };

static inline unsigned parse2(const char *p){return (unsigned)(p[0]-'0')*10u+(unsigned)(p[1]-'0');}

#ifdef LOCAL_TEST
static void solve_ntt(DuckInfo *di){
#else
static __attribute__((noreturn)) void solve_ntt(DuckInfo *di){
#endif
    unsigned mont_digit[101];
    for(int d=-50;d<=50;++d){unsigned x=d<0?(unsigned)(MOD+d):(unsigned)d;mont_digit[d+50]=(unsigned)((u64)x*301989884u%MOD);}
    const char *in=di->stdin_ptr;
    unsigned na=INPUT_LIMBS,nb=INPUT_LIMBS,carry=0;
    for(unsigned i=0;i<INPUT_LIMBS;++i){unsigned v=parse2(in+DECIMAL_DIGITS-2-2*i)+carry;carry=v>=50;if(carry)v-=100;fa[i]=mont_digit[(int)v+50];}
    if(carry)fa[na++]=mont_digit[51];
    carry=0;
    for(unsigned i=0;i<INPUT_LIMBS;++i){unsigned v=parse2(in+2*DECIMAL_DIGITS+1-2-2*i)+carry;carry=v>=50;if(carry)v-=100;fb[i]=mont_digit[(int)v+50];}
    if(carry)fb[nb++]=mont_digit[51];
    forward_both(fa,fb,NN);
    for(unsigned i=0;i<NN;i+=8){__m256i x=_mm256_loadu_si256((const __m256i*)(fa+i)),y=_mm256_loadu_si256((const __m256i*)(fb+i));_mm256_storeu_si256((__m256i*)(fa+i),mont_mul8(x,y));}
    inverse(fa,NN);
    long long cy=0;unsigned nc=na+nb;
    for(unsigned i=0;i<nc;++i){long long v=(long long)fa[i];if(v>(long long)MOD/2)v-=MOD;v+=cy;long long r=v%100;if(r<0)r+=100;cy=(v-r)/100;fa[i]=(unsigned)r;}
    while(cy){long long r=cy%100;if(r<0)r+=100;cy=(cy-r)/100;fa[nc++]=(unsigned)r;}
    while(nc>1&&!fa[nc-1])--nc;
    char *out=di->stdout_ptr,*p=out;p=duck_write_u64(p,fa[--nc]);
    while(nc){unsigned v=fa[--nc];*p++=(char)('0'+v/10);*p++=(char)('0'+v%10);}
    *p++='\n';di->stdout_size=p-out;
#ifdef LOCAL_TEST
    return;
#else
    duck_exit();
#endif
}

#ifndef LOCAL_TEST
extern "C" __attribute__((noreturn)) void __libc_start_main(void*,long argc,char **argv){solve_ntt(duck_info(argc,argv));}
int main(){}
#else
#include <stdio.h>
static char local_input_ntt[2000016],local_output_ntt[2000016];
int main(){DuckInfo di={};di.stdin_ptr=local_input_ntt;di.stdout_ptr=local_output_ntt;fread(local_input_ntt,1,sizeof(local_input_ntt),stdin);solve_ntt(&di);fwrite(local_output_ntt,1,di.stdout_size,stdout);}
#endif

评测结果