代码

#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 << 18 };
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);
    }
}

static const char digit_pairs[201] =
    "00010203040506070809"
    "10111213141516171819"
    "20212223242526272829"
    "30313233343536373839"
    "40414243444546474849"
    "50515253545556575859"
    "60616263646566676869"
    "70717273747576777879"
    "80818283848586878889"
    "90919293949596979899";

static __attribute__((always_inline)) inline char *small(char *p, u32 x) {
    if (x >= 1000u) {
        u32 a = x / 100u;
        *(unsigned short *)p = *(const unsigned short *)(digit_pairs + 2 * a);
        *(unsigned short *)(p + 2) =
            *(const unsigned short *)(digit_pairs + 2 * (x - a * 100u));
        return p + 4;
    }
    if (x >= 100u) {
        u32 a = x / 100u;
        *p++ = (char)('0' + a);
        *(unsigned short *)p =
            *(const unsigned short *)(digit_pairs + 2 * (x - a * 100u));
        return p + 2;
    }
    if (x >= 10u) {
        *(unsigned short *)p = *(const unsigned short *)(digit_pairs + 2 * x);
        return p + 2;
    }
    *p = (char)('0' + x);
    return p + 1;
}

static __attribute__((always_inline)) inline char *write_u32(char *p, u32 x) {
    if (x < 10000u) return small(p, x);
    u32 hi = x / 10000u;
    u32 lo = x - hi * 10000u;
    p = small(p, hi);
    u32 a = lo / 100u;
    *(unsigned short *)p = *(const unsigned short *)(digit_pairs + 2 * a);
    *(unsigned short *)(p + 2) =
        *(const unsigned short *)(digit_pairs + 2 * (lo - a * 100u));
    return p + 4;
}

static __attribute__((always_inline)) inline u32 read_u32(const char **pp) {
    const char *p = *pp;
    while ((unsigned)(*p - '0') > 9u) ++p;
    u32 x = 0;
    do { x = x * 10u + (u32)(*p++ - '0'); }
    while ((unsigned)(*p - '0') <= 9u);
    *pp = p;
    return x;
}

__attribute__((target("avx2")))
static void run(DuckInfo *info) {
    const char *p = info->stdin_ptr;
    int dn = (int)read_u32(&p), dm = (int)read_u32(&p);
    int an = dn + 1, bn = dm + 1;
    int need = an + bn - 1, size = 1;
    while (size < need) size <<= 1;

    static const u32 digit_mont[10] = {
        0u, 301989884u, 603979768u, 905969652u, 209715183u,
        511705067u, 813694951u, 117440482u, 419430366u, 721420250u
    };
    for (int i = 0; i < an; ++i) {
        while ((unsigned)(*p - '0') > 9u) ++p;
        fa[i] = digit_mont[(unsigned)(*p++ - '0')];
    }
    for (int i = an; i < size; ++i) fa[i] = 0;
    for (int i = 0; i < bn; ++i) {
        while ((unsigned)(*p - '0') > 9u) ++p;
        fb[i] = digit_mont[(unsigned)(*p++ - '0')];
    }
    for (int i = bn; i < size; ++i) fb[i] = 0;

    forward_both(fa, fb, size);
    for (int i = 0; i < size; i += 8) {
        __m256i a = _mm256_loadu_si256((const __m256i *)(fa + i));
        __m256i b = _mm256_loadu_si256((const __m256i *)(fb + i));
        _mm256_storeu_si256((__m256i *)(fa + i), mont_mul8(a, b));
    }
    inverse(fa, size);

    char *out = info->stdout_ptr;
#ifndef BENCH
    for (int i = 0; i < need; ++i) {
        out = write_u32(out, fa[i]);
        *out++ = i + 1 == need ? '\n' : ' ';
    }
#endif
    info->stdout_size = (duck_u64)(out - info->stdout_ptr);
}

#ifndef LOCAL
__attribute__((noreturn))
void __libc_start_main(void *unused, long argc, char **argv) {
    (void)unused;
    DuckInfo *info = duck_info(argc, argv);
    run(info);
    duck_exit();
}
int main(void) {}
#else
extern long read(int, void *, unsigned long);
extern long write(int, const void *, unsigned long);
static char local_in[4000000], local_out[4000000];
int main(void) {
    long n = read(0, local_in, sizeof(local_in));
    DuckInfo info = {0};
    info.stdin_ptr = local_in;
    info.stdin_size = (duck_u64)n;
    info.stdout_ptr = local_out;
    run(&info);
    write(1, local_out, info.stdout_size);
    return 0;
}
#endif

评测结果