代码

#define PROFILE_INVERSE
#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 C z[N];

static inline __m256d cmul(__m256d a,__m256d b){
    __m256d ar=_mm256_movedup_pd(a),ai=_mm256_permute_pd(a,15),bs=_mm256_permute_pd(b,5);
    return _mm256_fmaddsub_pd(ar,b,_mm256_mul_pd(ai,bs));
}
static inline __m256d conjv(__m256d a){
    return _mm256_xor_pd(a,_mm256_castsi256_pd(_mm256_setr_epi64x(0,0x8000000000000000ULL,0,0x8000000000000000ULL)));
}
static inline __m256d minus_i(__m256d a){
    return _mm256_xor_pd(_mm256_permute_pd(a,5),_mm256_castsi256_pd(_mm256_setr_epi64x(0,0x8000000000000000ULL,0,0x8000000000000000ULL)));
}
static inline __m256d plus_i(__m256d a){
    return _mm256_xor_pd(_mm256_permute_pd(a,5),_mm256_castsi256_pd(_mm256_setr_epi64x(0x8000000000000000ULL,0,0x8000000000000000ULL,0)));
}
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 void radix2_forward(){
    const double ang=-2.0*3.141592653589793238462643383279502884/N;
    C step={cos(ang),sin(ang)};
    for(unsigned j=0;j<H;j+=2){
        C w0,w1;
        if((j&63u)==0){double x=ang*j;w0={cos(x),sin(x)};}else w0={0,0};
        /* The previous pair is deliberately reconstructed only between resets. */
        static C w={1,0};
        if((j&63u)==0)w=w0; else w=mul(w,mul(step,step));
        w0=w;w1=mul(w0,step);
        __m256d a=_mm256_load_pd((double*)(z+j)),b=_mm256_load_pd((double*)(z+H+j));
        _mm256_store_pd((double*)(z+j),_mm256_add_pd(a,b));
        _mm256_store_pd((double*)(z+H+j),cmul(_mm256_sub_pd(a,b),_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i)));
    }
}

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+=2){
        if((j&63u)==0){double x=ang*j;w={cos(x),sin(x)};}
        C w0=w,w1=mul(w,step);w=mul(w1,step);
        __m256d a=_mm256_load_pd((double*)(z+j));
        __m256d b=cmul(_mm256_load_pd((double*)(z+H+j)),_mm256_setr_pd(w0.r,-w0.i,w1.r,-w1.i));
        _mm256_store_pd((double*)(z+j),_mm256_mul_pd(_mm256_add_pd(a,b),scale));
        _mm256_store_pd((double*)(z+H+j),_mm256_mul_pd(_mm256_sub_pd(a,b),scale));
    }
}

static void forward(){
    radix2_forward();
    for(unsigned len=N/2;len>=8;len>>=2){
        unsigned q=len>>2;
        double ang=-2.0*3.141592653589793238462643383279502884/len;
        C step={cos(ang),sin(ang)};
        for(unsigned p=0;p<N;p+=len){
            C w={1,0};
            for(unsigned j=0;j<q;j+=2){
                if((j&63u)==0){double x=ang*j;w={cos(x),sin(x)};}
                C w0=w,w1=mul(w,step);w=mul(w1,step);
                __m256d rw=_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i);
                __m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
                __m256d a=_mm256_load_pd((double*)(z+p+j));
                __m256d b=_mm256_load_pd((double*)(z+p+q+j));
                __m256d c=_mm256_load_pd((double*)(z+p+2*q+j));
                __m256d d=_mm256_load_pd((double*)(z+p+3*q+j));
                __m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c);
                __m256d t2=_mm256_add_pd(b,d),t3=minus_i(_mm256_sub_pd(b,d));
                _mm256_store_pd((double*)(z+p+j),_mm256_add_pd(t0,t2));
                _mm256_store_pd((double*)(z+p+q+j),cmul(_mm256_add_pd(t1,t3),rw));
                _mm256_store_pd((double*)(z+p+2*q+j),cmul(_mm256_sub_pd(t0,t2),rw2));
                _mm256_store_pd((double*)(z+p+3*q+j),cmul(_mm256_sub_pd(t1,t3),rw3));
            }
        }
    }
    for(unsigned p=0;p<N;p+=4){
        C a=z[p],b=z[p+1],c=z[p+2],d=z[p+3];
        C t0={a.r+c.r,a.i+c.i},t1={a.r-c.r,a.i-c.i};
        C t2={b.r+d.r,b.i+d.i},x={b.r-d.r,b.i-d.i},t3={x.i,-x.r};
        z[p]={t0.r+t2.r,t0.i+t2.i};z[p+1]={t1.r+t3.r,t1.i+t3.i};
        z[p+2]={t0.r-t2.r,t0.i-t2.i};z[p+3]={t1.r-t3.r,t1.i-t3.i};
    }
}

static void inverse(){
    for(unsigned p=0;p<N;p+=4){
        C a=z[p],b=z[p+1],c=z[p+2],d=z[p+3];
        C t0={a.r+c.r,a.i+c.i},t1={a.r-c.r,a.i-c.i};
        C t2={b.r+d.r,b.i+d.i},x={b.r-d.r,b.i-d.i},t3={-x.i,x.r};
        z[p]={t0.r+t2.r,t0.i+t2.i};z[p+1]={t1.r+t3.r,t1.i+t3.i};
        z[p+2]={t0.r-t2.r,t0.i-t2.i};z[p+3]={t1.r-t3.r,t1.i-t3.i};
    }
    for(unsigned len=16;len<=N/2;len<<=2){
        unsigned q=len>>2;
        double ang=-2.0*3.141592653589793238462643383279502884/len;
        C step={cos(ang),sin(ang)};
        for(unsigned p=0;p<N;p+=len){
            C w={1,0};
            for(unsigned j=0;j<q;j+=2){
                if((j&63u)==0){double x=ang*j;w={cos(x),sin(x)};}
                C w0=w,w1=mul(w,step);w=mul(w1,step);
                __m256d rw=conjv(_mm256_setr_pd(w0.r,w0.i,w1.r,w1.i));
                __m256d rw2=cmul(rw,rw),rw3=cmul(rw2,rw);
                __m256d a=_mm256_load_pd((double*)(z+p+j));
                __m256d b=cmul(_mm256_load_pd((double*)(z+p+q+j)),rw);
                __m256d c=cmul(_mm256_load_pd((double*)(z+p+2*q+j)),rw2);
                __m256d d=cmul(_mm256_load_pd((double*)(z+p+3*q+j)),rw3);
                __m256d t0=_mm256_add_pd(a,c),t1=_mm256_sub_pd(a,c);
                __m256d t2=_mm256_add_pd(b,d),t3=plus_i(_mm256_sub_pd(b,d));
                _mm256_store_pd((double*)(z+p+j),_mm256_add_pd(t0,t2));
                _mm256_store_pd((double*)(z+p+q+j),_mm256_add_pd(t1,t3));
                _mm256_store_pd((double*)(z+p+2*q+j),_mm256_sub_pd(t0,t2));
                _mm256_store_pd((double*)(z+p+3*q+j),_mm256_sub_pd(t1,t3));
            }
        }
    }
    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 unsigned logical_index(unsigned p){
    unsigned k=(p>>18)&1u;
    for(unsigned d=0;d<9;++d)k|=((p>>(2*d))&3u)<<(17-2*d);
    return k;
}
static inline unsigned physical_index(unsigned k){
    unsigned p=(k&1u)<<18;
    for(unsigned d=0;d<9;++d)p|=((k>>(17-2*d))&3u)<<(2*d);
    return p;
}

#ifdef LOCAL_TEST
static void solve(DuckInfo *di){
#else
static __attribute__((noreturn)) void solve(DuckInfo *di){
#endif
    const char *in=di->stdin_ptr;
#ifdef LOCAL_IMPULSE
    z[1].r=1;
#elif defined(LOCAL_PACKIMPULSE)
    z[1]={1,1};
#elif defined(LOCAL_ONES)
    for(unsigned i=0;i<LIMBS;++i)z[i]={1,1};
#elif defined(LOCAL_REALONES)
    for(unsigned i=0;i<LIMBS;++i)z[i].r=1;
#else
    for(unsigned i=0;i<LIMBS;++i){
        z[i].r=parse4(in+DIGITS-4-4*i);
        z[i].i=parse4(in+2*DIGITS+1-4-4*i);
    }
#endif
#ifdef PROFILE_PARSE
    di->stdout_size=0;duck_exit();
#endif
    forward();
#ifdef PROFILE_FORWARD
    di->stdout_size=0;duck_exit();
#endif
#ifdef INSPECT_FORWARD
    for(unsigned i=0;i<32;++i)fprintf(stderr,"%u %.9f %.9f\n",i,z[i].r,z[i].i);
    for(unsigned i=1;i<32;++i){unsigned p=i*(N/32);fprintf(stderr,"P %u %.9f %.9f\n",p,z[p].r,z[p].i);}
    return;
#endif
#ifdef LOCAL_REALONES
    for(unsigned p=0;p<N;++p)z[p]=mul(z[p],z[p]);
#else
#ifndef ROUNDTRIP
    for(unsigned p=0;p<N;++p){
        unsigned k=logical_index(p),q=physical_index((-k)&(N-1));
        if(p>q)continue;
        if(p==q){z[p]={z[p].r*z[p].i,0};continue;}
        C a=z[p],b={z[q].r,-z[q].i};
        C av={(a.r+b.r)*.5,(a.i+b.i)*.5};
        C bv={(a.i-b.i)*.5,(b.r-a.r)*.5};
        C v=mul(av,bv);z[p]=v;z[q]={v.r,-v.i};
    }
#endif
#endif
#ifdef PROFILE_PRODUCT
    di->stdout_size=0;duck_exit();
#endif
    inverse();
#ifdef PROFILE_INVERSE
    di->stdout_size=0;duck_exit();
#endif
#ifdef LOCAL_DEBUG
    for(unsigned i=0;i<16;++i)fprintf(stderr,"%u %.6f %.6f\n",i,z[i].r,z[i].i);
#endif
    unsigned long long carry=0;unsigned nc=2*LIMBS;
    for(unsigned i=0;i<nc;++i){unsigned long long v=(unsigned long long)(z[i].r+.5)+carry;carry=v/BASE;((unsigned long long*)z)[i]=v-carry*BASE;}
    while(carry){((unsigned long long*)z)[nc++]=carry%BASE;carry/=BASE;}
    while(nc>1&&!((unsigned long long*)z)[nc-1])--nc;
    char *out=di->stdout_ptr,*p=out;p=duck_write_u64(p,((unsigned long long*)z)[--nc]);
    while(nc){unsigned v=(unsigned)((unsigned long long*)z)[--nc];*p++='0'+v/1000;*p++='0'+v/100%10;*p++='0'+v/10%10;*p++='0'+v%10;}
    *p++='\n';di->stdout_size=p-out;
#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

评测结果