代码

// 1004: 1e6-digit multiply - double FFT, split-layout (no shuffles), real-packing.
// radix-2 DIF/DIT, base 1e4, N=2^19. DuckInfo direct IO.
#include <sys/auxv.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <immintrin.h>
#pragma GCC target("avx2,fma")
#pragma GCC optimize("O3,unroll-loops")

typedef uint32_t u32;
typedef uint64_t u64;

#define N (1<<19)
#define LOGBITS 19

static double re[N] __attribute__((aligned(32)));
static double im[N] __attribute__((aligned(32)));
static double tw_re[N] __attribute__((aligned(32)));
static double tw_im[N] __attribute__((aligned(32)));
static u32 outlimbs[N+16];
static char tab4[10000][4];

struct DuckInfo {
  uint64_t abi_version; const char *stdin_ptr; uint64_t stdin_size;
  char *stdout_ptr; uint64_t stdout_limit; uint64_t stdout_size;
  char *stderr_ptr; uint64_t stderr_limit; uint64_t stderr_size;
  const char *IB_ptr; uint64_t IB_limit; char *OB_ptr; uint64_t OB_limit;
  uint64_t tsc_frequency;
} __attribute__((packed));

static void build_tab4(void){
    for(int i=0;i<10000;i++){ int v=i; tab4[i][3]='0'+v%10; v/=10; tab4[i][2]='0'+v%10; v/=10; tab4[i][1]='0'+v%10; v/=10; tab4[i][0]='0'+v%10; }
}

// staged twiddle table: for stage len, twiddles w[j*step] for j=0..half-1, contiguous.
static void gen_tw(void){
    int off=0;
    for(int len=N; len>1; len>>=1){
        int half=len>>1, step=N/len;
        double d = 2.0*M_PI*(double)step/(double)N;
        double cd = cos(d), sd = sin(d);
        double c=1.0, s=0.0;
        for(int j=0;j<half;j++){
            if((j & 63)==0){
                double a = d*(double)j;
                c = cos(a); s = sin(a);
            }
            tw_re[off+j]=c; tw_im[off+j]=s;
            double nc = c*cd - s*sd; double ns = s*cd + c*sd;
            c=nc; s=ns;
        }
        off+=half;
    }
}

// forward DIF radix-2, split layout, natural -> bit-reversed
static void fft_fwd(void){
    int off=0;
    for(int len=N; len>1; len>>=1){
        int half=len>>1;
        for(int i=0;i<N;i+=len){
            int j=0;
            for(; j+3<half; j+=4){
                __m256d ur=_mm256_load_pd(&re[i+j]);
                __m256d ui=_mm256_load_pd(&im[i+j]);
                __m256d vr=_mm256_load_pd(&re[i+j+half]);
                __m256d vi=_mm256_load_pd(&im[i+j+half]);
                __m256d wr=_mm256_load_pd(&tw_re[off+j]);
                __m256d wi=_mm256_load_pd(&tw_im[off+j]);
                __m256d sr=_mm256_add_pd(ur,vr);
                __m256d si=_mm256_add_pd(ui,vi);
                __m256d tr=_mm256_sub_pd(ur,vr);
                __m256d ti=_mm256_sub_pd(ui,vi);
                __m256d dr=_mm256_fmsub_pd(tr,wr,_mm256_mul_pd(ti,wi));
                __m256d di=_mm256_fmadd_pd(tr,wi,_mm256_mul_pd(ti,wr));
                _mm256_store_pd(&re[i+j],sr);
                _mm256_store_pd(&im[i+j],si);
                _mm256_store_pd(&re[i+j+half],dr);
                _mm256_store_pd(&im[i+j+half],di);
            }
            for(; j<half; j++){
                double ur=re[i+j], ui=im[i+j], vr=re[i+j+half], vi=im[i+j+half];
                double wr=tw_re[off+j], wi=tw_im[off+j];
                re[i+j]=ur+vr; im[i+j]=ui+vi;
                re[i+j+half]=(ur-vr)*wr-(ui-vi)*wi;
                im[i+j+half]=(ur-vr)*wi+(ui-vi)*wr;
            }
        }
        off+=half;
    }
}

// inverse DIT radix-2, split layout, bit-reversed -> natural
static void fft_inv(void){
    int off=N-1;
    for(int len=2; len<=N; len<<=1){
        int half=len>>1;
        off-=half;
        for(int i=0;i<N;i+=len){
            int j=0;
            for(; j+3<half; j+=4){
                __m256d ur=_mm256_load_pd(&re[i+j]);
                __m256d ui=_mm256_load_pd(&im[i+j]);
                __m256d vr=_mm256_load_pd(&re[i+j+half]);
                __m256d vi=_mm256_load_pd(&im[i+j+half]);
                __m256d wr=_mm256_load_pd(&tw_re[off+j]);
                __m256d wi=_mm256_load_pd(&tw_im[off+j]);
                wi=_mm256_sub_pd(_mm256_setzero_pd(),wi);
                __m256d dr=_mm256_fmsub_pd(vr,wr,_mm256_mul_pd(vi,wi));
                __m256d di=_mm256_fmadd_pd(vr,wi,_mm256_mul_pd(vi,wr));
                _mm256_store_pd(&re[i+j],_mm256_add_pd(ur,dr));
                _mm256_store_pd(&im[i+j],_mm256_add_pd(ui,di));
                _mm256_store_pd(&re[i+j+half],_mm256_sub_pd(ur,dr));
                _mm256_store_pd(&im[i+j+half],_mm256_sub_pd(ui,di));
            }
            for(; j<half; j++){
                double ur=re[i+j], ui=im[i+j], vr=re[i+j+half], vi=im[i+j+half];
                double wr=tw_re[off+j], wi=-tw_im[off+j];
                double dr=vr*wr-vi*wi, di=vr*wi+vi*wr;
                re[i+j]=ur+dr; im[i+j]=ui+di;
                re[i+j+half]=ur-dr; im[i+j+half]=ui-di;
            }
        }
    }
    __m256d invn=_mm256_set1_pd(1.0/N);
    for(int i=0;i<N;i+=4){
        _mm256_store_pd(&re[i],_mm256_mul_pd(_mm256_load_pd(&re[i]),invn));
        _mm256_store_pd(&im[i],_mm256_mul_pd(_mm256_load_pd(&im[i]),invn));
    }
}

// in-place binary bit-reversal for 19 bits (re and im together)
static void bitrev(void){
    int j=0;
    for(int i=1;i<N-1;i++){
        int k=N>>1;
        while(k<=j){ j-=k; k>>=1; }
        j+=k;
        if(i<j){
            double t=re[i]; re[i]=re[j]; re[j]=t;
            t=im[i]; im[i]=im[j]; im[j]=t;
        }
    }
}

#ifdef LOCAL_TEST
extern uintptr_t jd_getauxval(uintptr_t);
int jd_main(){
    DuckInfo* di=(DuckInfo*)jd_getauxval(0x6b637564ull);
#else
int main(){
    DuckInfo* di=(DuckInfo*)getauxval(0x6b637564ull);
#endif
    build_tab4();
    gen_tw();

    const char* ibuf=di->stdin_ptr;
    int ilen=(int)di->stdin_size;
    int pa=0;
    while(pa<ilen && (ibuf[pa]=='\n'||ibuf[pa]=='\r'||ibuf[pa]==' '||ibuf[pa]=='\t')) pa++;
    int a_start=pa; while(pa<ilen && ibuf[pa]>='0'&&ibuf[pa]<='9') pa++; int a_end=pa;
    while(pa<ilen && (ibuf[pa]=='\n'||ibuf[pa]=='\r'||ibuf[pa]==' '||ibuf[pa]=='\t')) pa++;
    int b_start=pa; while(pa<ilen && ibuf[pa]>='0'&&ibuf[pa]<='9') pa++; int b_end=pa;

    int sa=a_start; while(sa<a_end-1 && ibuf[sa]=='0') sa++;
    int sb=b_start; while(sb<b_end-1 && ibuf[sb]=='0') sb++;

    int na=0, nb=0;
    { int pos=a_end; while(pos>=sa+4){ pos-=4; u32 v=(u32)(ibuf[pos]-'0')*1000+(u32)(ibuf[pos+1]-'0')*100+(u32)(ibuf[pos+2]-'0')*10+(u32)(ibuf[pos+3]-'0'); re[na++]=v; } if(pos>sa){ u32 v=0; for(int i=sa;i<pos;i++) v=v*10+(u32)(ibuf[i]-'0'); re[na++]=v; } }
    { int pos=b_end; while(pos>=sb+4){ pos-=4; u32 v=(u32)(ibuf[pos]-'0')*1000+(u32)(ibuf[pos+1]-'0')*100+(u32)(ibuf[pos+2]-'0')*10+(u32)(ibuf[pos+3]-'0'); im[nb++]=v; } if(pos>sb){ u32 v=0; for(int i=sb;i<pos;i++) v=v*10+(u32)(ibuf[i]-'0'); im[nb++]=v; } }
    for(int i=na;i<N;i++) re[i]=0.0;
    for(int i=nb;i<N;i++) im[i]=0.0;

    fft_fwd();
    bitrev();

    // pointwise packing in natural order: P[k] = A[k]*B[k]
    {
        double xr=re[0], xi=im[0];
        re[0]=xr*xi; im[0]=0.0;
        double hr=re[N/2], hi=im[N/2];
        re[N/2]=hr*hi; im[N/2]=0.0;
        for(int k=1;k<N/2;k++){
            double Xr=re[k], Xi=im[k];
            double Zr=re[N-k], Zi=im[N-k];
            double Sx=Xr+Zr, Sy=Xi-Zi;
            double Dx=Xr-Zr, Dy=Xi+Zi;
            double Pre=0.25*(Sx*Dy+Sy*Dx);
            double Pim=0.25*(Sy*Dy-Sx*Dx);
            re[k]=Pre; im[k]=Pim;
            re[N-k]=Pre; im[N-k]=-Pim;
        }
    }

    bitrev();
    fft_inv();

    int outlen=na+nb-1;
    u64 carry=0;
    int nout=0;
    for(int i=0;i<outlen;i++){
        u64 c=(u64)(int64_t)(re[i]+0.5)+carry;
        outlimbs[nout++]=(u32)(c%10000u);
        carry=c/10000u;
    }
    while(carry){ outlimbs[nout++]=(u32)(carry%10000u); carry/=10000u; }
    int hi=nout-1; while(hi>0 && outlimbs[hi]==0) hi--;

    char* out=di->stdout_ptr; char* o=out;
    {
        u32 v=outlimbs[hi]; char tmp[5]; int t=0;
        do{ tmp[t++]='0'+v%10; v/=10; }while(v);
        while(t>0) *o++=tmp[--t];
    }
    for(int i=hi-1;i>=0;i--){
        const char* p=tab4[outlimbs[i]];
        *o++=p[0]; *o++=p[1]; *o++=p[2]; *o++=p[3];
    }
    *o++='\n';
    di->stdout_size=(u64)(o-out);

#ifdef LOCAL_TEST
    return 0;
#else
    asm volatile("mov $60, %%eax; xor %%edi, %%edi; syscall" ::: "rax","rdi","rcx","r11","memory");
    __builtin_unreachable();
#endif
}

评测结果