代码

// 1004: 1e6-digit multiply - double FFT, six-step blocked (cache-resident), real-packing.
// split-layout (re/im separate, no shuffles), N=2^19 = 512 x 1024, base 1e4.
// 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 N1 512
#define N2 1024

static double re[N] __attribute__((aligned(32)));
static double im[N] __attribute__((aligned(32)));
static double re2[N] __attribute__((aligned(32)));
static double im2[N] __attribute__((aligned(32)));

// staged twiddles for inner radix-2 FFT (size N1, N2)
static double twr1[N1], twi1[N1];
static double twr2[N2], twi2[N2];
static double twr1i[N1], twi1i[N1];
static double twr2i[N2], twi2i[N2];
// twiddle-multiply table omega^{q*m}, q in [0,N1), m in [0,N2): index q*N2+m
static double twm_re[N], twm_im[N];

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 twiddles for radix-2 FFT of size m (sign: -1 forward, +1 inverse)
static void gen_tw_staged(double* tr, double* ti, int m, double sign){
    int off=0;
    for(int len=m; len>1; len>>=1){
        int half=len>>1, step=m/len;
        double d = 2.0*M_PI*(double)step/(double)m;
        double cd=cos(d), sd=sign*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=sign*sin(a); }
            tr[off+j]=c; ti[off+j]=s;
            double nc=c*cd-s*sd; double ns=s*cd+c*sd; c=nc; s=ns;
        }
        off+=half;
    }
}

// twiddle-multiply table omega^{q*m}, sign=-1 forward
static void gen_twm(void){
    for(int q=0;q<N1;q++){
        double d = 2.0*M_PI*(double)q/(double)N;
        double cd=cos(d), sd=-sin(d);
        double c=1.0, s=0.0;
        for(int m=0;m<N2;m++){
            if((m&63)==0){ double a=d*(double)m; c=cos(a); s=-sin(a); }
            twm_re[q*N2+m]=c; twm_im[q*N2+m]=s;
            double nc=c*cd-s*sd; double ns=s*cd+c*sd; c=nc; s=ns;
        }
    }
}

// inner radix-2 DIF + bitrev, on contiguous block re[off..off+m)
static void fft_inner(double* R, double* I, int off, int m, const double* tr, const double* ti){
    int toff=0;
    for(int len=m; len>1; len>>=1){
        int half=len>>1;
        for(int i=0;i<m;i+=len){
            int j=0;
            for(; j+3<half; j+=4){
                __m256d ur=_mm256_load_pd(&R[off+i+j]);
                __m256d ui=_mm256_load_pd(&I[off+i+j]);
                __m256d vr=_mm256_load_pd(&R[off+i+j+half]);
                __m256d vi=_mm256_load_pd(&I[off+i+j+half]);
                __m256d wr=_mm256_load_pd(&tr[toff+j]);
                __m256d wi=_mm256_load_pd(&ti[toff+j]);
                __m256d sr=_mm256_add_pd(ur,vr);
                __m256d si=_mm256_add_pd(ui,vi);
                __m256d trr=_mm256_sub_pd(ur,vr);
                __m256d ti2=_mm256_sub_pd(ui,vi);
                __m256d dr=_mm256_fmsub_pd(trr,wr,_mm256_mul_pd(ti2,wi));
                __m256d di=_mm256_fmadd_pd(trr,wi,_mm256_mul_pd(ti2,wr));
                _mm256_store_pd(&R[off+i+j],sr);
                _mm256_store_pd(&I[off+i+j],si);
                _mm256_store_pd(&R[off+i+j+half],dr);
                _mm256_store_pd(&I[off+i+j+half],di);
            }
            for(; j<half; j++){
                double ur=R[off+i+j], ui=I[off+i+j], vr=R[off+i+j+half], vi=I[off+i+j+half];
                double wr=tr[toff+j], wi=ti[toff+j];
                R[off+i+j]=ur+vr; I[off+i+j]=ui+vi;
                R[off+i+j+half]=(ur-vr)*wr-(ui-vi)*wi;
                I[off+i+j+half]=(ur-vr)*wi+(ui-vi)*wr;
            }
        }
        toff+=half;
    }
    int j=0;
    for(int i=1;i<m-1;i++){
        int k=m>>1;
        while(k<=j){ j-=k; k>>=1; }
        j+=k;
        if(i<j){ double t=R[off+i]; R[off+i]=R[off+j]; R[off+j]=t; t=I[off+i]; I[off+i]=I[off+j]; I[off+j]=t; }
    }
}

// blocked transpose: src rows x cols -> dst cols x rows (re/im)
static void transpose(double* dr, double* di, const double* sr, const double* si, int rows, int cols){
    const int B=64;
    for(int r0=0;r0<rows;r0+=B){
        for(int c0=0;c0<cols;c0+=B){
            int rmax = (r0+B<rows)?r0+B:rows;
            int cmax = (c0+B<cols)?c0+B:cols;
            for(int r=r0;r<rmax;r++){
                const double* sp = sr + (size_t)r*cols + c0;
                const double* sp2 = si + (size_t)r*cols + c0;
                double* dp = dr + (size_t)c0*rows + r;
                double* dp2 = di + (size_t)c0*rows + r;
                for(int c=c0;c<cmax;c++){
                    dp[(size_t)(c-c0)*rows] = sp[c-c0];
                    dp2[(size_t)(c-c0)*rows] = sp2[c-c0];
                }
            }
        }
    }
}

// transpose + twiddle multiply: src is n1 x n2 (A[q][m]), dst is n2 x n1 (dst[m][q] = src[q][m]*omega^{qm})
static void transpose_twiddle(double* dr, double* di, const double* sr, const double* si){
    const int B=64;
    for(int q0=0;q0<N1;q0+=B){
        for(int m0=0;m0<N2;m0+=B){
            int qmax=(q0+B<N1)?q0+B:N1;
            int mmax=(m0+B<N2)?m0+B:N2;
            for(int q=q0;q<qmax;q++){
                for(int m=m0;m<mmax;m++){
                    double wr=twm_re[q*N2+m], wi=twm_im[q*N2+m];
                    double sr2=sr[(size_t)q*N2+m], si2=si[(size_t)q*N2+m];
                    dr[(size_t)m*N1+q] = sr2*wr - si2*wi;
                    di[(size_t)m*N1+q] = sr2*wi + si2*wr;
                }
            }
        }
    }
}

// inverse transpose+twiddle: src is n2 x n1 (A[m][q]), dst is n1 x n2 (dst[q][m] = src[m][q]*conj(omega^{qm}))
static void transpose_twiddle_inv(double* dr, double* di, const double* sr, const double* si){
    const int B=64;
    for(int q0=0;q0<N1;q0+=B){
        for(int m0=0;m0<N2;m0+=B){
            int qmax=(q0+B<N1)?q0+B:N1;
            int mmax=(m0+B<N2)?m0+B:N2;
            for(int q=q0;q<qmax;q++){
                for(int m=m0;m<mmax;m++){
                    double wr=twm_re[q*N2+m], wi=-twm_im[q*N2+m];
                    double sr2=sr[(size_t)m*N1+q], si2=si[(size_t)m*N1+q];
                    dr[(size_t)q*N2+m] = sr2*wr - si2*wi;
                    di[(size_t)q*N2+m] = sr2*wi + si2*wr;
                }
            }
        }
    }
}

// six-step forward (sign<0) / inverse (sign>0)
static void fft6(double sign){
    if(sign<0){
        // 1. transpose re (N2 x N1) -> re2 (N1 x N2)
        transpose(re2, im2, re, im, N2, N1);
        // 2. N2-point FFT on each of N1 rows of re2
        for(int q=0;q<N1;q++) fft_inner(re2, im2, q*N2, N2, twr2, twi2);
        // 3. transpose + twiddle re2 (N1xN2) -> re (N2xN1)
        transpose_twiddle(re, im, re2, im2);
        // 4. N1-point FFT on each of N2 rows of re
        for(int m=0;m<N2;m++) fft_inner(re, im, m*N1, N1, twr1, twi1);
        // 5. transpose re (N2 x N1) -> re2 (N1 x N2) = natural order
        transpose(re2, im2, re, im, N2, N1);
    } else {
        // inverse
        // 1. undo fwd step5: re2 (N1xN2) -> re (N2xN1)
        transpose(re, im, re2, im2, N1, N2);
        // 2. inverse N1-point FFT on each of N2 rows of re
        for(int m=0;m<N2;m++) fft_inner(re, im, m*N1, N1, twr1i, twi1i);
        // 3. inverse transpose+twiddle re (N2xN1) -> re2 (N1xN2)
        transpose_twiddle_inv(re2, im2, re, im);
        // 4. inverse N2-point FFT on each of N1 rows of re2
        for(int q=0;q<N1;q++) fft_inner(re2, im2, q*N2, N2, twr2i, twi2i);
        // 5. undo fwd step1: re2 (N1xN2) -> re (N2xN1)
        transpose(re, im, re2, im2, N1, N2);
        __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));
        }
    }
}

#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_staged(twr1, twi1, N1, -1.0);
    gen_tw_staged(twr2, twi2, N2, -1.0);
    gen_tw_staged(twr1i, twi1i, N1, +1.0);
    gen_tw_staged(twr2i, twi2i, N2, +1.0);
    // inverse inner twiddles: conjugate of forward (reuse with negated sin)
    // note: fft_inner for inverse uses conjugate twiddles; we store both.
    // (we pass -1 twiddles to forward, and for inverse we need +sin; handle via separate tables below)
    // We'll regenerate: for simplicity, keep twr1/twi1/twr2/twi2 as FORWARD, and use negated sin for inverse.
    gen_twm();

    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;

    fft6(-1.0);  // forward, result in re2/im2 natural order

    // pointwise packing in natural order (operate on re2/im2)
    {
        double xr=re2[0], xi=im2[0];
        re2[0]=xr*xi; im2[0]=0.0;
        double hr=re2[N/2], hi=im2[N/2];
        re2[N/2]=hr*hi; im2[N/2]=0.0;
        for(int k=1;k<N/2;k++){
            double Xr=re2[k], Xi=im2[k];
            double Zr=re2[N-k], Zi=im2[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);
            re2[k]=Pre; im2[k]=Pim;
            re2[N-k]=Pre; im2[N-k]=-Pim;
        }
    }

    fft6(+1.0);  // inverse, result in re/im natural order

    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
}

评测结果