#ifndef FUSE16
#define FUSE16 1
#endif
#ifndef NTPARSE
#define NTPARSE 0
#endif
#ifndef NOMEMSET
#define NOMEMSET 1
#endif
#ifndef PFV
#define PFV 4
#endif
#ifndef PFW
#define PFW 2
#endif
#ifndef PFV
#define PFV 0
#endif
#ifndef PFW
#define PFW 0
#endif
#ifndef ROWPF
#define ROWPF 0
#endif
#ifndef NTPARSE
#define NTPARSE 0
#endif
#ifndef NTCARRY
#define NTCARRY 0
#endif
#ifndef FUSE16
#define FUSE16 0
#endif
#ifndef NOMEMSET
#define NOMEMSET 0
#endif
// CBW=32 variant of v26 (tiled column pass)
#pragma GCC optimize("O3","unroll-loops","rename-registers")
// 1004 -- multiply two 1,000,000-digit decimals. base 10^4 limbs, packed z=a+ib.
#include <cstdint>
#include <cstring>
#include <cmath>
#include <immintrin.h>
#ifndef CBW
#define CBW 16
#endif
#ifndef ZSKIP
#define ZSKIP 1
#endif
#ifndef CBW2
#define CBW2 8
#endif
#ifndef PF
#define PF 1
#endif
#define LOGN 19
#define NN (1<<LOGN)
#define NR (1<<10)
#define NC (1<<9)
#define STRIDE (NC+8)
#define ARRSZ ((size_t)NR*STRIDE)
typedef uint32_t u32; typedef uint64_t u64; typedef int64_t i64;
typedef double f64;
struct DI {
unsigned long abi;
const char *s; unsigned long sn;
char *o; unsigned long ol; unsigned long os;
char *e; unsigned long el; unsigned long es;
const char *IB; unsigned long IBl;
char *OB; unsigned long OBl;
unsigned long tsc;
} __attribute__((packed));
#define TGT __attribute__((target("avx2,fma")))
#define FT __attribute__((target("avx2,fma")))
// RE and IM must NOT be congruent mod 4096: the butterfly touches 4 rows x 8 lines of
// RE and the same of IM, and with a page-congruent IM those 64 lines fall into the same
// 8 L1 sets (8 ways -- the entire L1). Measured 29.4 -> 19.4 cyc/butterfly when the two
// arrays are 1 KB apart mod 4 KB instead of 0.
alignas(4096) static f64 RE_[2*ARRSZ + 256];
static f64 *const RE = RE_;
static f64 *const IM = RE_ + ARRSZ + 128; // +1024 B mod 4096
/* reformat18 twiddles: CJ0[j]=j/2 and CJ1[j]=256+j/2 for even j, so the destination is
plain sequential and only the DQ lookup is permuted. Precompute per kq so the vector
kernel loads four in a row. */
static f64 TWR[2][NC], TWI[2][NC];
alignas(64) static f64 RDr[NC], RDi[NC], RIr[NC], RIi[NC];
static f64 CDr[NR], CDi[NR], CIr[NR], CIi[NR];
static f64 DQr[NR], DQi[NR], DTr[NC], DTi[NC];
static u32 BREV[NR], BREVC[NC], BREV9[512], CJ0[NC], CJ1[NC];
static f64 DQ18[512], DQ18i[512], DT18[NC], DT18i[NC];
#define NN2 ((1<<18))
static u32 RVIX[NR], CVIX0[NC], CVIX1[NC];
static f64 LIN[(size_t)NN + 8];
static char OUTB[(size_t)2*1000000 + 2*1000000/8 + 4096];
static u32 D8[(1<<19)];
static char DIG2T[201] =
"00010203040506070809101112131415161718192021222324252627282930313233343536373839"
"40414243444546474849505152535455565758596061626364656667686970717273747576777879"
"8081828384858687888990919293949596979899";
FT static void build_r4(void);
FT static void build_all(void) {
/* base transcendentals: 2 048 sin/cos pairs (was 6 390) */
for (int q=0;q<NR;q++){ f64 a=-2.0*M_PI*(f64)(q*NC)/(f64)NN; DQr[q]=cos(a); DQi[q]=sin(a); }
for (int t=0;t<NC;t++){ f64 a=-2.0*M_PI*(f64)t/(f64)NN; DTr[t]=cos(a); DTi[t]=sin(a); }
for (int t=0;t<NC;t++){ f64 a=-2.0*M_PI*(f64)t/(f64)NN2; DT18[t]=cos(a); DT18i[t]=sin(a); }
build_r4(); /* T/RT/RDT/CDT all from DQr/DQi */
for (int i=0;i<NR;i++){ u32 r=0; for(int b=0;b<10;b++) if(i&(1u<<b)) r|=1u<<(9-b); BREV[i]=r; }
for (int i=0;i<512;i++){ u32 r=0; for(int b=0;b<9;b++) if(i&(1u<<b)) r|=1u<<(8-b); BREV9[i]=r; }
for (int q=0;q<512;q++){ DQ18[q]=DQr[2*q]; DQ18i[q]=DQi[2*q]; }
for (int j=0;j<NC;j++){ u32 r=0; for(int b=0;b<9;b++) if(j&(1u<<b)) r|=1u<<(8-b); BREVC[j]=r; }
for (int i=0;i<NR;i++){ u32 k1=BREV[i]; u32 k1n=(NR-k1)%NR; RVIX[i]=BREV[k1n]; }
for (int j=0;j<NC;j++){
u32 k2=BREVC[j];
CVIX0[j]=BREVC[(NC-k2)%NC];
CVIX1[j]=BREVC[(NC-k2+NC-1)%NC];
}
for (int j=0;j<NC;j+=2){
u32 k2 = BREVC[j];
for (int par=0; par<2; par++) {
u32 k2p = par + 2*k2;
u32 cc = BREVC[k2p & 511];
if (par) CJ1[j]=cc; else CJ0[j]=cc;
}
}
{ int n=NC, k=0;
for (int h=n>>1;h>=1;h>>=1){ for(int j=0;j<h;j++){ int q=512*j/h; RDr[k+j]=DQr[q]; RDi[k+j]=DQi[q];} k+=h; }
k=0;
for(int h=1;h<n;h<<=1){ for(int j=0;j<h;j++){ RIr[k+j]=RDr[k+j]; RIi[k+j]=-RDi[k+j]; } k+=h; } }
for (int kq=0;kq<2;kq++)
for (int c=0;c<NC;c++){ u32 eh = kq + 2u*BREVC[2*c]; TWR[kq][c]=DQr[eh]; TWI[kq][c]=DQi[eh]; }
{ int n=NR, k=0;
for (int h=n>>1;h>=1;h>>=1){ for(int j=0;j<h;j++){ int q=512*j/h; CDr[k+j]=DQr[q]; CDi[k+j]=DQi[q];} k+=h; }
k=0;
for(int h=1;h<n;h<<=1){ for(int j=0;j<h;j++){ CIr[k+j]=CDr[k+j]; CIi[k+j]=-CDi[k+j]; } k+=h; } }
}
FT static inline void low3_dif(f64*pr, f64*pi) {
__m256d r0=_mm256_load_pd(pr), r1=_mm256_load_pd(pr+4);
__m256d i0=_mm256_load_pd(pi), i1=_mm256_load_pd(pi+4);
{ const f64 twr[4]={1.0, 0.7071067811865475244, 0.0, -0.7071067811865475244};
const f64 twi[4]={0.0,-0.7071067811865475244,-1.0, -0.7071067811865475244};
__m256d wr=_mm256_load_pd(twr), wi=_mm256_load_pd(twi);
__m256d sr=_mm256_add_pd(r0,r1), si=_mm256_add_pd(i0,i1);
__m256d dr=_mm256_sub_pd(r0,r1), di=_mm256_sub_pd(i0,i1);
__m256d t2=_mm256_mul_pd(di,wi), drr=_mm256_fmsub_pd(dr,wr,t2);
__m256d t3=_mm256_mul_pd(dr,wi), dii=_mm256_fmadd_pd(di,wr,t3);
r0=sr;i0=si;r1=drr;i1=dii; }
#define H2D(V,W) do{ \
__m256d tv=_mm256_permute2f128_pd(V,V,1), tw=_mm256_permute2f128_pd(W,W,1); \
__m256d sV=_mm256_add_pd(V,tv), dV=_mm256_sub_pd(V,tv); \
__m256d sW=_mm256_add_pd(W,tw), dW=_mm256_sub_pd(W,tw); \
__m256d M=_mm256_permute2f128_pd(dV,dW,0x20); \
__m256d P1=_mm256_permute4x64_pd(M,0xC0); \
__m256d P2=_mm256_permute4x64_pd(M,_MM_SHUFFLE(1,2,0,0)); \
__m256d sg=_mm256_set_pd(-1.0,1.0,1.0,1.0); \
V=_mm256_blend_pd(sV,P1,0xC); \
W=_mm256_blend_pd(sW,_mm256_mul_pd(P2,sg),0xC); \
}while(0)
H2D(r0,i0); H2D(r1,i1);
#undef H2D
#define H1D(V,W) do{ \
__m256d tv=_mm256_permute_pd(V,5), tw=_mm256_permute_pd(W,5); \
V=_mm256_blend_pd(_mm256_add_pd(V,tv),_mm256_permute4x64_pd(_mm256_sub_pd(V,tv),0x80),0xA); \
W=_mm256_blend_pd(_mm256_add_pd(W,tw),_mm256_permute4x64_pd(_mm256_sub_pd(W,tw),0x80),0xA); \
}while(0)
H1D(r0,i0); H1D(r1,i1);
#undef H1D
_mm256_store_pd(pr,r0); _mm256_store_pd(pi,i0);
_mm256_store_pd(pr+4,r1); _mm256_store_pd(pi+4,i1);
}
FT static inline void low3_dit(f64*pr, f64*pi) {
__m256d r0=_mm256_load_pd(pr), r1=_mm256_load_pd(pr+4);
__m256d i0=_mm256_load_pd(pi), i1=_mm256_load_pd(pi+4);
#define H1I(V,W) do{ \
__m256d tv=_mm256_permute_pd(V,5), tw=_mm256_permute_pd(W,5); \
V=_mm256_blend_pd(_mm256_add_pd(V,tv),_mm256_permute4x64_pd(_mm256_sub_pd(V,tv),0x80),0xA); \
W=_mm256_blend_pd(_mm256_add_pd(W,tw),_mm256_permute4x64_pd(_mm256_sub_pd(W,tw),0x80),0xA); \
}while(0)
H1I(r0,i0); H1I(r1,i1);
#undef H1I
#define H2I(V,W) do{ \
__m256d tv=_mm256_permute2f128_pd(V,V,1), tw=_mm256_permute2f128_pd(W,W,1); \
__m256d Ur=_mm256_blend_pd(tv,_mm256_sub_pd(_mm256_setzero_pd(),tw),0x2); \
__m256d Ui=_mm256_blend_pd(tw,tv,0x2); \
__m256d sV=_mm256_add_pd(V,Ur), dV=_mm256_sub_pd(V,Ur); \
__m256d sW=_mm256_add_pd(W,Ui), dW=_mm256_sub_pd(W,Ui); \
V=_mm256_blend_pd(sV,_mm256_permute4x64_pd(dV,0x40),0xC); \
W=_mm256_blend_pd(sW,_mm256_permute4x64_pd(dW,0x40),0xC); \
}while(0)
H2I(r0,i0); H2I(r1,i1);
#undef H2I
{ const f64 twr[4]={1.0, 0.7071067811865475244, 0.0, -0.7071067811865475244};
const f64 twi[4]={0.0, 0.7071067811865475244, 1.0, 0.7071067811865475244};
__m256d wr=_mm256_load_pd(twr), wi=_mm256_load_pd(twi);
__m256d t2=_mm256_mul_pd(i1,wi), prr=_mm256_fmsub_pd(r1,wr,t2);
__m256d t3=_mm256_mul_pd(r1,wi), pii=_mm256_fmadd_pd(i1,wr,t3);
__m256d s=_mm256_add_pd(r0,prr), d=_mm256_sub_pd(r0,prr);
__m256d s2=_mm256_add_pd(i0,pii), d2=_mm256_sub_pd(i0,pii);
r0=s;i0=s2;r1=d;i1=d2; }
_mm256_store_pd(pr,r0); _mm256_store_pd(pi,i0);
_mm256_store_pd(pr+4,r1); _mm256_store_pd(pi+4,i1);
}
// Radix-4 DIF butterfly with outputs written as soon as each is complete. The
// monolithic form (compute all of y0,y1,y2,y3 then store) needs ~24 live YMM values
// against 16 architectural registers, so gcc spilled and the loop ran at IPC ~0.9.
#define R4BODY(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
{ __m256d y2r=_mm256_sub_pd(ar,cr), y2i=_mm256_sub_pd(ai,ci); \
__m256d u=_mm256_mul_pd(y2i,W2I), z2r=_mm256_fmsub_pd(y2r,W2R,u); \
__m256d v=_mm256_mul_pd(y2r,W2I), z2i=_mm256_fmadd_pd(y2i,W2R,v); \
_mm256_store_pd(R1+c,z2r); _mm256_store_pd(I1+c,z2i); } \
{ __m256d u1r=_mm256_add_pd(br,di), u1i=_mm256_sub_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u1i,W1I), y1r=_mm256_fmsub_pd(u1r,W1R,u); \
__m256d v=_mm256_mul_pd(u1r,W1I), y1i=_mm256_fmadd_pd(u1i,W1R,v); \
_mm256_store_pd(R2+c,y1r); _mm256_store_pd(I2+c,y1i); } \
{ __m256d u3r=_mm256_sub_pd(br,di), u3i=_mm256_add_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u3i,W3I), y3r=_mm256_fmsub_pd(u3r,W3R,u); \
__m256d v=_mm256_mul_pd(u3r,W3I), y3i=_mm256_fmadd_pd(u3i,W3R,v); \
_mm256_store_pd(R3+c,y3r); _mm256_store_pd(I3+c,y3i); } }while(0)
#define R4BODY0(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_setzero_pd(),x2i=x2r,x3r=x2r,x3i=x2r; /* rows 512..1023 are exactly +0.0 */ \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
{ __m256d y2r=_mm256_sub_pd(ar,cr), y2i=_mm256_sub_pd(ai,ci); \
__m256d u=_mm256_mul_pd(y2i,W2I), z2r=_mm256_fmsub_pd(y2r,W2R,u); \
__m256d v=_mm256_mul_pd(y2r,W2I), z2i=_mm256_fmadd_pd(y2i,W2R,v); \
_mm256_store_pd(R1+c,z2r); _mm256_store_pd(I1+c,z2i); } \
{ __m256d u1r=_mm256_add_pd(br,di), u1i=_mm256_sub_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u1i,W1I), y1r=_mm256_fmsub_pd(u1r,W1R,u); \
__m256d v=_mm256_mul_pd(u1r,W1I), y1i=_mm256_fmadd_pd(u1i,W1R,v); \
_mm256_store_pd(R2+c,y1r); _mm256_store_pd(I2+c,y1i); } \
{ __m256d u3r=_mm256_sub_pd(br,di), u3i=_mm256_add_pd(bi,dr); \
__m256d u=_mm256_mul_pd(u3i,W3I), y3r=_mm256_fmsub_pd(u3r,W3R,u); \
__m256d v=_mm256_mul_pd(u3r,W3I), y3i=_mm256_fmadd_pd(u3i,W3R,v); \
_mm256_store_pd(R3+c,y3r); _mm256_store_pd(I3+c,y3i); } }while(0)
// Radix-4 DIT butterfly, fusing the radix-2 DIT stages of half-sizes h and 2h.
// Rows are (s+j, s+j+h, s+j+2h, s+j+3h); W1 = exp(i*pi*j/h) is the stage-h twiddle
// (used twice, on x1 and x3); W2 = exp(i*pi*j/(2h)) is the stage-2h twiddle for the
// (A,C) pair and W3 = i*W2 the one for the (B,D) pair, whose within-group offset is
// j+h, i.e. exp(i*pi*j/(2h))*exp(i*pi/2).
#define R4DIT(R0,R1,R2,R3,I0,I1,I2,I3,W1R,W1I,W2R,W2I,W3R,W3I) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d a1r=_mm256_fmsub_pd(x1r,W1R,_mm256_mul_pd(x1i,W1I)); \
__m256d a1i=_mm256_fmadd_pd(x1r,W1I,_mm256_mul_pd(x1i,W1R)); \
__m256d a3r=_mm256_fmsub_pd(x3r,W1R,_mm256_mul_pd(x3i,W1I)); \
__m256d a3i=_mm256_fmadd_pd(x3r,W1I,_mm256_mul_pd(x3i,W1R)); \
__m256d Ar=_mm256_add_pd(x0r,a1r), Ai=_mm256_add_pd(x0i,a1i); \
__m256d Br=_mm256_sub_pd(x0r,a1r), Bi=_mm256_sub_pd(x0i,a1i); \
__m256d Cr=_mm256_add_pd(x2r,a3r), Ci=_mm256_add_pd(x2i,a3i); \
__m256d Dr=_mm256_sub_pd(x2r,a3r), Di=_mm256_sub_pd(x2i,a3i); \
__m256d c2r=_mm256_fmsub_pd(Cr,W2R,_mm256_mul_pd(Ci,W2I)); \
__m256d c2i=_mm256_fmadd_pd(Cr,W2I,_mm256_mul_pd(Ci,W2R)); \
__m256d d2r=_mm256_fmsub_pd(Dr,W3R,_mm256_mul_pd(Di,W3I)); \
__m256d d2i=_mm256_fmadd_pd(Dr,W3I,_mm256_mul_pd(Di,W3R)); \
_mm256_store_pd(R0+c,_mm256_add_pd(Ar,c2r)); \
_mm256_store_pd(I0+c,_mm256_add_pd(Ai,c2i)); \
_mm256_store_pd(R2+c,_mm256_sub_pd(Ar,c2r)); \
_mm256_store_pd(I2+c,_mm256_sub_pd(Ai,c2i)); \
_mm256_store_pd(R1+c,_mm256_add_pd(Br,d2r)); \
_mm256_store_pd(I1+c,_mm256_add_pd(Bi,d2i)); \
_mm256_store_pd(R3+c,_mm256_sub_pd(Br,d2r)); \
_mm256_store_pd(I3+c,_mm256_sub_pd(Bi,d2i)); }while(0)
#define R4BODYT(R0,R1,R2,R3,I0,I1,I2,I3) do{ \
__m256d x0r=_mm256_load_pd(R0+c),x0i=_mm256_load_pd(I0+c); \
__m256d x1r=_mm256_load_pd(R1+c),x1i=_mm256_load_pd(I1+c); \
__m256d x2r=_mm256_load_pd(R2+c),x2i=_mm256_load_pd(I2+c); \
__m256d x3r=_mm256_load_pd(R3+c),x3i=_mm256_load_pd(I3+c); \
__m256d ar=_mm256_add_pd(x0r,x2r), ai=_mm256_add_pd(x0i,x2i); \
__m256d br=_mm256_sub_pd(x0r,x2r), bi=_mm256_sub_pd(x0i,x2i); \
__m256d cr=_mm256_add_pd(x1r,x3r), ci=_mm256_add_pd(x1i,x3i); \
__m256d dr=_mm256_sub_pd(x1r,x3r), di=_mm256_sub_pd(x1i,x3i); \
_mm256_store_pd(R0+c,_mm256_add_pd(ar,cr)); \
_mm256_store_pd(I0+c,_mm256_add_pd(ai,ci)); \
_mm256_store_pd(R1+c,_mm256_sub_pd(ar,cr)); \
_mm256_store_pd(I1+c,_mm256_sub_pd(ai,ci)); \
_mm256_store_pd(R2+c,_mm256_add_pd(br,di)); \
_mm256_store_pd(I2+c,_mm256_sub_pd(bi,dr)); \
_mm256_store_pd(R3+c,_mm256_sub_pd(br,di)); \
_mm256_store_pd(I3+c,_mm256_add_pd(bi,dr)); }while(0)
#ifndef NST
#define NST 5
#endif
alignas(64) static f64 RT1r[176],RT1i[176],RT2r[176],RT2i[176],RT3r[176],RT3i[176];
alignas(64) static f64 RDT1r[168],RDT1i[168],RDT2r[168],RDT2i[168],RDT3r[168],RDT3i[168];
alignas(64) static f64 CDT1r[85],CDT1i[85],CDT2r[85],CDT2i[85],CDT3r[85],CDT3i[85];
FT static void row_dif(f64 *pr, f64 *pi) {
const int n=NC; int k=0;
// Radix-4 DIF: fuses the radix-2 stages (256,128),(64,32),(16,8) into three passes;
// low3_dif still covers (4,2,1), so the total is the same 9 radix-2 equivalents and the
// output ordering is bit-for-bit what the radix-2 kernel produced.
for (int h=n>>2;h>=8;k+=h,h>>=2)
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RT1r+k+j), w1i=_mm256_load_pd(RT1i+k+j);
__m256d w2r=_mm256_load_pd(RT2r+k+j), w2i=_mm256_load_pd(RT2i+k+j);
__m256d w3r=_mm256_load_pd(RT3r+k+j), w3i=_mm256_load_pd(RT3i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4BODY(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int s=0;s<n;s+=16) { low3_dif(pr+s,pi+s); low3_dif(pr+s+8,pi+s+8); }
}
static f64 T1r[400],T1i[400],T2r[400],T2i[400],T3r[400],T3i[400];
// Row-pass radix-4 twiddles. Fusing the radix-2 DIF stages (256,128),(64,32),(16,8) of the
// 512-point row transform needs W_t = exp(-2*pi*i*t*j/(2*h2)) with h2 the LARGER radix-2
// half-size, i.e. h_code=h2/2 and w=-2*pi/(4*h_code) -- the same convention build_r4 uses
// for the column pass, so the tables cannot be shared.
/* Twiddle tables are DERIVED from DQr/DQi instead of evaluating sin/cos again.
A bit-equality checker over all 11 732 entries reports 0 mismatches, so every
substitution below is exact, not approximate: the donor's argument is always
the target's argument times an exact power of two, and rounding commutes with
scaling by a power of two. build_all's transcendental calls drop 6 390 -> 2 048. */
FT static void build_r4(void){
int k=0;
for (int h=NR>>2; h>=1; h>>=2){ /* w = -2pi/(4h); donor q = 256*j*m/h */
for(int j=0;j<h;j++){
int q1=256*j/h, q2=512*j/h, q3=768*j/h;
T1r[k+j]=DQr[q1]; T1i[k+j]=DQi[q1];
T2r[k+j]=DQr[q2]; T2i[k+j]=DQi[q2];
T3r[k+j]=DQr[q3]; T3i[k+j]=DQi[q3];
} k+=h; }
k=0;
for (int h=NC>>2; h>=8; h>>=2){
for(int j=0;j<h;j++){
int q1=256*j/h, q2=512*j/h, q3=768*j/h;
RT1r[k+j]=DQr[q1]; RT1i[k+j]=DQi[q1];
RT2r[k+j]=DQr[q2]; RT2i[k+j]=DQi[q2];
RT3r[k+j]=DQr[q3]; RT3i[k+j]=DQi[q3];
} k+=h; }
{ int m=0;
for (int h=8; h<NC; h<<=2){ /* aa=pi/h, t=aa*j, u=t*0.5 */
for(int j=0;j<h;j++){ int q1=512*j/h, q2=256*j/h;
RDT1r[m+j]=DQr[q1]; RDT1i[m+j]=-DQi[q1];
RDT2r[m+j]=DQr[q2]; RDT2i[m+j]=-DQi[q2];
RDT3r[m+j]=DQi[q2]; RDT3i[m+j]=DQr[q2];
} m+=h; } }
{ int m=0;
for (int h=1; h*4<=512; h<<=2){
for(int j=0;j<h;j++){ int q1=512*j/h, q2=256*j/h;
CDT1r[m+j]=DQr[q1]; CDT1i[m+j]=-DQi[q1];
CDT2r[m+j]=DQr[q2]; CDT2i[m+j]=-DQi[q2];
CDT3r[m+j]=DQi[q2]; CDT3i[m+j]=DQr[q2];
} m+=h; } }
}
FT static void col_dif(int j0){ const int n=NR;
// Four-step column pass. Twiddle table offsets: the radix-4 passes fuse the radix-2
// half-sizes (512,256),(128,64),(32,16),(8,4),(2,1); build_r4 lays them out at
// cumulative offsets 0,256,320,336,340.
// The tail passes are TILED so that each pass works on the smallest window its
// butterflies need: h=64 over a 256-row window, h=16 over 64 rows, h=4 over 16 rows
// and h=1 over 4 rows -- those last two are L1-resident, so their re-reads never
// reach L2/L3.
{
const int h=n>>2; const int K=0;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K+j]),w1i=_mm256_set1_pd(T1i[K+j]);
__m256d w2r=_mm256_set1_pd(T2r[K+j]),w2i=_mm256_set1_pd(T2i[K+j]);
__m256d w3r=_mm256_set1_pd(T3r[K+j]),w3i=_mm256_set1_pd(T3i[K+j]);
const size_t st1=(size_t)h*STRIDE, st4=(size_t)4*h*STRIDE;
f64 *r0=RE+(size_t)j*STRIDE+j0,*i0=IM+(size_t)j*STRIDE+j0;
/* FIXED: the old block advanced char* p0 by st1 DOUBLES (byte distance 8x too
small) so p1/p2/p3 named lines far from any address the pass reads, and
`o<CBW` (16) prefetched only the first of the two 64 B lines of a row segment.
Corrected: byte stride, both lines, distance 4 (the loop body is ~60 cycles,
DRAM ~200). Judge-priced: -130.4 us of a 5809.4 us warm run (-2.245 %),
4 distances tested (1/2/4/8 all positive, 4 best), and disabling the
prefetch entirely is worth only -7.7 us -- so the win is the CORRECT ADDRESS,
not the mere presence of a prefetch. Byte-exact on the full 2e6-digit input. */
#if 1
if (j+4<h){
const char *p0=(const char*)(RE+(size_t)(j+4)*STRIDE+j0);
const char *q0=(const char*)(IM+(size_t)(j+4)*STRIDE+j0);
const size_t stp=(size_t)h*STRIDE*8;
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); } }
#endif
for(int s=0;s<n;s+=4*h){
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
#if ZSKIP
for(int c=0;c<CBW;c+=4) R4BODY0(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
#else
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
#endif
r0+=st4; i0+=st4; }
}
}
const int GW=n>>2, K64=n>>2, K16=K64+(n>>4), K4=K16+(n>>6);
for (int g=0; g<n; g+=GW) {
{ const int h=GW>>2; // 64
const size_t st1=(size_t)h*STRIDE, st4=(size_t)4*h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K64+j]),w1i=_mm256_set1_pd(T1i[K64+j]);
__m256d w2r=_mm256_set1_pd(T2r[K64+j]),w2i=_mm256_set1_pd(T2i[K64+j]);
__m256d w3r=_mm256_set1_pd(T3r[K64+j]),w3i=_mm256_set1_pd(T3i[K64+j]);
f64 *r0=RE+(size_t)(g+j)*STRIDE+j0,*i0=IM+(size_t)(g+j)*STRIDE+j0;
#if 1
if (j+1<h){
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(g+j+1)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(g+j+1)*STRIDE+j0);
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
#endif
for(int s=g;s<g+GW;s+=4*h){
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
r0+=st4; i0+=st4; }
}
}
for (int g2=g; g2<g+GW; g2+=64) {
{ const int h=16;
const size_t st1=(size_t)h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K16+j]),w1i=_mm256_set1_pd(T1i[K16+j]);
__m256d w2r=_mm256_set1_pd(T2r[K16+j]),w2i=_mm256_set1_pd(T2i[K16+j]);
__m256d w3r=_mm256_set1_pd(T3r[K16+j]),w3i=_mm256_set1_pd(T3i[K16+j]);
f64 *r0=RE+(size_t)(g2+j)*STRIDE+j0,*i0=IM+(size_t)(g2+j)*STRIDE+j0;
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
#if 1
if (j+1<h){
const size_t stp=(size_t)h*STRIDE*8;
const char*p0=(const char*)(RE+(size_t)(g2+j+1)*STRIDE+j0);
const char*q0=(const char*)(IM+(size_t)(g2+j+1)*STRIDE+j0);
for (int o=0;o<CBW*8;o+=64){
_mm_prefetch(p0+o,_MM_HINT_T0); _mm_prefetch(p0+stp+o,_MM_HINT_T0);
_mm_prefetch(p0+2*stp+o,_MM_HINT_T0); _mm_prefetch(p0+3*stp+o,_MM_HINT_T0);
_mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(q0+stp+o,_MM_HINT_T0);
_mm_prefetch(q0+2*stp+o,_MM_HINT_T0); _mm_prefetch(q0+3*stp+o,_MM_HINT_T0); }
}
#endif
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int g3=g2; g3<g2+64; g3+=16) {
{ const int h=4;
const size_t st1=(size_t)h*STRIDE;
for(int j=0;j<h;j++){
__m256d w1r=_mm256_set1_pd(T1r[K4+j]),w1i=_mm256_set1_pd(T1i[K4+j]);
__m256d w2r=_mm256_set1_pd(T2r[K4+j]),w2i=_mm256_set1_pd(T2i[K4+j]);
__m256d w3r=_mm256_set1_pd(T3r[K4+j]),w3i=_mm256_set1_pd(T3i[K4+j]);
f64 *r0=RE+(size_t)(g3+j)*STRIDE+j0,*i0=IM+(size_t)(g3+j)*STRIDE+j0;
f64 *r1=r0+st1,*i1=i0+st1, *r2=r1+st1,*i2=i1+st1, *r3=r2+st1,*i3=i2+st1;
for(int c=0;c<CBW;c+=4) R4BODY(r0,r1,r2,r3,i0,i1,i2,i3,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
for (int g4=g3; g4<g3+16; g4+=4) {
f64 *r0=RE+(size_t)g4*STRIDE+j0,*i0=IM+(size_t)g4*STRIDE+j0;
f64 *r1=r0+STRIDE,*i1=i0+STRIDE, *r2=r1+STRIDE,*i2=i1+STRIDE, *r3=r2+STRIDE,*i3=i2+STRIDE;
for(int c=0;c<CBW;c+=4) R4BODYT(r0,r1,r2,r3,i0,i1,i2,i3);
}
}
}
}
}
FT static void row_dit(f64 *pr, f64 *pi) {
const int n=NC;
for (int s=0;s<n;s+=16) { low3_dit(pr+s,pi+s); low3_dit(pr+s+8,pi+s+8); }
// Fused radix-4 DIT over the remaining stages (8,16),(32,64),(128,256).
int k=0;
for (int h=8; h*4<=n; h<<=2) {
for (int j=0;j<h;j+=4) {
__m256d w1r=_mm256_load_pd(RDT1r+k+j), w1i=_mm256_load_pd(RDT1i+k+j);
__m256d w2r=_mm256_load_pd(RDT2r+k+j), w2i=_mm256_load_pd(RDT2i+k+j);
__m256d w3r=_mm256_load_pd(RDT3r+k+j), w3i=_mm256_load_pd(RDT3i+k+j);
for (int s=0;s<n;s+=4*h) {
f64 *a=pr+s+j, *b=pi+s+j;
const int c=0;
R4DIT(a,a+h,a+2*h,a+3*h, b,b+h,b+2*h,b+3*h, w1r,w1i,w2r,w2i,w3r,w3i);
}
}
k+=h;
}
}
FT static void col_dit(int j0) {
const int n=NR; int k=n-1;
for (int h=1;h<n;h<<=1) { k-=h;
for (int s=0;s<n;s+=2*h)
for (int j=0;j<h;j++) {
f64 *ar=RE+(size_t)(s+j)*STRIDE+j0, *ai=IM+(size_t)(s+j)*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
__m256d wr=_mm256_set1_pd(CIr[k+j]), wi=_mm256_set1_pd(CIi[k+j]);
for (int c=0;c<128;c+=4) {
__m256d xr=_mm256_load_pd(ar+c), xi=_mm256_load_pd(ai+c);
__m256d yr=_mm256_load_pd(br+c), yi=_mm256_load_pd(bi+c);
__m256d t2=_mm256_mul_pd(yi,wi);
__m256d rr=_mm256_fmsub_pd(yr,wr,t2);
__m256d t3=_mm256_mul_pd(yr,wi);
__m256d ii=_mm256_fmadd_pd(yi,wr,t3);
_mm256_store_pd(ar+c,_mm256_add_pd(xr,rr));
_mm256_store_pd(ai+c,_mm256_add_pd(xi,ii));
_mm256_store_pd(br+c,_mm256_sub_pd(xr,rr));
_mm256_store_pd(bi+c,_mm256_sub_pd(xi,ii));
}
}
}
}
static inline void cw(u32 e, f64 &r, f64 &i, int sign) {
f64 rr = DQr[e>>9]*DTr[e&(NC-1)] - DQi[e>>9]*DTi[e&(NC-1)];
f64 ii = DQr[e>>9]*DTi[e&(NC-1)] + DQi[e>>9]*DTr[e&(NC-1)];
if (sign<0) ii = -ii;
r=rr; i=ii;
}
FT static inline void diagonal_row(int i, int sign, f64 sc) {
{
u32 k1 = BREV[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 rr,ri; cw(k1, rr, ri, sign);
f64 r2r=rr*rr-ri*ri, r2i=2.0*rr*ri;
f64 r3r=r2r*rr-r2i*ri, r3i=r2r*ri+r2i*rr;
f64 r4r=r2r*r2r-r2i*r2i, r4i=2.0*r2r*r2i;
__m256d Cr=_mm256_set_pd(r3r, r2r, rr, 1.0);
__m256d Ci=_mm256_set_pd(r3i, r2i, ri, 0.0);
/* four independent twiddle chains: the base advances by r16=r4^4 per 16 elements,
and chains 1..3 start at base*r4, base*r4^2, base*r4^3. A single chain was
latency-bound (one dependent complex multiply per 4 elements). */
f64 r8r=r4r*r4r-r4i*r4i, r8i=2.0*r4r*r4i; /* r4^2 */
f64 r16r=r8r*r8r-r8i*r8i, r16i=2.0*r8r*r8i; /* r4^4 */
for (int j=0;j<NC;j+=64) {
u32 be = (u32)(((u64)k1*(u32)j) & (u32)(NN-1));
f64 base_r, base_i; cw(be, base_r, base_i, sign);
base_r *= sc; base_i *= sc;
f64 b1r=base_r*r4r-base_i*r4i, b1i=base_r*r4i+base_i*r4r;
f64 b2r=b1r*r4r-b1i*r4i, b2i=b1r*r4i+b1i*r4r;
f64 b3r=b2r*r4r-b2i*r4i, b3i=b2r*r4i+b2i*r4r;
#if PFDIAG>0
if (j+64<NC) { _mm_prefetch((const char*)(pr+j+64),_MM_HINT_T0); _mm_prefetch((const char*)(pi+j+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pr+j+96),_MM_HINT_T0); _mm_prefetch((const char*)(pi+j+96),_MM_HINT_T0); }
#endif
for (int t=0;t<64;t+=16) {
f64 br_,bi_; int off;
#define DCHAIN(BR,BI,OFF) do{ \
__m256d br=_mm256_set1_pd(BR), bi=_mm256_set1_pd(BI); \
__m256d wr=_mm256_fmsub_pd(br,Cr,_mm256_mul_pd(bi,Ci)); \
__m256d wi=_mm256_fmadd_pd(br,Ci,_mm256_mul_pd(bi,Cr)); \
__m256d xr=_mm256_load_pd(pr+j+t+OFF), xi=_mm256_load_pd(pi+j+t+OFF); \
__m256d q2=_mm256_mul_pd(xi,wi); \
__m256d orr=_mm256_fmsub_pd(xr,wr,q2); \
__m256d q3=_mm256_mul_pd(xr,wi); \
__m256d oii=_mm256_fmadd_pd(xi,wr,q3); \
_mm256_store_pd(pr+j+t+OFF,orr); \
_mm256_store_pd(pi+j+t+OFF,oii); }while(0)
DCHAIN(base_r,base_i,0); DCHAIN(b1r,b1i,4);
DCHAIN(b2r,b2i,8); DCHAIN(b3r,b3i,12);
#undef DCHAIN
(void)br_;(void)bi_;(void)off;
f64 n0r=base_r*r16r-base_i*r16i, n0i=base_r*r16i+base_i*r16r;
f64 n1r=b1r*r16r-b1i*r16i, n1i=b1r*r16i+b1i*r16r;
f64 n2r=b2r*r16r-b2i*r16i, n2i=b2r*r16i+b2i*r16r;
f64 n3r=b3r*r16r-b3i*r16i, n3i=b3r*r16i+b3i*r16r;
base_r=n0r;base_i=n0i; b1r=n1r;b1i=n1i; b2r=n2r;b2i=n2i; b3r=n3r;b3i=n3i;
}
}
}
}
FT static void diagonal(int sign, f64 sc) {
for (int i=0;i<NR;i++) diagonal_row(i, sign, sc);
}
FT static inline void cw18(u32 e, f64 &r, f64 &i, int sign) {
r = DQ18[e>>9]*DT18[e&(NC-1)] - DQ18i[e>>9]*DT18i[e&(NC-1)];
i = DQ18[e>>9]*DT18i[e&(NC-1)] + DQ18i[e>>9]*DT18[e&(NC-1)];
if (sign<0) i = -i;
}
FT static inline void diagonal18_row(int i, int sign, f64 sc) {
{
u32 k1 = BREV9[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 rr,ri; cw18(k1, rr, ri, sign);
f64 r2r=rr*rr-ri*ri, r2i=2.0*rr*ri;
f64 r3r=r2r*rr-r2i*ri, r3i=r2r*ri+r2i*rr;
f64 r4r=r2r*r2r-r2i*r2i, r4i=2.0*r2r*r2i;
__m256d Cr=_mm256_set_pd(r3r, r2r, rr, 1.0);
__m256d Ci=_mm256_set_pd(r3i, r2i, ri, 0.0);
f64 r8r=r4r*r4r-r4i*r4i, r8i=2.0*r4r*r4i;
f64 r16r=r8r*r8r-r8i*r8i, r16i=2.0*r8r*r8i;
for (int j=0;j<NC;j+=64) {
u32 be = (u32)(((u64)k1*(u32)j) & (u32)(NN2-1));
f64 base_r, base_i; cw18(be, base_r, base_i, sign);
base_r *= sc; base_i *= sc;
f64 b1r=base_r*r4r-base_i*r4i, b1i=base_r*r4i+base_i*r4r;
f64 b2r=b1r*r4r-b1i*r4i, b2i=b1r*r4i+b1i*r4r;
f64 b3r=b2r*r4r-b2i*r4i, b3i=b2r*r4i+b2i*r4r;
for (int t=0;t<64;t+=16) {
#define ECHAIN(BR,BI,OFF) do{ \
__m256d br=_mm256_set1_pd(BR), bi=_mm256_set1_pd(BI); \
__m256d wr=_mm256_fmsub_pd(br,Cr,_mm256_mul_pd(bi,Ci)); \
__m256d wi=_mm256_fmadd_pd(br,Ci,_mm256_mul_pd(bi,Cr)); \
__m256d xr=_mm256_load_pd(pr+j+t+OFF), xi=_mm256_load_pd(pi+j+t+OFF); \
__m256d q2=_mm256_mul_pd(xi,wi); \
__m256d orr=_mm256_fmsub_pd(xr,wr,q2); \
__m256d q3=_mm256_mul_pd(xr,wi); \
__m256d oii=_mm256_fmadd_pd(xi,wr,q3); \
_mm256_store_pd(pr+j+t+OFF,orr); \
_mm256_store_pd(pi+j+t+OFF,oii); }while(0)
ECHAIN(base_r,base_i,0); ECHAIN(b1r,b1i,4);
ECHAIN(b2r,b2i,8); ECHAIN(b3r,b3i,12);
#undef ECHAIN
f64 n0r=base_r*r16r-base_i*r16i, n0i=base_r*r16i+base_i*r16r;
f64 n1r=b1r*r16r-b1i*r16i, n1i=b1r*r16i+b1i*r16r;
f64 n2r=b2r*r16r-b2i*r16i, n2i=b2r*r16i+b2i*r16r;
f64 n3r=b3r*r16r-b3i*r16i, n3i=b3r*r16i+b3i*r16r;
base_r=n0r;base_i=n0i; b1r=n1r;b1i=n1i; b2r=n2r;b2i=n2i; b3r=n3r;b3i=n3i;
}
}
}
}
FT static void diagonal18(int sign, f64 sc) {
for (int i=0;i<512;i++) diagonal18_row(i, sign, sc);
}
// Radix-4 DIT along the strided column axis: fuse the stage pairs (1,2),(4,8),(16,32),
// (64,128), leaving only h=256 as a radix-2 pass -- 5 passes instead of 9 and ~35% fewer
// instructions. The CDT tables (built for row_dit) have exactly the (h, j) layout needed.
// Unlike row_dit, the vector here spans 4 COLUMNS (broadcast twiddles), because a column
// block is many columns wide and only the position index j carries a twiddle.
FT static void col512_dit(int j0) {
const int n=512;
int m=0;
for (int h=1; h*4<=n; h<<=2) {
for (int j=0;j<h;j++) {
__m256d w1r=_mm256_set1_pd(CDT1r[m+j]), w1i=_mm256_set1_pd(CDT1i[m+j]);
__m256d w2r=_mm256_set1_pd(CDT2r[m+j]), w2i=_mm256_set1_pd(CDT2i[m+j]);
__m256d w3r=_mm256_set1_pd(CDT3r[m+j]), w3i=_mm256_set1_pd(CDT3i[m+j]);
for (int s=0;s<n;s+=4*h) {
f64 *ar=RE+(size_t)(s+j)*STRIDE+j0, *ai=IM+(size_t)(s+j)*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
f64 *cr=br+(size_t)h*STRIDE, *ci=bi+(size_t)h*STRIDE;
f64 *dr=cr+(size_t)h*STRIDE, *di=ci+(size_t)h*STRIDE;
for (int c=0;c<CBW2;c+=4)
R4DIT(ar,br,cr,dr,ai,bi,ci,di,w1r,w1i,w2r,w2i,w3r,w3i);
}
}
m+=h;
}
{ const int h=256; // last stage, twiddles live at RIr[0..255]
for (int j=0;j<h;j++) {
__m256d wr=_mm256_set1_pd(RIr[j]), wi=_mm256_set1_pd(RIi[j]);
f64 *ar=RE+(size_t)j*STRIDE+j0, *ai=IM+(size_t)j*STRIDE+j0;
f64 *br=ar+(size_t)h*STRIDE, *bi=ai+(size_t)h*STRIDE;
for (int c=0;c<CBW2;c+=4) {
__m256d xr=_mm256_load_pd(ar+c), xi=_mm256_load_pd(ai+c);
__m256d yr=_mm256_load_pd(br+c), yi=_mm256_load_pd(bi+c);
__m256d t2=_mm256_mul_pd(yi,wi);
__m256d rr=_mm256_fmsub_pd(yr,wr,t2);
__m256d t3=_mm256_mul_pd(yr,wi);
__m256d ii=_mm256_fmadd_pd(yi,wr,t3);
_mm256_store_pd(ar+c,_mm256_add_pd(xr,rr));
_mm256_store_pd(ai+c,_mm256_add_pd(xi,ii));
_mm256_store_pd(br+c,_mm256_sub_pd(xr,rr));
_mm256_store_pd(bi+c,_mm256_sub_pd(xi,ii));
}
}
}
}
alignas(64) static f64 SVr[2][NC], SVi[2][NC];
#ifndef R18VEC
#define R18VEC 1
#endif
#ifndef PFR18
#define PFR18 1
#endif
#ifndef PFDIAG
#define PFDIAG 1
#endif
#ifndef PFPAIR
#define PFPAIR 0
#endif
/* reformat18 twiddles: CJ0[j]=j/2, CJ1[j]=256+j/2 for even j, so the destination is
plain sequential and only the DQ lookup is permuted. Precompute it per kq so the
vector kernel can load four in a row. */
#if R18VEC>0
/* Vectorised reformat18. ccol is exactly j/2 (CJ0) or 256+j/2 (CJ1), so the writes are
sequential; the only irregular part is the DQ lookup, which TWR/TWI linearise. The
FMA contraction order is the one gcc-9 emits for the scalar form (verified from asm),
so this is bit-identical to it. */
FT static inline void reformat18_row(int i) {
const f64 *pr = RE+(size_t)i*STRIDE, *pi = IM+(size_t)i*STRIDE;
int r = i>>1, par = i&1;
f64 *wr = RE+(size_t)r*STRIDE + (par?256:0), *wi = IM+(size_t)r*STRIDE + (par?256:0);
u32 k1 = BREV[i]; u32 kq = k1>>9;
f64 dr = DTr[k1 & (NC-1)], di = DTi[k1 & (NC-1)];
if (i<2) { pr = SVr[i]; pi = SVi[i]; }
const f64 *TR = TWR[kq], *TI = TWI[kq];
const __m256d sg=_mm256_set1_pd(-0.0);
const __m256d Vr=_mm256_set1_pd(dr), Vi=_mm256_set1_pd(di);
for (int c=0;c<(NC>>1);c+=4) {
#if PFR18>0
if (c+32 < (NC>>1)) { _mm_prefetch((const char*)(pr+2*c+64),_MM_HINT_T0); _mm_prefetch((const char*)(pi+2*c+64),_MM_HINT_T0);
_mm_prefetch((const char*)(pr+2*c+96),_MM_HINT_T0); _mm_prefetch((const char*)(pi+2*c+96),_MM_HINT_T0); }
#endif
__m256d A=_mm256_load_pd(pr+2*c), B=_mm256_load_pd(pr+2*c+4);
// even/odd deinterleave with 2 permutes instead of 3 (identical values):
// P=perm4x64(A,0xD8)=[a0,a2,a1,a3]; the 128-halves of (P,Q) give even and odd.
__m256d P=_mm256_permute4x64_pd(A,0xD8), Q=_mm256_permute4x64_pd(B,0xD8);
__m256d z1r=_mm256_permute2f128_pd(P,Q,0x20);
__m256d z2r=_mm256_permute2f128_pd(P,Q,0x31);
A=_mm256_load_pd(pi+2*c); B=_mm256_load_pd(pi+2*c+4);
P=_mm256_permute4x64_pd(A,0xD8); Q=_mm256_permute4x64_pd(B,0xD8);
__m256d z1i=_mm256_permute2f128_pd(P,Q,0x20);
__m256d z2i=_mm256_permute2f128_pd(P,Q,0x31);
__m256d qr0=_mm256_load_pd(TR+c), qi0=_mm256_load_pd(TI+c);
__m256d er=_mm256_add_pd(z1r,z2r), dr2=_mm256_sub_pd(z1r,z2r);
__m256d ei=_mm256_add_pd(z1i,z2i), di2=_mm256_sub_pd(z1i,z2i);
__m256d t0=_mm256_mul_pd(qr0,Vr), tr=_mm256_fnmadd_pd(qi0,Vi,t0);
__m256d t1=_mm256_mul_pd(qr0,Vi), ti=_mm256_fnmadd_pd(qi0,Vr,_mm256_xor_pd(t1,sg));
__m256d t2=_mm256_mul_pd(tr,dr2), orr=_mm256_fnmadd_pd(ti,di2,t2);
__m256d t3=_mm256_mul_pd(tr,di2), oii=_mm256_fmadd_pd(ti,dr2,t3);
_mm256_store_pd(wr+c,_mm256_sub_pd(er,oii));
_mm256_store_pd(wi+c,_mm256_add_pd(ei,orr));
}
}
#else
FT static inline void reformat18_row_scalar(int i);
FT static inline void reformat18_row(int i) { reformat18_row_scalar(i); }
#endif
FT static inline void reformat18_row_scalar(int i) {
{
const f64 *pr = RE+(size_t)i*STRIDE, *pi = IM+(size_t)i*STRIDE;
int r = i>>1, par = i&1;
const u32 *CJ = par ? CJ1 : CJ0;
f64 *wr = RE+(size_t)r*STRIDE, *wi = IM+(size_t)r*STRIDE;
u32 k1 = BREV[i];
u32 kq = k1>>9;
f64 dr = DTr[k1 & (NC-1)], di = DTi[k1 & (NC-1)];
/* rows 0,1 were clobbered by the pointwise before their own pass, so they come from
the copies taken at entry; select the base pointers once instead of testing i<2
in the inner loop (the two branches have identical layout). */
if (i<2) { pr = SVr[i]; pi = SVi[i]; }
for (int j=0;j<NC;j+=2) {
f64 z1r,z1i,z2r,z2i;
{ z1r=pr[j]; z1i=pi[j]; z2r=pr[j+1]; z2i=pi[j+1]; }
u32 eh = kq + 2u*BREVC[j];
f64 qr0=DQr[eh], qi0=DQi[eh];
f64 tr = qr0*dr - qi0*di;
f64 ti = -(qr0*di + qi0*dr);
f64 er=z1r+z2r, ei=z1i+z2i;
f64 dr=z1r-z2r, di=z1i-z2i;
f64 orr=tr*dr-ti*di, oii=tr*di+ti*dr;
int ccol = (int)CJ[j];
wr[ccol] = er - oii;
wi[ccol] = ei + orr;
}
}
}
FT static void reformat18(void) {
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
for (int i=0;i<NR;i++) reformat18_row(i);
}
FT static int parse_digits2(const char *s, long len, f64 *dst, int stride) {
// v43 -> c2: the destination pointer is now advanced incrementally (wp += 8, plus
// a row jump when the column wraps). The old form recomputed
// (base>>9)*stride + (base&511) from scratch every 32 digits, which cost ~10 of
// the loop's 22 instructions (shift, imul by stride, and, two lea) for what is
// really a constant +8 step. Identical stores, identical order.
int m = (int)((len + 3) / 4);
long lo = len & 3;
if (lo) { u32 v=0; for(long j=0;j<lo;j++) v=v*10+(u32)(s[j]-'0'); long kk=m-1;
dst[(size_t)(kk>>9)*stride + (kk & 511)] = (f64)v; }
const __m256i c10 = _mm256_setr_epi8(10,1,10,1,10,1,10,1,10,1,10,1,10,1,10,1,
10,1,10,1,10,1,10,1,10,1,10,1,10,1,10,1);
const __m256i c100 = _mm256_setr_epi16(100,1,100,1,100,1,100,1,100,1,100,1,100,1,100,1);
const __m256i adj = _mm256_set1_epi32(53328);
const __m256i rv = _mm256_setr_epi32(7,6,5,4,3,2,1,0);
long i = len;
f64 *wp = dst;
int col = 0;
// two chunks (64 digits) per iteration: halves the loop overhead and lets the two
// (load, maddubs, madd, permute, cvt, cvt, store, store) chains overlap.
while (i - 64 >= lo) {
__m256i d0 = _mm256_loadu_si256((const __m256i*)(s+i-32));
__m256i d1 = _mm256_loadu_si256((const __m256i*)(s+i-64));
i -= 64;
// the scan runs BACKWARDS, which the hardware prefetchers largely ignore; pull
// the next few lines in by hand
{ long pa = i-320; if (pa < 0) pa = 0; _mm_prefetch((const char*)(s+pa), _MM_HINT_T0);
long pb = i-448; if (pb < 0) pb = 0; _mm_prefetch((const char*)(s+pb), _MM_HINT_T0); }
__m256i v0 = _mm256_sub_epi32(_mm256_madd_epi16(_mm256_maddubs_epi16(d0,c10),c100), adj);
__m256i v1 = _mm256_sub_epi32(_mm256_madd_epi16(_mm256_maddubs_epi16(d1,c10),c100), adj);
v0 = _mm256_permutevar8x32_epi32(v0, rv);
v1 = _mm256_permutevar8x32_epi32(v1, rv);
__m256d q0=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v0));
__m256d q1=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v0,1));
__m256d q2=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v1));
__m256d q3=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v1,1));
#if NTPARSE>0
_mm256_stream_pd(wp,q0); _mm256_stream_pd(wp+4,q1);
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
_mm256_stream_pd(wp,q2); _mm256_stream_pd(wp+4,q3);
#else
_mm256_store_pd(wp,q0); _mm256_store_pd(wp+4,q1);
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
_mm256_store_pd(wp,q2); _mm256_store_pd(wp+4,q3);
#endif
wp += 8; col += 8; if (col == 512) { col = 0; wp += stride - 512; }
}
while (i - 32 >= lo) {
i -= 32;
__m256i v = _mm256_sub_epi32(_mm256_madd_epi16(
_mm256_maddubs_epi16(_mm256_loadu_si256((const __m256i*)(s+i)),c10),c100), adj);
v = _mm256_permutevar8x32_epi32(v, rv);
__m256d q0=_mm256_cvtepi32_pd(_mm256_castsi256_si128(v));
__m256d q1=_mm256_cvtepi32_pd(_mm256_extracti128_si256(v,1));
#if NTPARSE>0
_mm256_stream_pd(wp,q0); _mm256_stream_pd(wp+4,q1);
#else
_mm256_store_pd(wp,q0); _mm256_store_pd(wp+4,q1);
#endif
wp += 8; col += 8;
if (col == 512) { col = 0; wp += stride - 512; }
}
#if NTPARSE>0
_mm_sfence();
#endif
int k = (int)((i - lo) / 4);
while (i - 4 >= lo) {
i -= 4;
long kk = --k;
dst[(size_t)(kk>>9)*stride + (kk & 511)] =
(f64)((((s[i]-'0')*10 + (s[i+1]-'0'))*10 + (s[i+2]-'0'))*10 + (s[i+3]-'0'));
}
return m;
}
/* Vectorised boundary scan. The scalar original ran at ~1.09 cycles per input
byte over the whole 2 MB input: measured judge-side at 2.18 M cycles = 8.5 % of
the run, and it was absent from every earlier phase map. AVX2 here: 32 bytes
per iteration, exact same semantics (including "no digit at all" -> 0 return).
ALIGNED loads only -- gcc-9 splits every unaligned 256-bit access into
vmovdqu xmm + vinserti128 (verified in the .s), which is what the note in
JUDGE_MEM warns about. A <=31-byte scalar prologue gets us to alignment. */
TGT static inline long kb_scan(const char *s, long i, long sn, int want_digit) {
const __m256i lo = _mm256_set1_epi8('0'-1), hi = _mm256_set1_epi8('9'+1);
if (i < sn) {
while ((((uintptr_t)(s + i)) & 31) != 0) {
if (i >= sn) break;
if (((s[i] >= '0' && s[i] <= '9') ? 1 : 0) == want_digit) return i;
i++;
}
long n = sn - 32;
if (want_digit) {
for (; i <= n; i += 32) {
__m256i v = _mm256_load_si256((const __m256i *)(s + i));
unsigned m = (unsigned)_mm256_movemask_epi8(
_mm256_and_si256(_mm256_cmpgt_epi8(v, lo), _mm256_cmpgt_epi8(hi, v)));
if (m) return i + __builtin_ctz(m);
}
while (i < sn && (s[i] < '0' || s[i] > '9')) i++;
} else {
for (; i <= n; i += 32) {
__m256i v = _mm256_load_si256((const __m256i *)(s + i));
unsigned m = (unsigned)_mm256_movemask_epi8(
_mm256_and_si256(_mm256_cmpgt_epi8(v, lo), _mm256_cmpgt_epi8(hi, v)));
if (~m) return i + __builtin_ctz(~m);
}
while (i < sn && s[i] >= '0' && s[i] <= '9') i++;
}
}
return i;
}
TGT static inline int split_lines(const char *s, long sn, const char **b0, long *n0,
const char **b1, long *n1) {
long i = kb_scan(s, 0, sn, 1);
long j = kb_scan(s, i, sn, 0);
if (j == i) return 0;
*b0 = s + i; *n0 = j - i;
long p = kb_scan(s, j, sn, 1);
long k = kb_scan(s, p, sn, 0);
if (k == p) return 0;
*b1 = s + p; *n1 = k - p;
return 1;
}
#define COEF(k) RE[(size_t)((k)>>9)*STRIDE + ((k) & 511)]
static u32 TAB4[10000];
FT static void build_tab4(void) {
for (u32 d = 0; d < 10000; d++)
TAB4[d] = (u32)('0'+d/1000) | ((u32)('0'+(d/100)%10) << 8)
| ((u32)('0'+(d/10)%10) << 16) | ((u32)('0'+d%10) << 24);
}
// Fused unpack + carry/format. Base-10^8 carry chain.
// NOTE: the digit group at position m is (C_m + carry_in) mod 1e8 and the new carry is
// (C_m + carry_in) div 1e8, so ONLY ONE 64x64->128 multiply is needed per step; the
// old form computed C/1e8 and v/1e8 separately. That removes two of the three
// rax:rdx-bound mulq's from the loop-carried recurrence (carry -> total -> mulhi -> carry).
// Fused unpack + carry/format. Base-10^8 carry chain: the digit group at position m
// is (C_m + carry_in) mod 1e8 and the new carry is (C_m + carry_in) div 1e8, so only ONE
// 64x64->128 multiply is needed per step and the loop-carried recurrence is
// carry -> total -> mulhi -> carry (~5 cycles).
// Changes vs the first version: (1) the per-iteration `2*m+1 < tot` test is hoisted by
// counting the iterations that take the IM term up front, (2) the row-wrap test
// `++cc==512` is replaced by an index loop over one row, (3) the two 32-bit TAB4
// stores become a single 64-bit store.
FT static long carry_format_m(int tot, char *o, char *end) {
// ---- v43 -> c1 changes (all provably exact integer transformations) ----
// (1) double->u64 via cvtsd2si (round-to-nearest, default MXCSR). The old
// (u64)(x+0.5) is floor(x+0.5); the two agree whenever |x-n| < 0.5 for the
// true integer n, which is exactly the condition under which the old form
// was correct, so this is bit-identical. It also drops the addsd and the
// branchy >=2^63 fixup gcc emits for a double->u64 conversion.
// (2) hi/lo are taken from t4 = total/1e4 rather than from dg = total mod 1e8.
// total/1e4 does not depend on ex, so BOTH magic divisions now hang off the
// loop-carried chain instead of one behind the other. Digit group is
// unchanged: (total mod 1e8) = ((total/1e4) mod 1e4)*1e4 + total mod 1e4.
// (3) two 32-bit ASCII stores instead of salq/orq + one 64-bit store.
// (4) the tail's 429497 (ceil) magic had a fixup that only handled the
// under-estimate; the over-estimate case (dg % 1e4 >= 9938) indexed TAB4
// with a wrapped u32 -> out-of-bounds read. The tail is now plain / and %.
const int nm = (tot + 1) >> 1;
int nfull = ((tot - 1) + 1) >> 1;
if (nfull > nm) nfull = nm;
if (nfull < 0) nfull = 0;
u64 carry = 0;
char *p = end;
const f64 *pr = RE, *pi = IM;
int m = 0;
while (m < nfull) {
int cnt = nfull - m; if (cnt > NC) cnt = NC;
const f64 *qr = pr, *qi = pi;
for (int k = 0; k < cnt; k++) {
u64 C = (u64)_mm_cvtsd_si64(_mm_load_sd((const f64*)(qr + k)))
+ (u64)_mm_cvtsd_si64(_mm_load_sd((const f64*)(qi + k))) * 10000ull;
u64 total = C + carry;
u64 t4 = total / 10000ull; // independent of the carry chain
u64 ex = total / 100000000ull; // the chain's next carry
carry = ex;
u32 lo = (u32)(total - t4 * 10000ull);
u32 hi = (u32)t4 - (u32)ex * 10000u;
p -= 8;
*(u32*)p = TAB4[hi];
*(u32*)(p + 4) = TAB4[lo];
}
pr += STRIDE; pi += STRIDE; m += cnt;
}
pr = RE + (size_t)(nfull/NC)*STRIDE + (nfull%NC);
pi = IM + (size_t)(nfull/NC)*STRIDE + (nfull%NC);
for (; m < nm; m++) { // tail: no IM term
u64 C = (u64)_mm_cvtsd_si64(_mm_load_sd(pr));
u64 total = C + carry;
u64 t4 = total / 10000ull;
u64 ex = total / 100000000ull;
carry = ex;
u32 lo = (u32)(total - t4 * 10000ull);
u32 hi = (u32)t4 - (u32)ex * 10000u;
p -= 8;
*(u32*)p = TAB4[hi];
*(u32*)(p + 4) = TAB4[lo];
pr++; pi++;
}
while (carry) { u32 d = (u32)(carry % 100000000ull); carry /= 100000000ull;
u32 hi = d / 10000u, lo = d % 10000u;
p -= 8; *(u32*)p = TAB4[hi]; *(u32*)(p+4) = TAB4[lo]; }
while (p < end && *p == '0') p++;
if (p == end) *(--p) = '0';
(void)o;
return (long)(end - p);
}
FT static inline void pointwise_row(int i);
static unsigned char FWD_DONE[NR];
/* Fused forward: the diagonal twiddle and the pointwise product are both per-row
operations, so they run while their rows are still L1/L2-resident straight out of
row_dif instead of as two extra full sweeps of the 8.5 MB array. Rows are visited in
pairs (i, RVIX[i]) -- RVIX is an involution -- so both operands of every pointwise
are transformed before it runs, and the pair is always entered at its smaller index,
which keeps the pointwise application order (and hence the rounding) identical to the
old separate loop. */
#if FUSE16>0
FT static void fft_forward2(int sign, f64 sc);
#endif
FT static void fft_forward(void) {
#ifdef REVCD
for (int j0=NC-CBW;j0>=0;j0-=CBW) col_dif(j0);
#else
for (int j0=0;j0<NC;j0+=CBW) col_dif(j0);
#endif
memset(FWD_DONE, 0, sizeof FWD_DONE);
for (int i=0;i<NR;i++) {
if (FWD_DONE[i]) continue;
int m = (int)RVIX[i];
#if ROWPF>0
{ int n1=i+1; if(n1<NR && !FWD_DONE[n1]){ const char*q0=(const char*)(RE+(size_t)n1*STRIDE); const char*w0=(const char*)(IM+(size_t)n1*STRIDE);
for(int o=0;o<NC*8;o+=64){ _mm_prefetch(q0+o,_MM_HINT_T0); _mm_prefetch(w0+o,_MM_HINT_T0); } }
int n2=(int)RVIX[i<NR-1?i+1:i]; if(n2!=i+1 && !FWD_DONE[n2]){ const char*q1=(const char*)(RE+(size_t)n2*STRIDE); const char*w1=(const char*)(IM+(size_t)n2*STRIDE);
for(int o=0;o<NC*8;o+=64){ _mm_prefetch(q1+o,_MM_HINT_T0); _mm_prefetch(w1+o,_MM_HINT_T0); } } }
#endif
diagonal_row(i, +1, 1.0);
row_dif(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
if (m != i) {
diagonal_row(m, +1, 1.0);
row_dif(RE+(size_t)m*STRIDE, IM+(size_t)m*STRIDE);
FWD_DONE[m] = 1;
}
FWD_DONE[i] = 1;
pointwise_row(i);
if (m != i) pointwise_row(m);
}
{ f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
}
FT static void fft_inverse(void) {
for (int i=0;i<NR;i++) row_dit(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
diagonal(-1, 1.0/(f64)NN);
for (int j0=0;j0<NC;j0+=128) col_dit(j0);
}
/* Fused inverse. reformat18 finalises destination row r at source index i=2r+1 and
nothing later in the sweep touches row r again, so row_dit and the diagonal18 twiddle
run on row r immediately, while it is still cache-resident, instead of as two extra
sweeps of the 4 MB half-array. */
#if FUSE16>0
/* Fused: reformat18 for source rows 2r,2r+1 finalises destination row r and nothing
later in the sweep touches row r again, so row_dit + diagonal18 run on it immediately
while it is still L1/L2-resident -- and because the fold writes row r only after both
of its source rows were pointwise-complete (and every later pair reads only rows >= i),
this is exactly as safe as the separate sweep. */
FT static void fused_row_loop(int sign, f64 sc) {
for (int i=0;i<NR;i++) {
if (!FWD_DONE[i]) {
int m = (int)RVIX[i];
diagonal_row(i, +1, 1.0);
row_dif(RE+(size_t)i*STRIDE, IM+(size_t)i*STRIDE);
if (m != i) { diagonal_row(m, +1, 1.0); row_dif(RE+(size_t)m*STRIDE, IM+(size_t)m*STRIDE); }
FWD_DONE[i] = 1; FWD_DONE[m] = 1;
pointwise_row(i);
if (m != i) pointwise_row(m);
}
if (i & 1) {
int r = i>>1;
if (i == 1) { f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
reformat18_row(i-1);
reformat18_row(i);
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, sign, sc);
}
}
}
#endif
FT static void fft_inverse18(int sign, f64 sc) {
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
for (int i=0;i<NR;i++) {
reformat18_row(i);
if (i & 1) { int r = i>>1;
row_dit(RE+(size_t)r*STRIDE, IM+(size_t)r*STRIDE);
diagonal18_row(r, sign, sc); }
}
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
}
#ifndef PWVEC
#define PWVEC 1
#endif
/* Vectorised pointwise. With U = row i at the 8 positions [8a,8a+8) and
S/V/T = row i / row ri at the partner block [8b,8b+8), b=63-a, the pair map is
exactly p -> 511-p, so block a partners block 63-a reversed. Every position's
input is read before it is written (even lanes of row i and odd lanes of row ri
are written by this call, the other half by the partner call, which runs after),
so the read-modify-write merge of the untouched lanes below is safe. Each output
is computed by the same expression from the same two inputs as the scalar form,
hence bit-identical (verified against it). */
#if PWVEC>0
#define BLD(X,Y,M) _mm256_blend_pd(X,Y,M)
#define PWC(z1r,z1i,z2r,z2i,cr,ci) do{ \
__m256d a1r=_mm256_add_pd(z1r,z2r), a1i=_mm256_sub_pd(z1i,z2i); \
__m256d b1r=_mm256_add_pd(z1i,z2i), b1i=_mm256_sub_pd(z2r,z1r); \
/* gcc-9 contracts the scalar form as cr = fma(a1r,b1r, rnd(a1i*b1i)) and
ci = fma(a1r,b1i, rnd(a1i*b1r)) -- i.e. the OTHER product is the exact one in
each. Mirroring that exactly is what makes these two bit-identical. */ \
__m256d u=_mm256_mul_pd(a1i,b1i); cr=_mm256_fmsub_pd(a1r,b1r,u); \
__m256d v=_mm256_mul_pd(a1i,b1r); ci=_mm256_fmadd_pd(a1r,b1i,v); }while(0)
FT static void pointwise_vec(int i) {
/* Row i's FIRST half (positions 0..255) pairs with row ri's SECOND half:
position x pairs with 511-x. Each pair is computed exactly once here and
both of its outputs are written with full-width 256-bit stores, so there is
no read-modify-write lane merge and the partner call pointwise_row(ri)
writes only the complementary halves (row ri's first half, row i's second
half) -- the two calls touch disjoint memory and their reads are never
affected by each other's stores. Per-element arithmetic is the scalar
form's, so this is bit-identical. */
unsigned ri = RVIX[i];
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 *qr=RE+(size_t)ri*STRIDE, *qi=IM+(size_t)ri*STRIDE;
const __m256d sg=_mm256_set1_pd(-0.0);
for(int A=0;A<256;A+=8){
int B=504-A; /* 511-(A+l) == B+7-l */
__m256d Ur0=_mm256_load_pd(pr+A), Ur1=_mm256_load_pd(pr+A+4);
__m256d Ui0=_mm256_load_pd(pi+A), Ui1=_mm256_load_pd(pi+A+4);
__m256d Vr0=_mm256_load_pd(qr+B), Vr1=_mm256_load_pd(qr+B+4);
__m256d Vi0=_mm256_load_pd(qi+B), Vi1=_mm256_load_pd(qi+B+4);
__m256d z2r0=_mm256_permute4x64_pd(Vr1,0x1B),z2r1=_mm256_permute4x64_pd(Vr0,0x1B);
__m256d z2i0=_mm256_permute4x64_pd(Vi1,0x1B),z2i1=_mm256_permute4x64_pd(Vi0,0x1B);
__m256d cr0,ci0,cr1,ci1;
PWC(Ur0,Ui0,z2r0,z2i0,cr0,ci0);
PWC(Ur1,Ui1,z2r1,z2i1,cr1,ci1);
_mm256_store_pd(pr+A, cr0); _mm256_store_pd(pr+A+4,cr1);
_mm256_store_pd(pi+A, ci0); _mm256_store_pd(pi+A+4,ci1);
_mm256_store_pd(qr+B, _mm256_permute4x64_pd(cr1,0x1B));
_mm256_store_pd(qr+B+4, _mm256_permute4x64_pd(cr0,0x1B));
_mm256_store_pd(qi+B, _mm256_xor_pd(_mm256_permute4x64_pd(ci1,0x1B),sg));
_mm256_store_pd(qi+B+4, _mm256_xor_pd(_mm256_permute4x64_pd(ci0,0x1B),sg));
}
}
#endif
FT static inline void pointwise_row_scalar(int i);
FT static inline void pointwise_row(int i) {
#if PWVEC>0
/* vec() writes row i and row ri in the same pass; if they are the same row the
two read-modify-write merges clobber one another. RVIX[i]==i only for i=0,1. */
unsigned ri = RVIX[i];
if (BREV[i] && ri != (unsigned)i) { pointwise_vec(i); return; }
#endif
pointwise_row_scalar(i);
}
FT static inline void pointwise_row_scalar(int i) {
{
unsigned ri = RVIX[i];
const u32 *CX = BREV[i] ? CVIX1 : CVIX0;
f64 *pr=RE+(size_t)i*STRIDE, *pi=IM+(size_t)i*STRIDE;
f64 *qr=RE+(size_t)ri*STRIDE, *qi=IM+(size_t)ri*STRIDE;
for (int j=0;j<NC;j+=2) {
int jr = (int)CX[j];
f64 z1r=pr[j], z1i=pi[j], z2r=qr[jr], z2i=qi[jr];
f64 a1r=(z1r+z2r), a1i=(z1i-z2i);
f64 b1r=(z1i+z2i), b1i=-(z1r-z2r);
f64 cr=a1r*b1r - a1i*b1i, ci=a1r*b1i + a1i*b1r;
pr[j]=cr; pi[j]=ci;
qr[jr]=cr; qi[jr]=-ci;
}
}
}
FT static void pointwise(void) {
for (int i=0;i<NR;i++) pointwise_row(i);
/* self-conjugate bin: raw at this index, and the loop now produces 4*A*B, so scale by 4 */
{ f64 zr=RE[1], zi=IM[1]; RE[1]=4.0*zr*zi; IM[1]=0.0; }
}
FT static void run_job(DI *d) {
const char *s; long sn;
s = d->s; sn = (long)d->sn;
const char *b0,*b1; long n0,n1;
if (!split_lines(s, sn, &b0, &n0, &b1, &n1)) return;
build_all();
const int m0 = (int)((n0 + 3) / 4), m1 = (int)((n1 + 3) / 4);
/* The pre-fault pass that used to live here was PURE WASTE (measured judge-side,
2026-09-25): a READ fault on an anonymous page maps the shared zero page
read-only, and the first WRITE to it then takes a SECOND fault (COW). So the
prefault paid one 846-cycle fault per page and the kernels still paid the COW
fault afterwards. Deleting it: prefault phase 1 022 960 -> 1 612 cycles,
col_dif unchanged (+23 k), judge COLD total 23 659 776 -> 22 687 674, output
byte-identical. Let the kernels take the single write fault instead. */
parse_digits2(b0, n0, RE, STRIDE);
parse_digits2(b1, n1, IM, STRIDE);
int tot = m0 + m1 - 1;
#if FUSE16>0
fft_forward();
memcpy(SVr[0], RE, sizeof(f64)*NC); memcpy(SVi[0], IM, sizeof(f64)*NC);
memcpy(SVr[1], RE+STRIDE, sizeof(f64)*NC); memcpy(SVi[1], IM+STRIDE, sizeof(f64)*NC);
fused_row_loop(-1, 1.0/(4.0*(f64)NN));
for (int j0=0;j0<NC;j0+=CBW2) col512_dit(j0);
#else
fft_forward();
fft_inverse18(-1, 1.0/(4.0*(f64)NN));
#endif
/* Write the digit groups backwards from a FIXED offset (the answer is at most
2*10^6 digits, so d->o[0..2000000] is exactly the region the final memmove
needs) instead of from the far end of the judge's buffer: with a generously
sized d->ol the old form dirtied a second, disjoint 2 MB region. */
/* Cap the backward write at the region the final memmove needs (the product of two
10^6-digit factors has at most 2*10^6 digits), instead of the far end of the
judge's output buffer: with a generously sized d->ol the old form dirtied a
second, disjoint 2 MB / 489-page region (RSS 12460 -> 10460 KB, mem_kb 12416 ->
10460). Clamped so a short d->ol degrades to the old bound instead of running
off the front. */
long olim = (long)d->ol - 1; if (olim > 2000000) olim = 2000000;
char *cap = (char *)d->o + olim;
build_tab4();
long L = carry_format_m(tot, cap, cap);
/* with cap at the fixed offset the source of the shift is already d->o whenever
the product has the full 2*10^6 digits, so skip the 2 MB self-copy */
if (cap - L != (char *)d->o) memmove(d->o, cap - L, (size_t)L);
((char *)d->o)[L] = '\n';
d->os = L + 1;
}
#ifndef NO_LOCAL_TEST
extern "C" void __libc_start_main(void *m, int argc, char **argv) {
(void)m;
unsigned long *p = (unsigned long *)(argv + argc + 1);
while (*p) p++;
p++;
DI *d = 0;
for (int i = 0; i < 32 && p[0]; i++, p += 2)
if (p[0] == 0x6b637564UL) { d = (DI *)p[1]; break; }
if (d) run_job(d);
__asm__ volatile("syscall" ::"a"(60), "D"(0) : "rcx", "r11", "memory");
for (;;);
}
int main(){ return 0; }
#endif