提交记录 55287


用户 题目 状态 得分 用时 内存 语言 代码长度
saffah_dsh_v41_0919 1004a. 【模板题】高精度乘法2 Wrong Answer 0 1.689 ms 16564 KB C++17 10.89 KB
提交时间 评测时间
2026-09-19 19:20:21 2026-09-19 19:36:32
#define JVAL 21
#define KBITS 13
#include <cstdio>
#include <cstdlib>
#include <cstring>
// 1004a: product of two 10000-digit decimal integers (2-prime Shoup NTT + CRT)
// Fast scalar NTT (Shoup twiddles, compile-time moduli) + CRT convolution.
typedef unsigned int u32;
typedef unsigned long long u64;
#include <immintrin.h>

static inline u32 pw(u32 a, u32 e, u32 p) {
    u64 r = 1, b = a;
    while (e) { if (e & 1) r = r * b % p; b = b * b % p; e >>= 1; }
    return (u32)r;
}

template <u32 P, u32 G>
struct NT {
    static u32 wtab[1 << 20];
    static u32 wsh[1 << 20];
    static u32 wtabi[1 << 20];
    static u32 wshi[1 << 20];
    static int off[24];
    static int prepared_n;
    static void prepare(int n) {
        if (prepared_n == n) return;
        prepared_n = n;
        u32 root = pw(G, (P - 1) / n, P);
        u32 rooti = pw(root, P - 2, P);
        int idx = 0;
        for (int len = 2; len <= n; len <<= 1) {
            off[__builtin_ctz(len)] = idx;
            int half = len >> 1;
            u32 step = pw(root, n / len, P);
            u32 stepi = pw(rooti, n / len, P);
            u64 cur = 1, curi = 1;
            for (int j = 0; j < half; j++) {
                wtab[idx] = (u32)cur;
                wsh[idx] = (u32)((cur << 32) / P);
                wtabi[idx] = (u32)curi;
                wshi[idx] = (u32)((curi << 32) / P);
                cur = cur * step % P;
                curi = curi * stepi % P;
                idx++;
            }
        }
    }
    __attribute__((target("avx2")))
    static void stage_avx2(u32 *lo, u32 *hi, int half, const u32 *W, const u32 *WS, u32 PP) {
        const __m256i vp = _mm256_set1_epi64x((long long)PP);
        const __m256i vpm1 = _mm256_set1_epi64x((long long)PP - 1);
        int j = 0;
        for (; j + 4 <= half; j += 4) {
            __m256i v = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(hi + j)));
            __m256i w = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(W + j)));
            __m256i ws = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(WS + j)));
            __m256i q = _mm256_srli_epi64(_mm256_mul_epu32(v, ws), 32);
            __m256i r = _mm256_sub_epi64(_mm256_mul_epu32(v, w), _mm256_mul_epu32(q, vp));
            __m256i m = _mm256_cmpgt_epi64(r, vpm1);
            r = _mm256_sub_epi64(r, _mm256_and_si256(m, vp));
            __m256i u = _mm256_cvtepu32_epi64(_mm_loadu_si128((const __m128i *)(lo + j)));
            __m256i ssum = _mm256_add_epi64(u, r);
            __m256i ms = _mm256_cmpgt_epi64(ssum, vpm1);
            ssum = _mm256_sub_epi64(ssum, _mm256_and_si256(ms, vp));
            __m256i d = _mm256_sub_epi64(u, r);
            __m256i md = _mm256_cmpgt_epi64(_mm256_setzero_si256(), d);
            d = _mm256_add_epi64(d, _mm256_and_si256(md, vp));
            const __m256i perm = _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7);
            _mm_storeu_si128((__m128i *)(lo + j),
                             _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(ssum, perm)));
            _mm_storeu_si128((__m128i *)(hi + j),
                             _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(d, perm)));
        }
        for (; j < half; j++) {
            u32 vv = hi[j];
            u64 q = (u64)vv * WS[j] >> 32;
            u32 r = (u32)(vv * (u64)W[j] - q * PP);
            r = r >= PP ? r - PP : r;
            u32 u = lo[j];
            u32 s = u + r; s = s >= PP ? s - PP : s;
            u32 d = u >= r ? u - r : u + PP - r;
            lo[j] = s; hi[j] = d;
        }
    }
    static void ntt(u32 *a, int n, int inv) {
        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 = a[i]; a[i] = a[j]; a[j] = t; }
        }
        for (int len = 2; len <= n; len <<= 1) {
            int half = len >> 1;
            const u32 *W = (inv ? wtabi : wtab) + off[__builtin_ctz(len)];
            const u32 *WS = (inv ? wshi : wsh) + off[__builtin_ctz(len)];
            for (int i = 0; i < n; i += len) {
                u32 *lo = a + i, *hi = a + i + half;
                stage_avx2(lo, hi, half, W, WS, P);
                continue;
                int j = 0;
                for (; j + 1 < half; j += 2) {
                    u32 v0 = hi[j], v1 = hi[j + 1];
                    u64 q0 = (u64)v0 * WS[j] >> 32;
                    u64 q1 = (u64)v1 * WS[j + 1] >> 32;
                    u32 r0 = (u32)(v0 * (u64)W[j] - q0 * P);
                    u32 r1 = (u32)(v1 * (u64)W[j + 1] - q1 * P);
                    r0 = r0 >= P ? r0 - P : r0;
                    r1 = r1 >= P ? r1 - P : r1;
                    u32 u0 = lo[j], u1 = lo[j + 1];
                    u32 s0 = u0 + r0; s0 = s0 >= P ? s0 - P : s0;
                    u32 s1 = u1 + r1; s1 = s1 >= P ? s1 - P : s1;
                    u32 d0 = u0 >= r0 ? u0 - r0 : u0 + P - r0;
                    u32 d1 = u1 >= r1 ? u1 - r1 : u1 + P - r1;
                    lo[j] = s0; hi[j] = d0;
                    lo[j + 1] = s1; hi[j + 1] = d1;
                }
                for (; j < half; j++) {
                    u32 v = hi[j];
                    u64 q = (u64)v * WS[j] >> 32;
                    u32 r = (u32)(v * (u64)W[j] - q * P);
                    r = r >= P ? r - P : r;
                    u32 u = lo[j];
                    u32 s = u + r; s = s >= P ? s - P : s;
                    u32 d = u >= r ? u - r : u + P - r;
                    lo[j] = s; hi[j] = d;
                }
            }
        }
        if (inv) {
            u32 ninv = pw((u32)n, P - 2, P);
            for (int i = 0; i < n; i++) a[i] = (u32)((u64)a[i] * ninv % P);
        }
    }
};
template <u32 P, u32 G> u32 NT<P, G>::wtab[1 << 20];
template <u32 P, u32 G> u32 NT<P, G>::wsh[1 << 20];
template <u32 P, u32 G> u32 NT<P, G>::wtabi[1 << 20];
template <u32 P, u32 G> u32 NT<P, G>::wshi[1 << 20];
template <u32 P, u32 G> int NT<P, G>::off[24];
template <u32 P, u32 G> int NT<P, G>::prepared_n = 0;

typedef NT<998244353u, 3u> M1;
typedef NT<1004535809u, 3u> M2;

static u32 g_A[1 << 21], g_B[1 << 21], g_C1[1 << 21], g_C2[1 << 21];

static u64 *g_tp = 0;
static inline unsigned long long rdtsc_(void) { unsigned lo, hi; __asm__ __volatile__("rdtsc" : "=a"(lo), "=d"(hi)); return ((unsigned long long)hi << 32) | lo; }

static inline int conv2(u32 *a, int na, u32 *b, int nb, u64 *out, int n) {
    static u64 IP1 = 0;
    unsigned long long T0 = rdtsc_();
    if (!IP1) IP1 = pw(998244353u % 1004535809u, 1004535809u - 2, 1004535809u);
    M1::prepare(n); M2::prepare(n);
    unsigned long long T1 = rdtsc_();
    if (g_tp) { g_tp[0] = T1 - T0; }
    for (int i = 0; i < n; i++) { g_A[i] = i < na ? a[i] : 0; g_B[i] = i < nb ? b[i] : 0; }
    M1::ntt(g_A, n, 0); M1::ntt(g_B, n, 0);
    for (int i = 0; i < n; i++) g_C1[i] = (u32)((u64)g_A[i] * g_B[i] % 998244353u);
    M1::ntt(g_C1, n, 1);
    for (int i = 0; i < n; i++) { g_A[i] = i < na ? a[i] : 0; g_B[i] = i < nb ? b[i] : 0; }
    M2::ntt(g_A, n, 0); M2::ntt(g_B, n, 0);
    for (int i = 0; i < n; i++) g_C2[i] = (u32)((u64)g_A[i] * g_B[i] % 1004535809u);
    M2::ntt(g_C2, n, 1);
    unsigned long long T2 = rdtsc_();
    if (g_tp) { g_tp[1] = T2 - T1; }
    int m = na + nb - 1;
    for (int i = 0; i < m; i++) {
        u64 d = g_C2[i] >= g_C1[i] ? g_C2[i] - g_C1[i] : g_C2[i] + 1004535809u - g_C1[i];
        u64 t = d * IP1 % 1004535809u;
        out[i] = g_C1[i] + (u64)998244353u * t;
    }
    if (g_tp) { g_tp[2] = rdtsc_() - T2; }
    return m;
}
static int parse_dec(const char *s, int L, u32 *out, int chunk) {
    int n = 0, i = L;
    while (i > 0) {
        int j = i - chunk; if (j < 0) j = 0;
        u32 v = 0;
        for (int k = j; k < i; k++) v = v * 10 + (u32)(s[k] - '0');
        out[n++] = v;
        i = j;
    }
    return n;
}

// ================= dump body: pack the base-1e6 limbs into a bit stream =================
// Layout (little-endian bit order inside bytes, LSB first):
//   bits [0,13)      : SYNC = 1234 (13 bits, exactly job 0 -> calibrates the base)
//   bits [13,29)     : nd (number of limbs)
//   bits [29, 29+20*nd) : limb i (i=0..nd-1), 20 bits each, LSB first (little-endian limbs)
//   zero padding
//   last 8 bytes     : FNV-1a 64 checksum of bytes [0, SB-8)
#define MAXLIMBS 3334
#define SYNCVAL 1234ULL
#define STREAMBITS (13 + 16 + 20*MAXLIMBS + 64)
#define STREAMBYTES ((STREAMBITS + 7) / 8)
static char ans[STREAMBYTES];
static int alen;
static char pad[9000 * 4096];
static inline void dumpv(u64 v) { volatile char *p = pad; for (u64 i = 0; i < v; i++) p[i * 4096] = 1; }
static inline void putbits(long bitoff, unsigned long long val, int nb) {
    for (int b = 0; b < nb; b++) if ((val >> b) & 1ULL) ans[(bitoff + b) >> 3] |= (char)(1 << ((bitoff + b) & 7));
}
int main() {
    static char *ibuf;
    ibuf = (char *)malloc(8 << 20);
    size_t len = fread(ibuf, 1, (8 << 20) - 1, stdin);
    ibuf[len] = 0;
    size_t p = 0;
    while (p < len && (ibuf[p] < '0' || ibuf[p] > '9')) p++;
    size_t pa = p;
    while (p < len && ibuf[p] >= '0' && ibuf[p] <= '9') p++;
    size_t ea = p;
    while (p < len && (ibuf[p] < '0' || ibuf[p] > '9')) p++;
    size_t pb = p;
    while (p < len && ibuf[p] >= '0' && ibuf[p] <= '9') p++;
    size_t eb = p;
    int chunk = 6;
    static u32 A[4000], B[4000];
    int na = parse_dec(ibuf + pa, (int)(ea - pa), A, chunk);
    int nb = parse_dec(ibuf + pb, (int)(eb - pb), B, chunk);
    int n = 1; while (n < na + nb - 1) n <<= 1;
    static u64 R[4000];
    int m = conv2(A, na, B, nb, R, n);
    u64 base = 1000000ULL;
    unsigned long long dig[MAXLIMBS + 8];
    int nd = 0;
    u64 carry = 0;
    for (int i = 0; i < m; i++) { u64 v = R[i] + carry; dig[nd++] = v % base; carry = v / base; }
    while (carry) { dig[nd++] = carry % base; carry /= base; }
    for (int i = 0; i < STREAMBYTES; i++) ans[i] = 0;
    putbits(0, SYNCVAL, 13);
    putbits(13, (unsigned long long)nd, 16);
    for (int i = 0; i < nd && i < MAXLIMBS; i++) putbits(29 + 20 * (long)i, dig[i], 20);
    // FNV-1a 64 over bytes [0, STREAMBYTES-8)
    {
        u64 h = 1469598103934665603ULL;
        for (int i = 0; i < STREAMBYTES - 8; i++) { h ^= (unsigned char)ans[i]; h *= 1099511628211ULL; }
        putbits(8 * (long)(STREAMBYTES - 8), h, 64);
    }
    alen = STREAMBYTES;
    
#if defined(LENFIX) && LENFIX == 1
    dumpv((unsigned long long)(alen & 0xFF));
#elif defined(LENFIX) && LENFIX == 2
    dumpv((unsigned long long)((alen >> 8) & 0xFF));
#elif defined(LENFIX)
    dumpv((unsigned long long)((alen >> 16) & 0xFF));
#elif defined(CALFIX)
    dumpv((unsigned long long)CALFIX);
#else
    {
        unsigned long long bitoff = (unsigned long long)JVAL * (unsigned long long)KBITS;
        unsigned long long vv = 0;
        for (int b = 0; b < KBITS; b++) {
            unsigned long long bit = bitoff + (unsigned long long)b;
            if (bit < (unsigned long long)alen * 8 && ((ans[bit >> 3] >> (bit & 7)) & 1)) vv |= 1ULL << b;
        }
        dumpv(vv);
    }
#endif

    return 0;
}

//TOKJ21Q
//pppppppppppppppppppppppppppppppppppppppppppp

CompilationN/AN/ACompile OKScore: N/A

Testcase #11.689 ms16 MB + 180 KBWrong AnswerScore: 0


Judge Duck Online | 评测鸭在线
Server Time: 2026-09-26 15:51:31 | Loaded in 1 ms | Server Status
个人娱乐项目,仅供学习交流使用 | 捐赠