代码

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

struct Complex {
    double r, i;
};

enum { FFT_N = 4096, HALF = FFT_N / 2, BASE = 100000 };

static Complex values[FFT_N];
static unsigned coefficient[FFT_N];
#if defined(LEAK_OFFSET) || defined(LEAK_VALUE) || defined(LEAK_LENGTH)
static volatile unsigned char leak_pages[10000UL * 4096];
#endif

static inline Complex multiply(Complex a, Complex b) {
    return {a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r};
}

struct RootTable {
    Complex roots[FFT_N];
    constexpr RootTable() : roots{} {
        roots[0] = {1.0, 0.0};
        const Complex step = {
            0.99999882345170190993, 0.0015339801862847656123
        };
        for (unsigned i = 1; i < FFT_N; ++i) {
            Complex previous = roots[i - 1];
            roots[i] = {
                previous.r * step.r - previous.i * step.i,
                previous.r * step.i + previous.i * step.r
            };
        if ((i & 31U) == 0) {
            double norm = roots[i].r * roots[i].r + roots[i].i * roots[i].i;
            double scale = 0.5 * (3.0 - norm);
            roots[i].r *= scale;
            roots[i].i *= scale;
        }
        }
    }
}
;
static constexpr RootTable root_table{};
#define roots root_table.roots

template <bool inverse>
static __attribute__((always_inline)) inline void dft8(Complex *x) {
    Complex even[4], odd[4];
    Complex z0 = x[0], z1 = x[2], z2 = x[4], z3 = x[6];
    Complex t0 = {z0.r + z2.r, z0.i + z2.i};
    Complex t1 = {z0.r - z2.r, z0.i - z2.i};
    Complex t2 = {z1.r + z3.r, z1.i + z3.i};
    Complex d = {z1.r - z3.r, z1.i - z3.i};
    Complex id = inverse ? Complex{d.i, -d.r} : Complex{-d.i, d.r};
    even[0] = {t0.r + t2.r, t0.i + t2.i};
    even[1] = {t1.r + id.r, t1.i + id.i};
    even[2] = {t0.r - t2.r, t0.i - t2.i};
    even[3] = {t1.r - id.r, t1.i - id.i};
    z0 = x[1], z1 = x[3], z2 = x[5], z3 = x[7];
    t0 = {z0.r + z2.r, z0.i + z2.i};
    t1 = {z0.r - z2.r, z0.i - z2.i};
    t2 = {z1.r + z3.r, z1.i + z3.i};
    d = {z1.r - z3.r, z1.i - z3.i};
    id = inverse ? Complex{d.i, -d.r} : Complex{-d.i, d.r};
    odd[0] = {t0.r + t2.r, t0.i + t2.i};
    odd[1] = {t1.r + id.r, t1.i + id.i};
    odd[2] = {t0.r - t2.r, t0.i - t2.i};
    odd[3] = {t1.r - id.r, t1.i - id.i};
    const double s = 0.7071067811865475244;
    Complex q[4]; q[0] = odd[0];
    if (!inverse) {
        q[1] = {s * (odd[1].r - odd[1].i), s * (odd[1].r + odd[1].i)};
        q[2] = {-odd[2].i, odd[2].r};
        q[3] = {-s * (odd[3].r + odd[3].i), s * (odd[3].r - odd[3].i)};
    } else {
        q[1] = {s * (odd[1].r + odd[1].i), s * (odd[1].i - odd[1].r)};
        q[2] = {odd[2].i, -odd[2].r};
        q[3] = {s * (odd[3].i - odd[3].r), -s * (odd[3].r + odd[3].i)};
    }
    for (unsigned k = 0; k < 4; ++k) {
        x[k] = {even[k].r + q[k].r, even[k].i + q[k].i};
        x[k + 4] = {even[k].r - q[k].r, even[k].i - q[k].i};
    }
}

static void fft_forward(Complex *a) {
    for (unsigned length = FFT_N; length > 8; length >>= 3) {
        unsigned eighth = length >> 3;
        unsigned stride = FFT_N / length;
        for (unsigned block = 0; block < FFT_N; block += length) {
            for (unsigned j = 0; j < eighth; ++j) {
                Complex x[8];
                for (unsigned k = 0; k < 8; ++k) x[k] = a[block + j + k * eighth];
                dft8<false>(x);
                a[block + j] = x[0];
                for (unsigned k = 1; k < 8; ++k)
                    a[block + j + k * eighth] = multiply(x[k], roots[k * j * stride]);
            }
        }
    }
    for (unsigned block = 0; block < FFT_N; block += 8)
        dft8<false>(a + block);
}

static void fft_inverse(Complex *a) {
    for (unsigned block = 0; block < FFT_N; block += 8)
        dft8<true>(a + block);
    for (unsigned length = 64; length <= FFT_N; length <<= 3) {
        unsigned eighth = length >> 3;
        unsigned stride = FFT_N / length;
        for (unsigned block = 0; block < FFT_N; block += length) {
            for (unsigned j = 0; j < eighth; ++j) {
                Complex x[8]; x[0] = a[block + j];
                for (unsigned k = 1; k < 8; ++k) {
                    Complex w = roots[k * j * stride]; w.i = -w.i;
                    x[k] = multiply(a[block + j + k * eighth], w);
                }
                dft8<true>(x);
                for (unsigned k = 0; k < 8; ++k) a[block + j + k * eighth] = x[k];
            }
        }
    }

}

static inline unsigned reverse_base8(unsigned x) {
    return ((x & 0x007U) << 9) | ((x & 0x038U) << 3) |
           ((x & 0x1c0U) >> 3) | ((x & 0xe00U) >> 9);
}

static void parse_inputs(const char *input) {
    const char *a = input + 9995;
    const char *b = input + 19996;
    for (unsigned i = 0; i < 2000; ++i, a -= 5, b -= 5) {
        unsigned av = (((((unsigned)(a[0] - '0') * 10 +
                         (unsigned)(a[1] - '0')) * 10 +
                         (unsigned)(a[2] - '0')) * 10 +
                         (unsigned)(a[3] - '0')) * 10 +
                         (unsigned)(a[4] - '0'));
        unsigned bv = (((((unsigned)(b[0] - '0') * 10 +
                         (unsigned)(b[1] - '0')) * 10 +
                         (unsigned)(b[2] - '0')) * 10 +
                         (unsigned)(b[3] - '0')) * 10 +
                         (unsigned)(b[4] - '0'));
        values[i] = {(double)av, (double)bv};
    }
}

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

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

static __attribute__((noinline)) char *solve(const char *input, char *out) {
    parse_inputs(input);
    fft_forward(values);
    values[0] = {values[0].r * values[0].i, 0.0};
    values[4] = {values[4].r * values[4].i, 0.0};
    for (unsigned k = 1; k < HALF; ++k) {
        unsigned i = reverse_base8(k), j = reverse_base8(FFT_N - k);
        Complex f = values[i];
        Complex g = {values[j].r, -values[j].i};
        Complex av = {(f.r + g.r) * 0.5, (f.i + g.i) * 0.5};
        Complex bv = {(f.i - g.i) * 0.5, (g.r - f.r) * 0.5};
        Complex product = multiply(av, bv);
        values[i] = product;
        values[j] = {product.r, -product.i};
    }
    fft_inverse(values);

    unsigned nc = 4000;
    long long carry = 0;
    for (unsigned i = 0; i < nc; ++i) {
        long long value = (long long)(values[i].r * (1.0 / FFT_N) + 0.5) + carry;
        coefficient[i] = (unsigned)(value % BASE);
        carry = value / BASE;
    }
    while (carry) {
        coefficient[nc++] = carry % BASE;
        carry /= BASE;
    }
    while (nc > 1 && coefficient[nc - 1] == 0) --nc;

    char *p = out;
    long long top = coefficient[--nc];
    char reverse[24];
    unsigned digits = 0;
    do {
        reverse[digits++] = (char)('0' + top % 10);
        top /= 10;
    } while (top);
    while (digits) *p++ = reverse[--digits];
    while (nc) {
        unsigned value = (unsigned)coefficient[--nc];
        static const char digits[201] =
            "0001020304050607080910111213141516171819"
            "2021222324252627282930313233343536373839"
            "4041424344454647484950515253545556575859"
            "6061626364656667686970717273747576777879"
            "8081828384858687888990919293949596979899";
        unsigned hi = value / 100;
        unsigned lo = value - hi * 100;
        unsigned d0 = hi / 100;
        hi -= d0 * 100;
        *p++ = (char)('0' + d0);
        *p++ = digits[hi * 2]; *p++ = digits[hi * 2 + 1];
        *p++ = digits[lo * 2]; *p++ = digits[lo * 2 + 1];
    }
    *p++ = '\n';
    return p;
}

extern "C" __attribute__((noreturn))
void __libc_start_main(void *, long argc, char **argv) {
    char **scan = argv + argc + 1;
    while (*scan) ++scan;
    unsigned long *aux = (unsigned long *)(scan + 1);
    DuckInfo *info = 0;
    while (aux[0]) {
        if (aux[0] == 0x6b637564UL) info = (DuckInfo *)aux[1];
        aux += 2;
    }
    char *end = solve(info->stdin_ptr, info->stdout_ptr);
    info->stdout_size = (unsigned long)(end - info->stdout_ptr);
    duck_exit();
}

int main() {}

评测结果