代码

// 1004a 高精度乘法2 - single 64-bit NTT (Goldilocks prime 2^64-2^32+1), base 1e8
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>

typedef uint64_t u64;
typedef int64_t i64;
typedef __uint128_t u128;
typedef __int128_t i128;

static const u64 P = 18446744069414584321ULL; // 2^64 - 2^32 + 1
static const i128 P128 = (i128)P;
static const u64 W4096 = 0xF2C35199959DFCB6ULL; // primitive 4096-th root of unity
static const u64 WINV   = 0x9AF01E431FBD6EA0ULL; // W4096^{-1}

// (a*b) mod P, a,b in [0,P)  (branchless)
static inline u64 reduce(u64 a, u64 b){
    u128 x = (u128)a * b;
    u64 x0 = (u64)x, x1 = (u64)(x >> 64);
    u64 h = x1 >> 32, l = x1 & 0xFFFFFFFFULL;
    u64 s = x0 - h;   u64 b1 = (x0 < h);
    u64 t = s - l;    u64 b2 = (s < l);
    u64 r = t + (l << 32);  u64 c = (r < t);
    i64 kk = (i64)c - (i64)b1 - (i64)b2;
    i128 res = (i128)r + (i128)kk * 0xFFFFFFFFLL;
    u64 lo = (u64)res;
    i64 hi = (i64)(res >> 64);
    u64 lo2 = lo + P;
    u64 lo_nonneg = (hi < 0) ? lo2 : lo;
    i64 hi2 = (hi < 0) ? 0 : hi;
    u64 ge = (u64)(hi2 > 0) | (u64)(lo_nonneg >= P);
    u64 lo3 = lo_nonneg - P;
    return ge ? lo3 : lo_nonneg;
}
// (u+v) mod P, u,v in [0,P)
static inline u64 addmod(u64 u, u64 v){
    u64 s = u + v; u64 t = s - P;
    u64 m = -(u64)((s < u) | (s >= P));
    return (t & m) | (s & ~m);
}
// (u-v) mod P, u,v in [0,P)
static inline u64 submod(u64 u, u64 v){
    u64 d = u - v; u64 m = -(u64)(u < v);
    return (d & ~m) | ((d + P) & m);
}

#define N 4096
static u64 A[N], B[N];
static u64 roots[N], iroots[N]; // roots[k] = w^k, iroots[k] = w^{-k}

static void build_roots(void){
    roots[0] = 1;
    for(int i=1;i<N;i++) roots[i] = reduce(roots[i-1], W4096);
    iroots[0] = 1;
    for(int i=1;i<N;i++) iroots[i] = reduce(iroots[i-1], WINV);
}

// forward DIF NTT, in-place, output bit-reversed
static void ntt_fwd(u64* a){
    for(int len=N; len>1; len>>=1){
        int half = len>>1;
        int step = N / len;
        const u64* w = roots;
        for(int i=0;i<N;i+=len){
            u64* base = a + i;
            int j = 0;
            for(; j+3 < half; j += 4){
                u64 u0=base[j], v0=base[j+half];
                u64 u1=base[j+1], v1=base[j+1+half];
                u64 u2=base[j+2], v2=base[j+2+half];
                u64 u3=base[j+3], v3=base[j+3+half];
                base[j]      = addmod(u0,v0);
                base[j+half]   = reduce(submod(u0,v0), w[(j)*step]);
                base[j+1]      = addmod(u1,v1);
                base[j+1+half] = reduce(submod(u1,v1), w[(j+1)*step]);
                base[j+2]      = addmod(u2,v2);
                base[j+2+half] = reduce(submod(u2,v2), w[(j+2)*step]);
                base[j+3]      = addmod(u3,v3);
                base[j+3+half] = reduce(submod(u3,v3), w[(j+3)*step]);
            }
            for(; j < half; j++){
                u64 u=base[j], v=base[j+half];
                base[j] = addmod(u,v);
                base[j+half] = reduce(submod(u,v), w[j*step]);
            }
        }
    }
}

// inverse DIT NTT, in-place, input bit-reversed, output natural
static void ntt_inv(u64* a){
    for(int len=2; len<=N; len<<=1){
        int half = len>>1;
        int step = N / len;
        const u64* w = iroots;
        for(int i=0;i<N;i+=len){
            u64* base = a + i;
            int j = 0;
            for(; j+3 < half; j += 4){
                u64 u0=base[j], v0=reduce(base[j+half], w[(j)*step]);
                u64 u1=base[j+1], v1=reduce(base[j+1+half], w[(j+1)*step]);
                u64 u2=base[j+2], v2=reduce(base[j+2+half], w[(j+2)*step]);
                u64 u3=base[j+3], v3=reduce(base[j+3+half], w[(j+3)*step]);
                base[j]      = addmod(u0,v0);
                base[j+half] = submod(u0,v0);
                base[j+1]      = addmod(u1,v1);
                base[j+1+half] = submod(u1,v1);
                base[j+2]      = addmod(u2,v2);
                base[j+2+half] = submod(u2,v2);
                base[j+3]      = addmod(u3,v3);
                base[j+3+half] = submod(u3,v3);
            }
            for(; j < half; j++){
                u64 u=base[j], v=reduce(base[j+half], w[j*step]);
                base[j] = addmod(u,v);
                base[j+half] = submod(u,v);
            }
        }
    }
    // scale by N^{-1} mod P
    u128 r=1, b=N, ee=P-2;
    while(ee){ if(ee&1) r=(u128)r*b%P; b=(u128)b*b%P; ee>>=1; }
    u64 ninv=(u64)r;
    for(int i=0;i<N;i++) a[i] = reduce(a[i], ninv);
}

static char ibuf[30000];
static int ilen = 0;

static char obuf[21000];
static int olen = 0;

static char tab3[1000][4];
static void build_tab3(void){
    for(int i=0;i<1000;i++){
        int v=i;
        tab3[i][2] = '0' + v%10; v/=10;
        tab3[i][1] = '0' + v%10; v/=10;
        tab3[i][0] = '0' + v%10;
    }
}

int main(){
    build_roots();
    build_tab3();

    ilen = fread(ibuf, 1, sizeof(ibuf), stdin);

    int pa = 0;
    while(pa < ilen && (ibuf[pa]==' '||ibuf[pa]=='\n'||ibuf[pa]=='\r'||ibuf[pa]=='\t')) pa++;
    int a_start = pa;
    while(pa < ilen && ibuf[pa]!=' '&&ibuf[pa]!='\n'&&ibuf[pa]!='\r'&&ibuf[pa]!='\t') pa++;
    int a_end = pa;
    while(pa < ilen && (ibuf[pa]==' '||ibuf[pa]=='\n'||ibuf[pa]=='\r'||ibuf[pa]=='\t')) pa++;
    int b_start = pa;
    while(pa < ilen && ibuf[pa]!=' '&&ibuf[pa]!='\n'&&ibuf[pa]!='\r'&&ibuf[pa]!='\t') pa++;
    int b_end = pa;

    int sa = a_start; while(sa < a_end-1 && ibuf[sa]=='0') sa++;
    int sb = b_start; while(sb < b_end-1 && ibuf[sb]=='0') sb++;

    int na=0, nb=0;
    {
        int pos = a_end;
        while(pos > sa){
            int start = pos-8; if(start < sa) start = sa;
            u64 v=0;
            for(int i=start;i<pos;i++) v = v*10 + (ibuf[i]-'0');
            A[na++] = v; pos = start;
        }
    }
    {
        int pos = b_end;
        while(pos > sb){
            int start = pos-8; if(start < sb) start = sb;
            u64 v=0;
            for(int i=start;i<pos;i++) v = v*10 + (ibuf[i]-'0');
            B[nb++] = v; pos = start;
        }
    }

    ntt_fwd(A);
    ntt_fwd(B);
    for(int i=0;i<N;i++) A[i] = reduce(A[i], B[i]);
    ntt_inv(A);

    u64 carry = 0;
    int outlen = na + nb - 1;
    for(int i=0;i<outlen;i++){
        u64 c = A[i] + carry;
        A[i] = c % 100000000ULL;
        carry = c / 100000000ULL;
    }
    while(carry){
        A[outlen++] = carry % 100000000ULL;
        carry /= 100000000ULL;
    }
    int hi = outlen-1;
    while(hi>0 && A[hi]==0) hi--;

    {
        u64 v = A[hi];
        char tmp[9];
        int t=0;
        do { tmp[t++] = '0' + v%10; v/=10; } while(v);
        while(t>0) obuf[olen++] = tmp[--t];
    }
    for(int i=hi-1;i>=0;i--){
        u64 v = A[i];
        u64 g2 = v / 1000000;
        u64 g1 = (v/1000) % 1000;
        u64 g0 = v % 1000;
        obuf[olen++] = '0' + g2/10;
        obuf[olen++] = '0' + g2%10;
        const char* p = tab3[g1];
        obuf[olen++]=p[0]; obuf[olen++]=p[1]; obuf[olen++]=p[2];
        p = tab3[g0];
        obuf[olen++]=p[0]; obuf[olen++]=p[1]; obuf[olen++]=p[2];
    }
    obuf[olen++] = '\n';
    fwrite(obuf, 1, olen, stdout);
    return 0;
}

评测结果