代码

#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));
    __m256i ws = _mm256_loadu_si256((const __m256i *)(twiddle_shoup + j));
    _mm256_storeu_si256((__m256i *)(a + p + j), sum);
    _mm256_storeu_si256((__m256i *)(a + p + j + half),
                        shoup_mul8(diff, w, ws));
}

__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));
}

__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]);
        __m256i ws = _mm256_setr_epi32(
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[0], (int)twiddle_shoup[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 = shoup_mul8(reduce_diff(x, y), w, ws);
            _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]);
        __m256i ws = _mm256_setr_epi32(
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[2], (int)twiddle_shoup[3],
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[2], (int)twiddle_shoup[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 = shoup_mul8(reduce_diff(x, y), w, ws);
            _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;
    for (int len = n; len >= 2; len >>= 1) {
        int half = len >> 1;
        u32 step = mod_pow(3, (MOD - 1u) / (u32)len);
        u32 wgen = 1;
        for (int j = 0; j < half; ++j) {
            twiddle[j] = wgen;
            twiddle_shoup[j] = (u32)(((u64)wgen << 32) / MOD);
            wgen = (u32)((u64)wgen * step % MOD);
        }
        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_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));
    __m256i ws = _mm256_loadu_si256((const __m256i *)(twiddle_shoup + j));
    y = shoup_mul8(y, w, ws);
    __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]);
        __m256i ws = _mm256_setr_epi32(
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[0], (int)twiddle_shoup[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 = shoup_mul8(_mm256_permute4x64_epi64(v, 0xf5), w, ws);
            __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]);
        __m256i ws = _mm256_setr_epi32(
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[2], (int)twiddle_shoup[3],
            (int)twiddle_shoup[0], (int)twiddle_shoup[1],
            (int)twiddle_shoup[2], (int)twiddle_shoup[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 = shoup_mul8(
                _mm256_permute2x128_si256(v, v, 0x11), w, ws);
            __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;
    for (int len = 2; len <= n; len <<= 1) {
        int half = len >> 1;
        u32 root = mod_pow(3, (MOD - 1u) / (u32)len);
        u32 step = mod_pow(root, MOD - 2u);
        u32 wgen = 1;
        for (int j = 0; j < half; ++j) {
            twiddle[j] = wgen;
            twiddle_shoup[j] = (u32)(((u64)wgen << 32) / MOD);
            wgen = (u32)((u64)wgen * step % MOD);
        }
        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;
            }
        }
    }
    u32 inv_n = mod_pow((u32)n, MOD - 2u);
    for (int i = 0; i < n; ++i) {
        u32 x = mont_mul(a[i], inv_n);
        if (x >= MOD) x -= MOD;
        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;

    for (int i = 0; i < an; ++i) {
        while ((unsigned)(*p - '0') > 9u) ++p;
        fa[i] = mont_mul((u32)(*p++ - '0'), 932051910u);
    }
    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] = mont_mul((u32)(*p++ - '0'), 932051910u);
    }
    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

评测结果