代码

#define PROFILE_PRODUCT
#pragma GCC optimize("O3,unroll-loops,omit-frame-pointer")
#pragma GCC target("arch=skylake")

#ifndef DUCK_FASTIO_H
#define DUCK_FASTIO_H

typedef unsigned long duck_u64;
typedef long duck_i64;

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

static __attribute__((always_inline)) inline DuckInfo *duck_info(long argc, char **argv) {
    char **p = argv + argc + 1;
    while (*p) ++p;
    duck_u64 *aux = (duck_u64 *)(p + 1);
    while (aux[0]) {
        if (aux[0] == 0x6b637564UL) return (DuckInfo *)aux[1];
        aux += 2;
    }
    return (DuckInfo *)0;
}

static __attribute__((always_inline)) inline duck_u64 duck_read_u64(const char **cursor) {
    const char *p = *cursor;
    while ((unsigned char)(*p - '0') > 9) ++p;
    duck_u64 value = 0;
    do {
        value = value * 10 + (unsigned char)(*p - '0');
        ++p;
    } while ((unsigned char)(*p - '0') <= 9);
    *cursor = p;
    return value;
}

static __attribute__((always_inline)) inline duck_i64 duck_read_i64(const char **cursor) {
    const char *p = *cursor;
    while (*p != '-' && (unsigned char)(*p - '0') > 9) ++p;
    int negative = *p == '-';
    p += negative;
    duck_u64 value = 0;
    do {
        value = value * 10 + (unsigned char)(*p - '0');
        ++p;
    } while ((unsigned char)(*p - '0') <= 9);
    *cursor = p;
    return negative ? -(duck_i64)value : (duck_i64)value;
}

static __attribute__((always_inline)) inline char *duck_write_u64(char *out, duck_u64 value) {
    char tmp[24];
    unsigned n = 0;
    do {
        tmp[n++] = (char)('0' + value % 10);
        value /= 10;
    } while (value);
    do *out++ = tmp[--n]; while (n);
    return out;
}

static __attribute__((always_inline)) inline char *duck_write_i64(char *out, duck_i64 value) {
    if (value < 0) {
        *out++ = '-';
        return duck_write_u64(out, (duck_u64)(-value));
    }
    return duck_write_u64(out, (duck_u64)value);
}

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

#endif

#include <immintrin.h>
#include <math.h>
#ifdef LOCAL_TEST
#include <stdio.h>
#endif

struct C { double r, i; };
enum { N=1<<19, H=N>>1, DIGITS=1000000, LIMBS=DIGITS/4, BASE=10000 };
alignas(32) static double zr[N], zi[N];
alignas(32) static unsigned digit4[10000];

static inline C mul(C a,C b){return {a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r};}
static inline C norm(C a){double s=1.5-.5*(a.r*a.r+a.i*a.i);return {a.r*s,a.i*s};}
static inline void vmul(__m256d ar,__m256d ai,__m256d br,__m256d bi,__m256d&rr,__m256d&ri){
    rr=_mm256_fmsub_pd(ar,br,_mm256_mul_pd(ai,bi));
    ri=_mm256_fmadd_pd(ar,bi,_mm256_mul_pd(ai,br));
}

static void radix2_forward(){
    const double ang=-2.0*3.141592653589793238462643383279502884/N;
    C step={cos(ang),sin(ang)},w={1,0};
    for(unsigned j=0;j<H;j+=4){
        if(j&&(j&63u)==0)w=norm(w);
        C w0=w,w1=mul(w0,step),w2=mul(w1,step),w3=mul(w2,step);w=mul(w3,step);
        __m256d wr=_mm256_setr_pd(w0.r,w1.r,w2.r,w3.r),wi=_mm256_setr_pd(w0.i,w1.i,w2.i,w3.i);
        __m256d ar=_mm256_load_pd(zr+j),ai=_mm256_load_pd(zi+j);
        __m256d br=_mm256_load_pd(zr+H+j),bi=_mm256_load_pd(zi+H+j),dr,di;
        _mm256_store_pd(zr+j,_mm256_add_pd(ar,br));_mm256_store_pd(zi+j,_mm256_add_pd(ai,bi));
        vmul(_mm256_sub_pd(ar,br),_mm256_sub_pd(ai,bi),wr,wi,dr,di);
        _mm256_store_pd(zr+H+j,dr);_mm256_store_pd(zi+H+j,di);
    }
}

static void radix2_inverse(){
    const double ang=-2.0*3.141592653589793238462643383279502884/N;
    C step={cos(ang),sin(ang)},w={1,0};
    const __m256d scale=_mm256_set1_pd(1.0/N);
    for(unsigned j=0;j<H;j+=4){
        if(j&&(j&63u)==0)w=norm(w);
        C w0=w,w1=mul(w0,step),w2=mul(w1,step),w3=mul(w2,step);w=mul(w3,step);
        __m256d wr=_mm256_setr_pd(w0.r,w1.r,w2.r,w3.r),wi=_mm256_setr_pd(-w0.i,-w1.i,-w2.i,-w3.i);
        __m256d ar=_mm256_load_pd(zr+j),ai=_mm256_load_pd(zi+j),br,bi;
        vmul(_mm256_load_pd(zr+H+j),_mm256_load_pd(zi+H+j),wr,wi,br,bi);
        _mm256_store_pd(zr+j,_mm256_mul_pd(_mm256_add_pd(ar,br),scale));
        _mm256_store_pd(zi+j,_mm256_mul_pd(_mm256_add_pd(ai,bi),scale));
        _mm256_store_pd(zr+H+j,_mm256_mul_pd(_mm256_sub_pd(ar,br),scale));
        _mm256_store_pd(zi+H+j,_mm256_mul_pd(_mm256_sub_pd(ai,bi),scale));
    }
}

template<unsigned M> static inline void forward_block(double *r,double *i){
    constexpr unsigned q=M/4;
    constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
    C step={cos(ang),sin(ang)},w={1,0};
    for(unsigned j=0;j<q;j+=4){
        if(j&&(j&63u)==0)w=norm(w);
        C w0=w,w1=mul(w0,step),w2=mul(w1,step),w3=mul(w2,step);w=mul(w3,step);
        __m256d wr=_mm256_setr_pd(w0.r,w1.r,w2.r,w3.r),wi=_mm256_setr_pd(w0.i,w1.i,w2.i,w3.i),w2r,w2i,w3r,w3i;
        vmul(wr,wi,wr,wi,w2r,w2i);vmul(w2r,w2i,wr,wi,w3r,w3i);
        __m256d ar=_mm256_load_pd(r+j),ai=_mm256_load_pd(i+j);
        __m256d br=_mm256_load_pd(r+q+j),bi=_mm256_load_pd(i+q+j);
        __m256d cr=_mm256_load_pd(r+2*q+j),ci=_mm256_load_pd(i+2*q+j);
        __m256d dr=_mm256_load_pd(r+3*q+j),di=_mm256_load_pd(i+3*q+j);
        __m256d t0r=_mm256_add_pd(ar,cr),t0i=_mm256_add_pd(ai,ci),t1r=_mm256_sub_pd(ar,cr),t1i=_mm256_sub_pd(ai,ci);
        __m256d t2r=_mm256_add_pd(br,dr),t2i=_mm256_add_pd(bi,di),vr=_mm256_sub_pd(br,dr),vi=_mm256_sub_pd(bi,di);
        __m256d o1r,o1i,o2r,o2i,o3r,o3i;
        vmul(_mm256_add_pd(t1r,vi),_mm256_sub_pd(t1i,vr),wr,wi,o1r,o1i);
        vmul(_mm256_sub_pd(t0r,t2r),_mm256_sub_pd(t0i,t2i),w2r,w2i,o2r,o2i);
        vmul(_mm256_sub_pd(t1r,vi),_mm256_add_pd(t1i,vr),w3r,w3i,o3r,o3i);
        _mm256_store_pd(r+j,_mm256_add_pd(t0r,t2r));_mm256_store_pd(i+j,_mm256_add_pd(t0i,t2i));
        _mm256_store_pd(r+q+j,o1r);_mm256_store_pd(i+q+j,o1i);
        _mm256_store_pd(r+2*q+j,o2r);_mm256_store_pd(i+2*q+j,o2i);
        _mm256_store_pd(r+3*q+j,o3r);_mm256_store_pd(i+3*q+j,o3i);
    }
}

template<unsigned M> static inline void inverse_block(double *r,double *i){
    constexpr unsigned q=M/4;
    constexpr double ang=-2.0*3.141592653589793238462643383279502884/M;
    C step={cos(ang),sin(ang)},w={1,0};
    for(unsigned j=0;j<q;j+=4){
        if(j&&(j&63u)==0)w=norm(w);
        C w0=w,w1=mul(w0,step),w2=mul(w1,step),w3=mul(w2,step);w=mul(w3,step);
        __m256d wr=_mm256_setr_pd(w0.r,w1.r,w2.r,w3.r),wi=_mm256_setr_pd(-w0.i,-w1.i,-w2.i,-w3.i),w2r,w2i,w3r,w3i;
        vmul(wr,wi,wr,wi,w2r,w2i);vmul(w2r,w2i,wr,wi,w3r,w3i);
        __m256d ar=_mm256_load_pd(r+j),ai=_mm256_load_pd(i+j),br,bi,cr,ci,dr,di;
        vmul(_mm256_load_pd(r+q+j),_mm256_load_pd(i+q+j),wr,wi,br,bi);
        vmul(_mm256_load_pd(r+2*q+j),_mm256_load_pd(i+2*q+j),w2r,w2i,cr,ci);
        vmul(_mm256_load_pd(r+3*q+j),_mm256_load_pd(i+3*q+j),w3r,w3i,dr,di);
        __m256d t0r=_mm256_add_pd(ar,cr),t0i=_mm256_add_pd(ai,ci),t1r=_mm256_sub_pd(ar,cr),t1i=_mm256_sub_pd(ai,ci);
        __m256d t2r=_mm256_add_pd(br,dr),t2i=_mm256_add_pd(bi,di),vr=_mm256_sub_pd(br,dr),vi=_mm256_sub_pd(bi,di);
        _mm256_store_pd(r+j,_mm256_add_pd(t0r,t2r));_mm256_store_pd(i+j,_mm256_add_pd(t0i,t2i));
        _mm256_store_pd(r+q+j,_mm256_sub_pd(t1r,vi));_mm256_store_pd(i+q+j,_mm256_add_pd(t1i,vr));
        _mm256_store_pd(r+2*q+j,_mm256_sub_pd(t0r,t2r));_mm256_store_pd(i+2*q+j,_mm256_sub_pd(t0i,t2i));
        _mm256_store_pd(r+3*q+j,_mm256_add_pd(t1r,vi));_mm256_store_pd(i+3*q+j,_mm256_sub_pd(t1i,vr));
    }
}

static inline void forward4(double*r,double*i){
    double ar=r[0],ai=i[0],br=r[1],bi=i[1],cr=r[2],ci=i[2],dr=r[3],di=i[3];
    double t0r=ar+cr,t0i=ai+ci,t1r=ar-cr,t1i=ai-ci,t2r=br+dr,t2i=bi+di,vr=br-dr,vi=bi-di;
    r[0]=t0r+t2r;i[0]=t0i+t2i;r[1]=t1r+vi;i[1]=t1i-vr;r[2]=t0r-t2r;i[2]=t0i-t2i;r[3]=t1r-vi;i[3]=t1i+vr;
}
static inline void inverse4(double*r,double*i){
    double ar=r[0],ai=i[0],br=r[1],bi=i[1],cr=r[2],ci=i[2],dr=r[3],di=i[3];
    double t0r=ar+cr,t0i=ai+ci,t1r=ar-cr,t1i=ai-ci,t2r=br+dr,t2i=bi+di,vr=br-dr,vi=bi-di;
    r[0]=t0r+t2r;i[0]=t0i+t2i;r[1]=t1r-vi;i[1]=t1i+vr;r[2]=t0r-t2r;i[2]=t0i-t2i;r[3]=t1r+vi;i[3]=t1i-vr;
}

template<unsigned M> struct RecFFT{
    static void forward(double*r,double*i){forward_block<M>(r,i);constexpr unsigned q=M/4;RecFFT<q>::forward(r,i);RecFFT<q>::forward(r+q,i+q);RecFFT<q>::forward(r+2*q,i+2*q);RecFFT<q>::forward(r+3*q,i+3*q);}
    static void inverse(double*r,double*i){constexpr unsigned q=M/4;RecFFT<q>::inverse(r,i);RecFFT<q>::inverse(r+q,i+q);RecFFT<q>::inverse(r+2*q,i+2*q);RecFFT<q>::inverse(r+3*q,i+3*q);inverse_block<M>(r,i);}
};
template<> struct RecFFT<256>{
    static void forward(double*r,double*i){forward_block<256>(r,i);for(unsigned p=0;p<256;p+=64)forward_block<64>(r+p,i+p);for(unsigned p=0;p<256;p+=16)forward_block<16>(r+p,i+p);for(unsigned p=0;p<256;p+=4)forward4(r+p,i+p);}
    static void inverse(double*r,double*i){for(unsigned p=0;p<256;p+=4)inverse4(r+p,i+p);for(unsigned p=0;p<256;p+=16)inverse_block<16>(r+p,i+p);for(unsigned p=0;p<256;p+=64)inverse_block<64>(r+p,i+p);inverse_block<256>(r,i);}
};
static void forward(){radix2_forward();RecFFT<H>::forward(zr,zi);RecFFT<H>::forward(zr+H,zi+H);}
static void inverse(){RecFFT<H>::inverse(zr,zi);RecFFT<H>::inverse(zr+H,zi+H);radix2_inverse();}

static inline unsigned parse4(const char *p){return (unsigned)(p[0]-'0')*1000u+(unsigned)(p[1]-'0')*100u+(unsigned)(p[2]-'0')*10u+(unsigned)(p[3]-'0');}
static inline void spectral_pair(unsigned p,unsigned q){
    if(p==q){zr[p]*=zi[p];zi[p]=0;return;}
    double ar=zr[p],ai=zi[p],br=zr[q],bi=-zi[q];
    double xr=(ar+br)*.5,xi=(ai+bi)*.5,yr=(ai-bi)*.5,yi=(br-ar)*.5;
    double vr=xr*yr-xi*yi,vi=xr*yi+xi*yr;zr[p]=vr;zi[p]=vi;zr[q]=vr;zi[q]=-vi;
}
static void spectral_range(unsigned p,unsigned q){
    const __m256d half=_mm256_set1_pd(.5),sign=_mm256_set1_pd(-0.0);
    while(p+6<q){
        __m256d ar=_mm256_load_pd(zr+p),ai=_mm256_load_pd(zi+p);
        __m256d br=_mm256_permute4x64_pd(_mm256_loadu_pd(zr+q-3),0x1b);
        __m256d bi=_mm256_xor_pd(_mm256_permute4x64_pd(_mm256_loadu_pd(zi+q-3),0x1b),sign);
        __m256d xr=_mm256_mul_pd(_mm256_add_pd(ar,br),half),xi=_mm256_mul_pd(_mm256_add_pd(ai,bi),half);
        __m256d yr=_mm256_mul_pd(_mm256_sub_pd(ai,bi),half),yi=_mm256_mul_pd(_mm256_sub_pd(br,ar),half),vr,vi;
        vmul(xr,xi,yr,yi,vr,vi);_mm256_store_pd(zr+p,vr);_mm256_store_pd(zi+p,vi);
        _mm256_storeu_pd(zr+q-3,_mm256_permute4x64_pd(vr,0x1b));
        _mm256_storeu_pd(zi+q-3,_mm256_xor_pd(_mm256_permute4x64_pd(vi,0x1b),sign));
        p+=4;q-=4;
    }
    while(p<=q){spectral_pair(p,q);++p;--q;}
}

#ifdef LOCAL_TEST
static void solve(DuckInfo *di){
#else
static __attribute__((noreturn)) void solve(DuckInfo *di){
#endif
    const char *in=di->stdin_ptr;
    for(unsigned j=0;j<LIMBS;++j){zr[j]=parse4(in+DIGITS-4-4*j);zi[j]=parse4(in+2*DIGITS+1-4-4*j);}
#ifdef PROFILE_PARSE
    di->stdout_size=0;duck_exit();
#endif
    forward();
#ifdef PROFILE_FORWARD
    di->stdout_size=0;duck_exit();
#endif
    spectral_pair(0,0);for(unsigned h=1;h<H;h<<=2)spectral_range(h,4*h-1);spectral_range(H,N-1);
#ifdef PROFILE_PRODUCT
    di->stdout_size=0;duck_exit();
#endif
    inverse();
#ifdef PROFILE_INVERSE
    di->stdout_size=0;duck_exit();
#endif
    for(unsigned v=0;v<10000;++v)digit4[v]=(unsigned)('0'+v/1000)|((unsigned)('0'+v/100%10)<<8)|((unsigned)('0'+v/10%10)<<16)|((unsigned)('0'+v%10)<<24);
    char *out=di->stdout_ptr;unsigned long long carry=0;
    for(unsigned j=0;j<2*LIMBS;++j){unsigned long long v=(unsigned long long)(zr[j]+.5)+carry;carry=v/BASE;((unsigned*)out)[2*LIMBS-1-j]=digit4[v-carry*BASE];}
    unsigned skip=0;while(skip<3&&out[skip]=='0')++skip;
    if(skip)for(unsigned j=0;j<2*DIGITS-skip;++j)out[j]=out[j+skip];
    out[2*DIGITS-skip]='\n';di->stdout_size=2*DIGITS+1-skip;
#ifdef LOCAL_TEST
    return;
#else
    duck_exit();
#endif
}

#ifndef LOCAL_TEST
extern "C" __attribute__((noreturn)) void __libc_start_main(void*,long argc,char **argv){solve(duck_info(argc,argv));}
int main(){}
#else
static char local_input[2000016],local_output[2000016];
int main(){DuckInfo di={};di.stdin_ptr=local_input;di.stdout_ptr=local_output;fread(local_input,1,sizeof(local_input),stdin);solve(&di);fwrite(local_output,1,di.stdout_size,stdout);}
#endif

评测结果