代码

#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 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};
    x[0] = {t0.r + t2.r, t0.i + t2.i};
    x[1] = {t1.r + id.r, t1.i + id.i};
    x[2] = {t0.r - t2.r, t0.i - t2.i};
    x[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};
    x[4] = {t0.r + t2.r, t0.i + t2.i};
    x[5] = {t1.r + id.r, t1.i + id.i};
    x[6] = {t0.r - t2.r, t0.i - t2.i};
    x[7] = {t1.r - id.r, t1.i - id.i};
    const double s = 0.7071067811865475244;
#define COMBINE(K, Q) do { \
        Complex e = x[K], q = (Q); \
        x[K] = {e.r + q.r, e.i + q.i}; \
        x[(K) + 4] = {e.r - q.r, e.i - q.i}; \
    } while (0)
    COMBINE(0, x[4]);
    if (!inverse) {
        COMBINE(1, (Complex{s * (x[5].r - x[5].i), s * (x[5].r + x[5].i)}));
        COMBINE(2, (Complex{-x[6].i, x[6].r}));
        COMBINE(3, (Complex{-s * (x[7].r + x[7].i), s * (x[7].r - x[7].i)}));
    } else {
        COMBINE(1, (Complex{s * (x[5].r + x[5].i), s * (x[5].i - x[5].r)}));
        COMBINE(2, (Complex{x[6].i, -x[6].r}));
        COMBINE(3, (Complex{s * (x[7].i - x[7].r), -s * (x[7].r + x[7].i)}));
    }
#undef COMBINE
}

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

static void fft_inverse(Complex *a) {
    for (unsigned length = 8; 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];
            }
        }
    }
    const double scale = 1.0 / FFT_N;
    for (unsigned i = 0; i < FFT_N; ++i) {
        a[i].r *= scale;
        a[i].i *= scale;
    }
}

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 + 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];
        *p++ = (char)('0' + value / 10000);
        *p++ = (char)('0' + value / 1000 % 10);
        *p++ = (char)('0' + value / 100 % 10);
        *p++ = (char)('0' + value / 10 % 10);
        *p++ = (char)('0' + value % 10);
    }
    *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() {}

评测结果