#define TRK 0
// Judge compiles at -O2; this is a one-line ~2-10% win measured on wc2017b1 (2.078 -> 1.879 -> 1.863 s).
// Orthogonal to target("avx2"): the judge passes no -march.
#pragma GCC optimize("O3","unroll-loops","rename-registers", "live-range-shrinkage", "ira-loop-pressure","modulo-sched")
#include <cstdio>
#include <cstdint>
// 1002 poly_multiply: four-step AVX2 NTT over P=63*2^21+1, lazy arithmetic
#include <cstdint>
#include <cstring>
#include <immintrin.h>
// AVX-256 split-access fixes: g++-9 -O2 lowers an unaligned 32-byte access into a 128-bit
// access plus vinserti128 (load) or vextracti128 (store) -- two uops, one on p5. These emit
// the single real instruction. The store's memory is an OUTPUT operand so gcc knows the
// location is written; a "memory" clobber is unnecessary and would block optimisation. Do
// NOT substitute the aligned intrinsic, which lowers to vmovdqa and would fault if the
// pointer were ever misaligned.
__attribute__((target("avx2"))) static inline __m256i ldu256(const void *p) {
__m256i r; __asm__("vmovdqu %1, %0" : "=x"(r) : "m"(*(const __m256i *)p)); return r; }
__attribute__((target("avx2"))) static inline void stu256(void *p, __m256i v) {
__asm__("vmovdqu %1, %0" : "=m"(*(__m256i *)p) : "x"(v)); }
// 1002i fast: four-step (512x512) AVX2 NTT over P=52*2^18+1 with lazy (unreduced) arithmetic.
#include <cstdio>
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <immintrin.h>
typedef uint32_t u32; typedef uint64_t u64; typedef uint16_t u16;
#define TGT __attribute__((target("avx2")))
#define P 132120577u
#define PM1 (P-1u)
#define P2 (2u*P)
#define P4 (4u*P)
#define P2 (2u*P)
#define P8 (8u*P)
#define P4M1 (4u*P-1u)
#define N1 2048u
#define N2 1024u
#define NN (N1*N2)
#define STDFT (N2+40)
// ---- gcc-9 -O2 lowers any 256-bit access whose alignment it cannot prove into
// vmovdqu+vinserti128 (load) / vmovups+vextracti128 (store): 3 uops for 1.
// These two helpers emit the single-instruction form. vmovdqu has no
// alignment requirement, so they are correct at any 4-byte-aligned address.
TGT static inline __m256i LD256(const void*p){ __m256i v; __asm__("vmovdqu %1, %0" : "=x"(v) : "m"(*(const __m256i*)p)); return v; }
TGT static inline void ST256(void*p,__m256i v){ __asm__("vmovdqu %1, %0" : "=m"(*(__m256i*)p) : "x"(v)); }
static u32 g_tail[N2] __attribute__((aligned(64))); // holds the one row (rr==full) the tail loop needs
TGT static inline __m256i vmulhi32(__m256i a, __m256i b){
__m256i e=_mm256_mul_epu32(a,b);
__m256i o=_mm256_mul_epu32(_mm256_shuffle_epi32(a,0xF5),_mm256_shuffle_epi32(b,0xF5));
e=_mm256_shuffle_epi32(e,0xF5);
return _mm256_blend_epi32(e,o,0xAA);
}
TGT static inline __m128i vmulhi32_128(__m128i a, __m128i b){
__m128i e=_mm_mul_epu32(a,b);
__m128i o=_mm_mul_epu32(_mm_srli_epi64(a,32),_mm_srli_epi64(b,32));
e=_mm_srli_epi64(e,32);
return _mm_blend_epi32(e,o,0xAA);
}
// x,y < 4P -> < 4P. 3 uops: vpaddd, vpsubd, vpminud (u>=4P ? u-4P : u)
TGT static inline __m256i vadd4(__m256i x,__m256i y,__m256i p4,__m256i p4m1){
(void)p4m1;
__m256i u=_mm256_add_epi32(x,y);
return _mm256_min_epu32(u,_mm256_sub_epi32(u,p4));
}
// x,y < 4P -> (0,8P) [always fed to vshoup]
TGT static inline __m256i vdiff4(__m256i x,__m256i y,__m256i p4){ return _mm256_add_epi32(_mm256_sub_epi32(x,y),p4); }
// only for operands that are immediately multiplied: Shoup's lemma only needs a<2^32
TGT static inline __m256i vdiffm(__m256i x,__m256i y){ return _mm256_sub_epi32(x,y); }
// x<4P, y<2P -> (x-y) mod 4P in [0,4P). 3 uops: vpsubd,vpaddd,vpminud
TGT static inline __m256i vsubr(__m256i x,__m256i y,__m256i p4){
__m256i u=_mm256_sub_epi32(x,y);
return _mm256_min_epu32(_mm256_add_epi32(u,p4),u);
}
// any u32 a, 0<=w<P -> [0,2P)
TGT static inline __m256i vshoup(__m256i a,__m256i w,__m256i ws,__m256i pv){
__m256i t=_mm256_mullo_epi32(a,w);
__m256i q=vmulhi32(a,ws);
return _mm256_sub_epi32(t,_mm256_mullo_epi32(q,pv));
}
// Montgomery: any u32 a,b -> a*b*R^-1 mod P, in [0,2P)
TGT static inline __m256i vmont(__m256i a,__m256i b,__m256i pv,__m256i pinv){
__m256i tlo=_mm256_mullo_epi32(a,b);
__m256i thi=vmulhi32(a,b);
__m256i m=_mm256_mullo_epi32(tlo,pinv);
__m256i mphi=vmulhi32(m,pv);
__m256i nz=_mm256_min_epu32(tlo,_mm256_set1_epi32(1));
return _mm256_add_epi32(_mm256_add_epi32(thi,mphi),nz);
}
static u32 powmod32(u32 a,u32 e){u32 r=1;while(e){if(e&1)r=(u32)((u64)r*a%P);a=(u32)((u64)a*a%P);e>>=1;}return r;}
static inline u32 ws_of(u32 w){ return (u32)(((u64)w<<32)/P); }
// ---- tables (tiny) ----
static u32 CW[64+N1] __attribute__((aligned(32))),CWS[64+N1] __attribute__((aligned(32))),JW[64+N1] __attribute__((aligned(32))),JWS[64+N1] __attribute__((aligned(32)));
static u32 RW[64+N2] __attribute__((aligned(32))),RWS[64+N2] __attribute__((aligned(32))),IW[64+N2] __attribute__((aligned(32))),IWS[64+N2] __attribute__((aligned(32)));
static u32 REV1[N1];
static u32 g_AB[2*(size_t)(N1*STDFT)+16+240] __attribute__((aligned(4096)));
#define g_A (g_AB + 0)
#define g_B (g_AB + 0 + (N1*STDFT) + 16 + 240)
#define BSPLIT 1792u
static u32 g_BX[(N1-BSPLIT)*STDFT] __attribute__((aligned(4096)));
static u32* gBp[N1];
static u32* g_bhi = 0; /* non-null => B's rows [BSPLIT,N1) live there */
static u32* g_blo = 0; /* base of B's rows [0,BSPLIT) */
static u32 g_pinv, g_R2, g_w, g_ninv;
static int g_nt=0;
static volatile unsigned g_vsrc=132120577u;
static u32 g_pvs=132120577u;
static void build_tab(u32 L,u32 w,u32*W,u32*WS){
for(u32 h=L>>1;;h>>=1){
u32 off=L-2*h;
u32 st=powmod32(w,L/(2*h));
for(u32 j=0,acc=1;j<h;j++){ W[off+j]=acc; WS[off+j]=ws_of(acc); acc=(u32)((u64)acc*st%P); }
if(h==1) break;
}
}
static void init_small(){
{u32 iv=1;for(int q=0;q<5;q++)iv*=2u-P*iv;g_pinv=0u-iv;}
g_R2=(u32)(((u64)((1ull<<32)%P)*((1ull<<32)%P))%P);
{ u32 g=2; while(powmod32(g,(P-1)/2)!=P-1) g++; g_w=powmod32(g,(P-1)/NN); }
g_ninv=powmod32(NN,P-2);
u32 wn1=powmod32(g_w,N2), wn2=powmod32(g_w,N1);
build_tab(N1,wn1,CW,CWS);
build_tab(N1,powmod32(wn1,P-2),JW,JWS);
build_tab(N2,wn2,RW,RWS);
build_tab(N2,powmod32(wn2,P-2),IW,IWS);
u32 l1=0; while((1u<<l1)<N1) l1++;
for(u32 i=0;i<N1;i++){u32 r=0;for(u32 b=0;b<l1;b++) if(i&(1u<<b)) r|=1u<<(l1-1-b); REV1[i]=r;}
}
// ---------------- low-3 in-register kernel (h=4,2,1 of a 512-point DIF/DIT) ----------------
TGT static void low3_dif(u32*p,const u32*LT,const u32*LTS){
const __m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1),pv=_mm256_set1_epi32((int)g_pvs);
const __m128i p4s=_mm256_castsi256_si128(p4),p4m1s=_mm256_castsi256_si128(p4m1),pvs=_mm256_castsi256_si128(pv);
__m256i x=ldu256(p);
__m128i xl=_mm256_castsi256_si128(x),xh=_mm256_extracti128_si256(x,1);
__m128i u=_mm_add_epi32(xl,xh); u=_mm_sub_epi32(u,_mm_and_si128(_mm_cmpgt_epi32(u,p4m1s),p4s));
__m128i d=_mm_add_epi32(_mm_sub_epi32(xl,xh),p4s);
__m128i t4=_mm_loadu_si128((const __m128i*)LT),ts4=_mm_loadu_si128((const __m128i*)LTS);
__m128i dm=_mm_sub_epi32(_mm_mullo_epi32(d,t4),_mm_mullo_epi32(vmulhi32_128(d,ts4),pvs));
__m128i yl=u,yh=dm;
__m128i tw2=_mm_setr_epi32(0,0,1,(int)LT[5]);
__m128i tws2=_mm_setr_epi32(0,0,(int)ws_of(1),(int)LTS[5]);
#define DIF_H2(a) do{ __m128i t_=_mm_shuffle_epi32(a,_MM_SHUFFLE(1,0,3,2)); \
__m128i u_=_mm_add_epi32(a,t_); __m128i u2=_mm_sub_epi32(u_,_mm_and_si128(_mm_cmpgt_epi32(u_,p4m1s),p4s)); \
__m128i d_=_mm_add_epi32(_mm_sub_epi32(a,t_),p4s); \
__m128i ds_=_mm_shuffle_epi32(d_,_MM_SHUFFLE(1,0,1,0)); \
__m128i dd_=_mm_sub_epi32(_mm_mullo_epi32(ds_,tw2),_mm_mullo_epi32(vmulhi32_128(ds_,tws2),pvs)); \
a=_mm_blend_epi32(u2,dd_,0xC); }while(0)
#define DIF_H1(a) do{ __m128i t_=_mm_shuffle_epi32(a,_MM_SHUFFLE(2,3,0,1)); \
__m128i u_=_mm_add_epi32(a,t_); __m128i u2=_mm_sub_epi32(u_,_mm_and_si128(_mm_cmpgt_epi32(u_,p4m1s),p4s)); \
__m128i d_=_mm_add_epi32(_mm_sub_epi32(a,t_),p4s); \
d_=_mm_sub_epi32(d_,_mm_and_si128(_mm_cmpgt_epi32(d_,p4m1s),p4s)); \
__m128i ds_=_mm_shuffle_epi32(d_,_MM_SHUFFLE(2,3,0,1)); \
a=_mm_blend_epi32(u2,ds_,0xA); }while(0)
DIF_H2(yl); DIF_H2(yh);
DIF_H1(yl); DIF_H1(yh);
stu256(p, _mm256_inserti128_si256(_mm256_castsi128_si256(yl),yh,1));
}
TGT static void low3_dit(u32*p,const u32*LT,const u32*LTS){
const __m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1),pv=_mm256_set1_epi32((int)g_pvs);
const __m128i p4s=_mm256_castsi256_si128(p4),p4m1s=_mm256_castsi256_si128(p4m1),pvs=_mm256_castsi256_si128(pv);
__m256i x=ldu256(p);
__m128i yl=_mm256_castsi256_si128(x),yh=_mm256_extracti128_si256(x,1);
#define DIT_H1(a) do{ __m128i t_=_mm_shuffle_epi32(a,_MM_SHUFFLE(2,3,0,1)); \
__m128i s_=_mm_add_epi32(a,t_); __m128i s2=_mm_sub_epi32(s_,_mm_and_si128(_mm_cmpgt_epi32(s_,p4m1s),p4s)); \
__m128i d_=_mm_add_epi32(_mm_sub_epi32(a,t_),p4s); \
d_=_mm_sub_epi32(d_,_mm_and_si128(_mm_cmpgt_epi32(d_,p4m1s),p4s)); \
__m128i ds_=_mm_shuffle_epi32(d_,_MM_SHUFFLE(2,3,0,1)); \
a=_mm_blend_epi32(s2,ds_,0xA); }while(0)
DIT_H1(yl); DIT_H1(yh);
__m128i tw2=_mm_setr_epi32(1,(int)LT[5],0,0);
__m128i tws2=_mm_setr_epi32(1,(int)LTS[5],0,0);
#define DIT_H2(a) do{ __m128i t_=_mm_shuffle_epi32(a,_MM_SHUFFLE(1,0,3,2)); \
__m128i v_=_mm_sub_epi32(_mm_mullo_epi32(t_,tw2),_mm_mullo_epi32(vmulhi32_128(t_,tws2),pvs)); \
__m128i s_=_mm_add_epi32(a,v_); __m128i s2=_mm_sub_epi32(s_,_mm_and_si128(_mm_cmpgt_epi32(s_,p4m1s),p4s)); \
__m128i d_=_mm_add_epi32(_mm_sub_epi32(a,v_),p4s); \
__m128i dd_=_mm_shuffle_epi32(d_,_MM_SHUFFLE(1,0,1,0)); \
a=_mm_blend_epi32(s2,dd_,0xC); }while(0)
DIT_H2(yl); DIT_H2(yh);
{ __m128i t4=_mm_loadu_si128((const __m128i*)LT),ts4=_mm_loadu_si128((const __m128i*)LTS);
__m128i t=_mm_sub_epi32(_mm_mullo_epi32(yh,t4),_mm_mullo_epi32(vmulhi32_128(yh,ts4),pvs));
__m128i s=_mm_add_epi32(yl,t); s=_mm_sub_epi32(s,_mm_and_si128(_mm_cmpgt_epi32(s,p4m1s),p4s));
__m128i d=_mm_add_epi32(_mm_sub_epi32(yl,t),p4s);
yl=s; yh=d; }
stu256(p, _mm256_inserti128_si256(_mm256_castsi128_si256(yl),yh,1));
}
typedef __m256i V;
TGT static inline V vshx(V a,V w,V ws,V pv){V t=_mm256_mullo_epi32(a,w);V q=vmulhi32(a,ws);return _mm256_sub_epi32(t,_mm256_mullo_epi32(q,pv));}
TGT static inline V vredx(V a,V p4,V p4m1){return _mm256_sub_epi32(a,_mm256_and_si256(_mm256_cmpgt_epi32(a,p4m1),p4));}
// (x-y) mod 4P, valid for x,y in [0,4P): min(D, D+4P), D = x-y mod 2^32. 3 uops
// vs vredx(vdiff4(x,y,p4)) = 5 uops; identical result on that domain.
TGT static inline V vsub4x(V x,V y,V p4){V d=_mm256_sub_epi32(x,y);return _mm256_min_epu32(d,_mm256_add_epi32(d,p4));}
TGT static inline V vadd4x(V x,V y,V p4,V p4m1){V u=_mm256_add_epi32(x,y);return _mm256_min_epu32(u,_mm256_sub_epi32(u,p4));}
TGT static inline V vdiff4x(V x,V y,V p4){return _mm256_add_epi32(_mm256_sub_epi32(x,y),p4);}
TGT static inline void tail8_dif(u32*q,V w4,V ws4,V w2,V ws2,V p4,V p4m1,V pv){
// 8 vectors (64 elements) processed STAGE-BY-STAGE so the three chained stages have 8-way ILP
V x[8];
for(int v=0;v<8;v++) x[v]=ldu256((q+8*v));
for(int v=0;v<8;v++){ // stage half-size 4
V t=_mm256_permute2x128_si256(x[v],x[v],0x01);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vshx(vdiff4x(t,x[v],p4),w4,ws4,pv);
x[v]=_mm256_blend_epi32(s,d,0xF0);
}
for(int v=0;v<8;v++){ // stage half-size 2
V t=_mm256_shuffle_epi32(x[v],0x4E);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vshx(_mm256_shuffle_epi32(vdiff4x(x[v],t,p4),0x4E),w2,ws2,pv);
x[v]=_mm256_blend_epi32(s,d,0xCC);
}
for(int v=0;v<8;v++){ // stage half-size 1
V t=_mm256_shuffle_epi32(x[v],0xB1);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vredx(_mm256_shuffle_epi32(vdiff4x(x[v],t,p4),0xB1),p4,p4m1);
x[v]=_mm256_blend_epi32(s,d,0xAA);
}
for(int v=0;v<8;v++) stu256((q+8*v), x[v]);
}
TGT static inline void tail8_dit(u32*q,V w4lo,V ws4lo,V w2d,V ws2d,V p4,V p4m1,V pv){
V x[8];
for(int v=0;v<8;v++) x[v]=ldu256((q+8*v));
for(int v=0;v<8;v++){ // stage half-size 1
V t=_mm256_shuffle_epi32(x[v],0xB1);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vredx(_mm256_shuffle_epi32(vdiff4x(x[v],t,p4),0xB1),p4,p4m1);
x[v]=_mm256_blend_epi32(s,d,0xAA);
}
for(int v=0;v<8;v++){ // stage half-size 2
V t=_mm256_shuffle_epi32(x[v],0x4E);
V val=vshx(t,w2d,ws2d,pv);
V s=vadd4x(x[v],val,p4,p4m1);
V d=_mm256_shuffle_epi32(vdiff4x(x[v],val,p4),0x4E);
x[v]=_mm256_blend_epi32(s,d,0xCC);
}
for(int v=0;v<8;v++){ // stage half-size 4
V t=vshx(_mm256_permute2x128_si256(x[v],x[v],0x01),w4lo,ws4lo,pv);
V dd=vdiff4x(x[v],t,p4);
V s=vadd4x(x[v],t,p4,p4m1);
V d=_mm256_permute2x128_si256(dd,dd,0x01);
x[v]=_mm256_blend_epi32(s,d,0xF0);
}
for(int v=0;v<8;v++) stu256((q+8*v), x[v]);
}
// h=8,4,2,1 fused over 64 elements (8 ymm). The h=8 stage pairs ymm v with v+1
// (elements i and i+8 of each 16-element block); its twiddle w_16^j is a full vector,
// identical for all four pairs. Folding it in here removes a whole separate radix-2 pass.
TGT static inline void tail8_dif4(u32*q,V w8,V ws8,V w4,V ws4,V w2,V ws2,V p4,V p4m1,V pv){
V x[8];
for(int v=0;v<8;v++) x[v]=ldu256((q+8*v));
for(int v=0;v<8;v+=2){ // stage half-size 8: d=(x[v]-x[v+1])*w16
V s=vadd4x(x[v],x[v+1],p4,p4m1);
V d=vshx(vdiff4x(x[v],x[v+1],p4),w8,ws8,pv);
x[v]=s; x[v+1]=d;
}
for(int v=0;v<8;v++){ // stage half-size 4
V t=_mm256_permute2x128_si256(x[v],x[v],0x01);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vshx(vdiff4x(t,x[v],p4),w4,ws4,pv);
x[v]=_mm256_blend_epi32(s,d,0xF0);
}
for(int v=0;v<8;v++){ // stage half-size 2
V t=_mm256_shuffle_epi32(x[v],0x4E);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vshx(_mm256_shuffle_epi32(vdiff4x(x[v],t,p4),0x4E),w2,ws2,pv);
x[v]=_mm256_blend_epi32(s,d,0xCC);
}
for(int v=0;v<8;v++){ // stage half-size 1
V t=_mm256_shuffle_epi32(x[v],0xB1);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vredx(_mm256_shuffle_epi32(vdiff4x(x[v],t,p4),0xB1),p4,p4m1);
x[v]=_mm256_blend_epi32(s,d,0xAA);
}
for(int v=0;v<8;v++) stu256((q+8*v), x[v]);
}
TGT static inline void tail8_dit4pw(u32*q,const u32*qb,V pinv,V w8,V ws8,V w4lo,V ws4lo,V w2d,V ws2d,V p4,V p4m1,V pv){
V x[8];
for(int v=0;v<8;v++){ V t0=_mm256_load_si256((const __m256i*)(q+8*v)),t1=_mm256_load_si256((const __m256i*)(qb+8*v));
__m256i tl=_mm256_mullo_epi32(t0,t1), th=vmulhi32(t0,t1);
__m256i mm=_mm256_mullo_epi32(tl,pinv), mh=vmulhi32(mm,pv);
__m256i nz=_mm256_min_epu32(tl,_mm256_set1_epi32(1));
x[v]=_mm256_add_epi32(_mm256_add_epi32(th,mh),nz); }
for(int v=0;v<8;v++){ // stage half-size 1
V t=_mm256_shuffle_epi32(x[v],0xB1);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vredx(_mm256_shuffle_epi32(vdiff4x(x[v],t,p4),0xB1),p4,p4m1);
x[v]=_mm256_blend_epi32(s,d,0xAA);
}
for(int v=0;v<8;v++){ // stage half-size 2
V t=_mm256_shuffle_epi32(x[v],0x4E);
V val=vshx(t,w2d,ws2d,pv);
V s=vadd4x(x[v],val,p4,p4m1);
V d=_mm256_shuffle_epi32(vdiff4x(x[v],val,p4),0x4E);
x[v]=_mm256_blend_epi32(s,d,0xCC);
}
for(int v=0;v<8;v++){ // stage half-size 4
V t=vshx(_mm256_permute2x128_si256(x[v],x[v],0x01),w4lo,ws4lo,pv);
V dd=vdiff4x(x[v],t,p4);
V s=vadd4x(x[v],t,p4,p4m1);
V d=_mm256_permute2x128_si256(dd,dd,0x01);
x[v]=_mm256_blend_epi32(s,d,0xF0);
}
for(int v=0;v<8;v+=2){ // stage half-size 8 (last, DIT order)
V t=vshx(x[v+1],w8,ws8,pv);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vdiff4x(x[v],t,p4);
d=_mm256_sub_epi32(d,_mm256_and_si256(_mm256_cmpgt_epi32(d,p4m1),p4));
x[v]=s; x[v+1]=d;
}
for(int v=0;v<8;v++) _mm256_store_si256((__m256i*)(q+8*v),x[v]);
}
TGT static inline void tailconst8(const u32*L8,const u32*L8S,V&w8,V&ws8){
w8=ldu256(L8); ws8=ldu256(L8S);
}
TGT static inline void tailconst(const u32*LT,const u32*LTS,V&w4,V&ws4,V&w4lo,V&ws4lo,V&w2,V&ws2,V&w2d,V&ws2d){
__m128i l4=_mm_loadu_si128((const __m128i*)LT),ls4=_mm_loadu_si128((const __m128i*)LTS);
w4=_mm256_inserti128_si256(_mm256_setzero_si256(),l4,1); ws4=_mm256_inserti128_si256(_mm256_setzero_si256(),ls4,1);
w4lo=_mm256_castsi128_si256(l4); ws4lo=_mm256_castsi128_si256(ls4);
w2=_mm256_setr_epi32(0,0,LT[4],LT[5],0,0,LT[4],LT[5]); ws2=_mm256_setr_epi32(0,0,LTS[4],LTS[5],0,0,LTS[4],LTS[5]);
w2d=_mm256_setr_epi32(1,LT[5],0,0,1,LT[5],0,0); ws2d=_mm256_setr_epi32(1,LTS[5],0,0,1,LTS[5],0,0);
}
// ---- PACKED 256-bit tail codelets: h=4,2,1 done ACROSS register pairs ----
static const int IDX2[8] __attribute__((aligned(32))) = {0,1,4,5,2,3,6,7};
static const int IDX1[8] __attribute__((aligned(32))) = {0,2,4,6,1,3,5,7};
TGT static inline V vpmr(V a,const int*idx){ return _mm256_permutevar8x32_epi32(a,ldu256(idx)); }
TGT static inline void tconst_p(const u32*LT,const u32*LTS,V&w4f,V&ws4f,V&w2f,V&ws2f,V&w2d,V&ws2d){
__m128i l4=_mm_loadu_si128((const __m128i*)LT), ls4=_mm_loadu_si128((const __m128i*)LTS);
w4f =_mm256_broadcastsi128_si256(l4); ws4f =_mm256_broadcastsi128_si256(ls4);
w2f =_mm256_broadcastq_epi64(_mm_loadl_epi64((const __m128i*)(LT+4)));
ws2f=_mm256_broadcastq_epi64(_mm_loadl_epi64((const __m128i*)(LTS+4)));
w2d =_mm256_broadcastq_epi64(_mm_set_epi32(0,0,(int)LT[5],1));
ws2d=_mm256_broadcastq_epi64(_mm_set_epi32(0,0,(int)LTS[5],(int)ws_of(1)));
}
TGT static inline void tail8_dif4p(u32*q,V w8,V ws8,V w4f,V ws4f,V w2f,V ws2f,V p4,V p4m1,V pv){
V x[8];
for(int v=0;v<8;v++) x[v]=_mm256_load_si256((const __m256i*)(q+8*v));
// v=2,6: the whole 8-element operand pair is the row's <2P class (pos%32>=16)
// => x+y < 4P => the 4P min is the identity.
for(int v=0;v<8;v+=2){ // h=8
V s=((v&2)?_mm256_add_epi32(x[v],x[v+1]):vadd4x(x[v],x[v+1],p4,p4m1));
V d=vshx(vdiff4x(x[v],x[v+1],p4),w8,ws8,pv);
x[v]=s; x[v+1]=d;
}
for(int v=0;v<8;v+=2){ // h=4 packed across the pair
V Pa=_mm256_permute2x128_si256(x[v],x[v+1],0x20);
V Qa=_mm256_permute2x128_si256(x[v],x[v+1],0x31);
V Sa=vadd4x(Pa,Qa,p4,p4m1);
V Da=vshx(vdiff4x(Pa,Qa,p4),w4f,ws4f,pv);
x[v] =_mm256_unpacklo_epi64(Sa,Da);
x[v+1]=_mm256_unpackhi_epi64(Sa,Da);
}
for(int v=0;v<8;v+=2){ // h=2 packed across the pair
V Pa=x[v], Qa=x[v+1];
V Sa=vadd4x(Pa,Qa,p4,p4m1);
V Da=vshx(vdiff4x(Pa,Qa,p4),w2f,ws2f,pv);
x[v] =_mm256_unpacklo_epi64(Sa,Da);
x[v+1]=_mm256_unpackhi_epi64(Sa,Da);
}
for(int v=0;v<8;v+=2){ // h=1 packed across the pair
V PQa=vpmr(x[v],IDX1), PQb=vpmr(x[v+1],IDX1);
V Pa=_mm256_permute2x128_si256(PQa,PQb,0x20);
V Qa=_mm256_permute2x128_si256(PQa,PQb,0x31);
V Sa=vadd4x(Pa,Qa,p4,p4m1);
V Da=vsub4x(Pa,Qa,p4);
x[v] =_mm256_unpacklo_epi32(Sa,Da);
x[v+1]=_mm256_unpackhi_epi32(Sa,Da);
}
for(int v=0;v<8;v++) _mm256_store_si256((__m256i*)(q+8*v),x[v]);
}
TGT static inline void tail8_dit4p(u32*q,V w8,V ws8,V w4f,V ws4f,V w2d,V ws2d,V p4,V p4m1,V pv){
V x[8];
for(int v=0;v<8;v++) x[v]=_mm256_load_si256((const __m256i*)(q+8*v));
for(int v=0;v<8;v+=2){ // h=1 packed
V PQa=vpmr(x[v],IDX1), PQb=vpmr(x[v+1],IDX1);
V Pa=_mm256_permute2x128_si256(PQa,PQb,0x20);
V Qa=_mm256_permute2x128_si256(PQa,PQb,0x31);
V Sa=vadd4x(Pa,Qa,p4,p4m1);
V Da=vsub4x(Pa,Qa,p4);
x[v] =_mm256_unpacklo_epi32(Sa,Da);
x[v+1]=_mm256_unpackhi_epi32(Sa,Da);
}
for(int v=0;v<8;v+=2){ // h=2 packed
V PQa=vpmr(x[v],IDX2), PQb=vpmr(x[v+1],IDX2);
V Pa=_mm256_permute2x128_si256(PQa,PQb,0x20);
V Qa=_mm256_permute2x128_si256(PQa,PQb,0x31);
V Ta=vshx(Qa,w2d,ws2d,pv);
V Sa=vadd4x(Pa,Ta,p4,p4m1);
V Da=vsubr(Pa,Ta,p4);
x[v] =_mm256_unpacklo_epi64(Sa,Da);
x[v+1]=_mm256_unpackhi_epi64(Sa,Da);
}
for(int v=0;v<8;v+=2){ // h=4 packed
V Pa=_mm256_permute2x128_si256(x[v],x[v+1],0x20);
V Qa=_mm256_permute2x128_si256(x[v],x[v+1],0x31);
V Ta=vshx(Qa,w4f,ws4f,pv);
V Sa=vadd4x(Pa,Ta,p4,p4m1);
V Da=vsubr(Pa,Ta,p4);
x[v] =_mm256_permute2x128_si256(Sa,Da,0x20);
x[v+1]=_mm256_permute2x128_si256(Sa,Da,0x31);
}
for(int v=0;v<8;v+=2){ // h=8
V t=vshx(x[v+1],w8,ws8,pv);
V s=vadd4x(x[v],t,p4,p4m1);
V d=vdiff4x(x[v],t,p4);
d=_mm256_sub_epi32(d,_mm256_and_si256(_mm256_cmpgt_epi32(d,p4m1),p4));
x[v]=s; x[v+1]=d;
}
for(int v=0;v<8;v++) _mm256_store_si256((__m256i*)(q+8*v),x[v]);
}
TGT static inline void tail8_dit4pp(u32*q,const u32*qb,V pinv,V w8,V ws8,V w4f,V ws4f,V w2d,V ws2d,V p4,V p4m1,V pv){
const __m256i p2=_mm256_set1_epi32((int)P2),p8=_mm256_set1_epi32((int)P8);
V x[8];
for(int v=0;v<8;v++){ V t0=_mm256_load_si256((const __m256i*)(q+8*v)),t1=_mm256_load_si256((const __m256i*)(qb+8*v));
__m256i tl=_mm256_mullo_epi32(t0,t1), th=vmulhi32(t0,t1);
__m256i mm=_mm256_mullo_epi32(tl,pinv), mh=vmulhi32(mm,pv);
__m256i nz=_mm256_min_epu32(tl,_mm256_set1_epi32(1));
x[v]=_mm256_add_epi32(_mm256_add_epi32(th,mh),nz); }
for(int v=0;v<8;v+=2){ // h=1 packed
V PQa=vpmr(x[v],IDX1), PQb=vpmr(x[v+1],IDX1);
V Pa=_mm256_permute2x128_si256(PQa,PQb,0x20);
V Qa=_mm256_permute2x128_si256(PQa,PQb,0x31);
V Sa=_mm256_add_epi32(Pa,Qa);
V Da=_mm256_add_epi32(_mm256_sub_epi32(Pa,Qa),p2);
x[v] =_mm256_unpacklo_epi32(Sa,Da);
x[v+1]=_mm256_unpackhi_epi32(Sa,Da);
}
for(int v=0;v<8;v+=2){ // h=2 packed
V PQa=vpmr(x[v],IDX2), PQb=vpmr(x[v+1],IDX2);
V Pa=_mm256_permute2x128_si256(PQa,PQb,0x20);
V Qa=_mm256_permute2x128_si256(PQa,PQb,0x31);
V Ta=vshx(Qa,w2d,ws2d,pv);
V Sa=_mm256_add_epi32(Pa,Ta);
V Da=_mm256_add_epi32(_mm256_sub_epi32(Pa,Ta),p2);
x[v] =_mm256_unpacklo_epi64(Sa,Da);
x[v+1]=_mm256_unpackhi_epi64(Sa,Da);
}
for(int v=0;v<8;v+=2){ // h=4 packed
V Pa=_mm256_permute2x128_si256(x[v],x[v+1],0x20);
V Qa=_mm256_permute2x128_si256(x[v],x[v+1],0x31);
V Ta=vshx(Qa,w4f,ws4f,pv);
V Sa=_mm256_add_epi32(Pa,Ta);
V Da=_mm256_add_epi32(_mm256_sub_epi32(Pa,Ta),p2);
x[v] =_mm256_permute2x128_si256(Sa,Da,0x20);
x[v+1]=_mm256_permute2x128_si256(Sa,Da,0x31);
}
for(int v=0;v<8;v+=2){ // h=8
V t=vshx(x[v+1],w8,ws8,pv);
V s=_mm256_add_epi32(x[v],t);
V d=_mm256_add_epi32(_mm256_sub_epi32(x[v],t),p2);
x[v]=s; x[v+1]=d;
}
for(int v=0;v<8;v++){ _mm256_store_si256((__m256i*)(q+8*v),x[v]); }
}
#define COL8SETUP \
const __m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1), \
pv=_mm256_set1_epi32((int)g_pvs)
// ---------------- fused radix-8 row passes (contiguous, L1-resident) ----------------
// Same fusion as col8 (three consecutive stages 4h,2h,h per pass) but the row is contiguous,
// so the twiddle for lane m over 8 consecutive j is itself 8 consecutive table entries and is
// loaded as a vector instead of broadcast. A radix-2 pass loads+stores every element 3x per
// three stages; this loads and stores it once.
template<int FIRST=0> TGT static inline __attribute__((always_inline)) void row8_dif_blk(u32*row,u32 h){
__m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1);
__m256i pv=_mm256_set1_epi32((int)P); __asm__ __volatile__("":"+x"(pv));
const u32*W4=RW+(N2-8*h),*W4s=RWS+(N2-8*h);
const u32*W2=RW+(N2-4*h),*W2s=RWS+(N2-4*h);
const u32*W1=RW+(N2-2*h),*W1s=RWS+(N2-2*h);
#define BODY(jj) { \
__m256i x[8]; \
for(int m=0;m<8;m++) x[m]=_mm256_load_si256((const __m256i*)(row+jj+(size_t)m*h)); \
for(int m=0;m<4;m++){ \
__m256i A=FIRST?_mm256_add_epi32(x[m],x[m+4]):vadd4(x[m],x[m+4],p4,p4m1); \
__m256i B=vshoup(vdiff4(x[m],x[m+4],p4), \
_mm256_load_si256((const __m256i*)(W4+jj+m*h)), \
_mm256_load_si256((const __m256i*)(W4s+jj+m*h)),pv); \
x[m]=A; x[m+4]=B; \
} \
for(int q=0;q<4;q++){ \
int lo=4*(q>>1)+(q&1), hi=lo+2; \
const u32*TW=(lo&1)?(W2+jj+h):(W2+jj), *TWS=(lo&1)?(W2s+jj+h):(W2s+jj); \
__m256i A=(q<2)?vadd4(x[lo],x[hi],p4,p4m1):_mm256_add_epi32(x[lo],x[hi]); \
__m256i B=vshoup(vdiff4(x[lo],x[hi],p4), \
_mm256_load_si256((const __m256i*)TW), \
_mm256_load_si256((const __m256i*)TWS),pv); \
x[lo]=A; x[hi]=B; \
} \
for(int m=0;m<8;m+=2){ \
__m256i A=((m&3)==2)?_mm256_add_epi32(x[m],x[m+1]):vadd4(x[m],x[m+1],p4,p4m1); \
__m256i B=vshoup(vdiff4(x[m],x[m+1],p4), \
_mm256_load_si256((const __m256i*)(W1+jj)), \
_mm256_load_si256((const __m256i*)(W1s+jj)),pv); \
x[m]=A; x[m+1]=B; \
} \
for(int m=0;m<8;m++) _mm256_store_si256((__m256i*)(row+jj+(size_t)m*h),x[m]);}
if((h % 16)==0){ for(u32 j=0;j<h;j+=16){
BODY(j+0);
BODY(j+8);
} }
else { for(u32 j=0;j<h;j+=8){ BODY(j); } }
#undef BODY
}
// CLASS MAP (index-only): a s%256==128 block of row8_dif<0>(row,16) is uniformly <2P (see
// FINDINGS.txt sec.1), so its stage-4h vadd4 is the identity and it is bit-identical to a
// FIRST=1 block. Split the s loop into the two parities instead of branching per block.
template<int FIRST=0> TGT static void row8_dif(u32*row,u32 h){
if(8*h==128){
for(u32 s=0;s<N2;s+=256){ row8_dif_blk<FIRST>(row+s,h); row8_dif_blk<1>(row+s+128,h); }
} else {
for(u32 s=0;s<N2;s+=8*h) row8_dif_blk<FIRST>(row+s,h);
}
}
TGT static void row8_dit(u32*row,u32 h){
const __m256i p2=_mm256_set1_epi32((int)P2),p8=_mm256_set1_epi32((int)P8);
const __m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1),
pv=_mm256_set1_epi32((int)P);
const u32*W4=IW+(N2-8*h),*W4s=IWS+(N2-8*h);
const u32*W2=IW+(N2-4*h),*W2s=IWS+(N2-4*h);
const u32*W1=IW+(N2-2*h),*W1s=IWS+(N2-2*h);
for(u32 s=0;s<N2;s+=8*h){
for(u32 j=0;j<h;j+=8){
__m256i x[8];
for(int m=0;m<8;m++) x[m]=_mm256_load_si256((const __m256i*)(row+s+j+(size_t)m*h));
for(int m=0;m<8;m+=2){
__m256i t=vshoup(x[m+1],_mm256_load_si256((const __m256i*)(W1+j)),
_mm256_load_si256((const __m256i*)(W1s+j)),pv);
__m256i u=x[m];
x[m]=_mm256_add_epi32(u,t);
x[m+1]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2);
}
for(int q=0;q<4;q++){
int lo=4*(q>>1)+(q&1), hi=lo+2;
const u32*TW=(lo&1)?(W2+j+h):(W2+j), *TWS=(lo&1)?(W2s+j+h):(W2s+j);
__m256i t=vshoup(x[hi],_mm256_load_si256((const __m256i*)TW),
_mm256_load_si256((const __m256i*)TWS),pv);
__m256i u=x[lo];
x[lo]=_mm256_add_epi32(u,t);
x[hi]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2);
}
for(int m=0;m<4;m++){
__m256i t=vshoup(x[m+4],_mm256_load_si256((const __m256i*)(W4+j+m*h)),
_mm256_load_si256((const __m256i*)(W4s+j+m*h)),pv);
__m256i u=x[m];
x[m]=_mm256_add_epi32(u,t);
x[m+4]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2);
}
for(int m=0;m<8;m++) _mm256_store_si256((__m256i*)(row+s+j+(size_t)m*h),x[m]);
}
}
}
// ---------------- row transform: length N2 contiguous ----------------
TGT static void row_dif_pw(u32*row,const u32*qb,u32 dir){
const __m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1),pv=_mm256_set1_epi32((int)P);
const u32*W = dir? IW : RW; const u32*WS = dir? IWS : RWS;
if(!dir){
row8_dif<1>(row,128); row8_dif<0>(row,16); // stages 512..128 and 64..16
{ const u32*LT=W+(N2-8); const u32*LTS=WS+(N2-8);
V w4f,ws4f,w2f,ws2f,w2dd,ws2dd; tconst_p(LT,LTS,w4f,ws4f,w2f,ws2f,w2dd,ws2dd);
V w8,ws8; tailconst8(W+(N2-16),WS+(N2-16),w8,ws8);
for(u32 i=0;i<N2;i+=64) tail8_dif4p(row+i,w8,ws8,w4f,ws4f,w2f,ws2f,p4,p4m1,pv); }
} else {
{ const u32*LT=W+(N2-8); const u32*LTS=WS+(N2-8);
V w4f,ws4f,w2f,ws2f,w2dd,ws2dd; tconst_p(LT,LTS,w4f,ws4f,w2f,ws2f,w2dd,ws2dd);
V w8,ws8; tailconst8(W+(N2-16),WS+(N2-16),w8,ws8);
V pinv=_mm256_set1_epi32((int)g_pinv);
for(u32 i=0;i<N2;i+=64) tail8_dit4pp(row+i,(dir?qb:row)+i,pinv,w8,ws8,w4f,ws4f,w2dd,ws2dd,p4,p4m1,pv); }
row8_dit(row,16); row8_dit(row,128); // stages 16..64 and 128..512 (DIT ascends)
}
}
TGT static void row_dif_pw_z(u32*row,u32 dir){ row_dif_pw(row,row,dir); }
// ---------------- fused radix-8 column passes ----------------
// Fuses the three consecutive DIF (or DIT) stages with half-sizes 4h,2h,h (block = 8h) into
// ONE pass over the data. The twiddles come from the SAME per-stage tables the generic loop
// uses (offsets N1-8h / N1-4h / N1-2h), and the butterflies are applied in the same order, so
// the arithmetic and the output permutation are bit-for-bit those of the three separate
// passes -- no reindexing anywhere. Lane m of group (s,j) is the position s+j+m*h, i.e. the
// eight values are one stride-h row apart, so the pass touches 8 rows per 8 columns instead
// of 2 rows per 8 columns of the radix-2 loop: 5 passes instead of 11 over the same 176 MB.
template<u32 ST,u32 PFK=128,int D4=0> TGT static void col8_dif(u32*a,u32 h,u32 nrows=N1){
COL8SETUP;
const u32*W4=CW+(N1-8*h),*W4s=CWS+(N1-8*h);
const u32*W2=CW+(N1-4*h),*W2s=CWS+(N1-4*h);
const u32*W1=CW+(N1-2*h),*W1s=CWS+(N1-2*h);
for(u32 s=0;s<nrows;s+=8*h){
u32*rb=a+(size_t)s*ST;
for(u32 j=0;j<h;j++){
u32*r=rb+(size_t)j*ST;
__m256i t1v[4],t1sv[4],t20,t20s,t21,t21s,t3v,t3sv;
for(int m=0;m<4;m++){ t1v[m]=_mm256_set1_epi32((int)W4[j+m*h]); t1sv[m]=_mm256_set1_epi32((int)W4s[j+m*h]); }
t20=_mm256_set1_epi32((int)W2[j]); t20s=_mm256_set1_epi32((int)W2s[j]);
t21=_mm256_set1_epi32((int)W2[j+h]); t21s=_mm256_set1_epi32((int)W2s[j+h]);
t3v=_mm256_set1_epi32((int)W1[j]); t3sv=_mm256_set1_epi32((int)W1s[j]);
for(u32 c=0;c<N2;c+=8){
__m256i x[8];
for(int m=0;m<8;m++) _mm_prefetch((const char*)(r+(size_t)m*h*ST+c+PFK),_MM_HINT_T0);
for(int m=0;m<8;m++) x[m]=_mm256_load_si256((const __m256i*)(r+(size_t)m*h*ST+c));
for(int m=0;m<4;m++){ // stage 4h: (m, m+4)
// D4: both operands are in the row's <2P class (the (m,m+4) pair sits
// inside one 8h-block) => x+y < 4P => the 4P min is the identity.
__m256i A=D4?_mm256_add_epi32(x[m],x[m+4]):vadd4(x[m],x[m+4],p4,p4m1);
__m256i B=vshoup(vdiff4(x[m],x[m+4],p4),t1v[m],t1sv[m],pv);
__asm__("":"+x"(A),"+x"(B));
x[m]=A; x[m+4]=B;
}
for(int q=0;q<4;q++){ // stage 2h: (0,2)(1,3)(4,6)(5,7)
int lo=4*(q>>1)+(q&1), hi=lo+2;
__m256i tw=(lo&1)?t21:t20, tws=(lo&1)?t21s:t20s;
__m256i A=(q<2)?vadd4(x[lo],x[hi],p4,p4m1):_mm256_add_epi32(x[lo],x[hi]);
__m256i B=vshoup(vdiff4(x[lo],x[hi],p4),tw,tws,pv);
__asm__("":"+x"(A),"+x"(B));
x[lo]=A; x[hi]=B;
}
for(int m=0;m<8;m+=2){ // stage h: (m, m+1)
__m256i A=((m&3)==2)?_mm256_add_epi32(x[m],x[m+1]):vadd4(x[m],x[m+1],p4,p4m1);
__m256i B=vshoup(vdiff4(x[m],x[m+1],p4),t3v,t3sv,pv);
__asm__("":"+x"(A),"+x"(B));
x[m]=A; x[m+1]=B;
}
for(int m=0;m<8;m++) _mm256_store_si256((__m256i*)(r+(size_t)m*h*ST+c),x[m]);
}
}
}
}
template<u32 ST,int PS> TGT static void col8_dif1s(u32*dst,const u32*src,const u32*scr,u32 nrows,u32 h){
COL8SETUP;
const u32*W2=CW+(N1-4*h),*W2s=CWS+(N1-4*h);
const u32*W1=CW+(N1-2*h),*W1s=CWS+(N1-2*h);
u32 Rm=(u32)((1ull<<32)%P);
const __m256i Rv=_mm256_set1_epi32((int)Rm), Rs=_mm256_set1_epi32((int)(((u64)Rm<<32)/P));
for(u32 s=0;s<N1;s+=8*h){
for(u32 j=0;j<h;j++){
__m256i t1v[4],t1sv[4],t20,t20s,t21,t21s,t3v,t3sv;
for(int m=0;m<4;m++){ t1v[m]=_mm256_set1_epi32((int)CW[N1-8*h+j+m*h]); t1sv[m]=_mm256_set1_epi32((int)CWS[N1-8*h+j+m*h]); }
t20=_mm256_set1_epi32((int)W2[j]); t20s=_mm256_set1_epi32((int)W2s[j]);
t21=_mm256_set1_epi32((int)W2[j+h]); t21s=_mm256_set1_epi32((int)W2s[j+h]);
t3v=_mm256_set1_epi32((int)W1[j]); t3sv=_mm256_set1_epi32((int)W1s[j]);
for(u32 c=0;c<N2;c+=8){
__m256i x[8];
for(int m=0;m<4;m++) _mm_prefetch((const char*)(src+(size_t)(s+j+m*h)*N2+c+26),_MM_HINT_T0);
for(int m=0;m<4;m++){
u32 r=s+j+m*h; __m256i v;
if(r+1<nrows) v=LD256(src+(size_t)r*N2+c);
else if(r+1==nrows) v=LD256(scr+c);
else v=_mm256_setzero_si256();
if(PS) v=vshoup(v,Rv,Rs,pv);
if(PS) __asm__("":"+x"(v));
x[m]=v;
}
for(int m=0;m<4;m++) x[m+4]=vshoup(x[m],t1v[m],t1sv[m],pv);
// PS=1: x[0..3] are vshoup(raw) < 2P (Shoup) => stage 2h's q=0,1 pairs are
// B+B < 4P, where the 4P min is the identity. PS=0 keeps the raw input,
// which is only known < 4P, so those stay live.
for(int q=0;q<4;q++){
int lo=4*(q>>1)+(q&1), hi=lo+2;
__m256i tw=(lo&1)?t21:t20, tws=(lo&1)?t21s:t20s;
__m256i A=_mm256_add_epi32(x[lo],x[hi]);
__m256i B=vshoup(vdiff4(x[lo],x[hi],p4),tw,tws,pv);
__asm__("":"+x"(A),"+x"(B));
x[lo]=A; x[hi]=B;
}
for(int m=0;m<8;m+=2){
__m256i A=((m&3)==2)?_mm256_add_epi32(x[m],x[m+1]):vadd4(x[m],x[m+1],p4,p4m1);
__m256i B=vshoup(vdiff4(x[m],x[m+1],p4),t3v,t3sv,pv);
__asm__("":"+x"(A),"+x"(B));
x[m]=A; x[m+1]=B;
}
for(int m=0;m<8;m++){
u32 row=s+j+m*h;
u32* dp=(PS && g_bhi && row>=BSPLIT) ? (g_bhi+(size_t)(row-BSPLIT)*ST)
: (dst +(size_t)row*ST);
_mm256_stream_si256((__m256i*)(dp+c),x[m]); }
}
}
}
}
template<u32 ST> TGT static void col8_dif1(u32*a,u32 h){
COL8SETUP;
const u32*W4=CW+(N1-8*h),*W4s=CWS+(N1-8*h);
const u32*W2=CW+(N1-4*h),*W2s=CWS+(N1-4*h);
const u32*W1=CW+(N1-2*h),*W1s=CWS+(N1-2*h);
for(u32 s=0;s<N1;s+=8*h){
u32*rb=a+(size_t)s*ST;
for(u32 j=0;j<h;j++){
u32*r=rb+(size_t)j*ST;
__m256i t1v[4],t1sv[4],t20,t20s,t21,t21s,t3v,t3sv;
for(int m=0;m<4;m++){ t1v[m]=_mm256_set1_epi32((int)W4[j+m*h]); t1sv[m]=_mm256_set1_epi32((int)W4s[j+m*h]); }
t20=_mm256_set1_epi32((int)W2[j]); t20s=_mm256_set1_epi32((int)W2s[j]);
t21=_mm256_set1_epi32((int)W2[j+h]); t21s=_mm256_set1_epi32((int)W2s[j+h]);
t3v=_mm256_set1_epi32((int)W1[j]); t3sv=_mm256_set1_epi32((int)W1s[j]);
for(u32 c=0;c<N2;c+=8){
__m256i x[8];
for(int m=0;m<8;m++) x[m]=ldu256((r+(size_t)m*h*ST+c));
for(int m=0;m<4;m++){ // stage 4h: upper operand is 0
x[m+4]=vshoup(x[m],t1v[m],t1sv[m],pv);
}
for(int q=0;q<4;q++){ // stage 2h: (0,2)(1,3)(4,6)(5,7)
int lo=4*(q>>1)+(q&1), hi=lo+2;
__m256i tw=(lo&1)?t21:t20, tws=(lo&1)?t21s:t20s;
__m256i A=(q<2)?vadd4(x[lo],x[hi],p4,p4m1):_mm256_add_epi32(x[lo],x[hi]);
__m256i B=vshoup(vdiff4(x[lo],x[hi],p4),tw,tws,pv);
x[lo]=A; x[hi]=B;
}
for(int m=0;m<8;m+=2){ // stage h: (m, m+1)
__m256i A=((m&3)==2)?_mm256_add_epi32(x[m],x[m+1]):vadd4(x[m],x[m+1],p4,p4m1);
__m256i B=vshoup(vdiff4(x[m],x[m+1],p4),t3v,t3sv,pv);
x[m]=A; x[m+1]=B;
}
for(int m=0;m<8;m++) stu256((r+(size_t)m*h*ST+c), x[m]);
}
}
}
}
template<u32 ST,u32 PFK=256> TGT static void col8_dit(u32*a,u32 h,u32 nrows=N1){
const __m256i p2=_mm256_set1_epi32((int)P2),p8=_mm256_set1_epi32((int)P8);
COL8SETUP;
const u32*W4=JW+(N1-8*h),*W4s=JWS+(N1-8*h);
const u32*W2=JW+(N1-4*h),*W2s=JWS+(N1-4*h);
const u32*W1=JW+(N1-2*h),*W1s=JWS+(N1-2*h);
for(u32 s=0;s<nrows;s+=8*h){
u32*rb=a+(size_t)s*ST;
for(u32 j=0;j<h;j++){
u32*r=rb+(size_t)j*ST;
__m256i t1v[4],t1sv[4],t20,t20s,t21,t21s,t3v,t3sv;
for(int m=0;m<4;m++){ t1v[m]=_mm256_set1_epi32((int)W4[j+m*h]); t1sv[m]=_mm256_set1_epi32((int)W4s[j+m*h]); }
t20=_mm256_set1_epi32((int)W2[j]); t20s=_mm256_set1_epi32((int)W2s[j]);
t21=_mm256_set1_epi32((int)W2[j+h]); t21s=_mm256_set1_epi32((int)W2s[j+h]);
t3v=_mm256_set1_epi32((int)W1[j]); t3sv=_mm256_set1_epi32((int)W1s[j]);
for(u32 c=0;c<N2;c+=16){
__m256i x[8],y[8];
for(int m=0;m<8;m++) _mm_prefetch((const char*)(r+(size_t)m*h*ST+c+PFK),_MM_HINT_T0);
for(int m=0;m<8;m++){ x[m]=_mm256_load_si256((const __m256i*)(r+(size_t)m*h*ST+c));
y[m]=_mm256_load_si256((const __m256i*)(r+(size_t)m*h*ST+c+8)); }
for(int m=0;m<8;m+=2){ // stage h first (DIT order)
{ __m256i t=vshoup(x[m+1],t3v,t3sv,pv); __m256i u=x[m];
__asm__("":"+x"(t),"+x"(u)); x[m]=_mm256_add_epi32(u,t); x[m+1]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2); }
{ __m256i t=vshoup(y[m+1],t3v,t3sv,pv); __m256i u=y[m];
__asm__("":"+x"(t),"+x"(u)); y[m]=_mm256_add_epi32(u,t); y[m+1]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2); }
}
for(int q=0;q<4;q++){ // stage 2h
int lo=4*(q>>1)+(q&1), hi=lo+2;
__m256i tw=(lo&1)?t21:t20, tws=(lo&1)?t21s:t20s;
{ __m256i t=vshoup(x[hi],tw,tws,pv); __m256i u=x[lo];
__asm__("":"+x"(t),"+x"(u)); x[lo]=_mm256_add_epi32(u,t); x[hi]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2); }
{ __m256i t=vshoup(y[hi],tw,tws,pv); __m256i u=y[lo];
__asm__("":"+x"(t),"+x"(u)); y[lo]=_mm256_add_epi32(u,t); y[hi]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2); }
}
for(int m=0;m<4;m++){ // stage 4h last
{ __m256i t=vshoup(x[m+4],t1v[m],t1sv[m],pv); __m256i u=x[m];
__asm__("":"+x"(t),"+x"(u)); x[m]=_mm256_add_epi32(u,t); x[m+4]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2); }
{ __m256i t=vshoup(y[m+4],t1v[m],t1sv[m],pv); __m256i u=y[m];
__asm__("":"+x"(t),"+x"(u)); y[m]=_mm256_add_epi32(u,t); y[m+4]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2); }
}
for(int m=0;m<8;m++){ _mm256_store_si256((__m256i*)(r+(size_t)m*h*ST+c),x[m]); }
for(int m=0;m<8;m++){ _mm256_store_si256((__m256i*)(r+(size_t)m*h*ST+c+8),y[m]); }
}
}
}
}
// ---- the EMIT. The judge's output buffer c is 8 mod 32, so the old guarded
// _mm256_stream_si256 never fired and the storeu path was lowered by gcc-9 into
// vmovups %xmm + vextracti128 on every one of the 250k output stores. Here the
// 32-byte alignment is CONSTRUCTED: the column loop is shifted by `shift` so that
// (dst + rr*N2 + c) is 32-byte aligned for every row (N2=1024 is a multiple of 8,
// so the shift is the same for every row), the head/tail columns are covered by
// two extra storeu iterations, and the bulk becomes one vmovntdq per 8 elements.
// Unaligned vmovdqu loads read the (32-byte aligned) source rows at the shifted c.
// The rows rr>full are no longer written back into 'a': nothing ever read them
// (inv_pad's tail loop breaks at rr=full+1) and that write-back aliased the rows
// this split-pass loop still had to read. Only row 'full' is needed, into g_tail.
TGT static inline u32* OPQ_P(u32* p){ __asm__ __volatile__("" : "+r"(p)); return p; }
template<u32 ST> TGT static void col8_dit_dst(u32*a,u32 h,u32*dst,u32 full){
const __m256i p2=_mm256_set1_epi32((int)P2),p8=_mm256_set1_epi32((int)P8);
COL8SETUP;
const u32*W4=JW+(N1-8*h),*W4s=JWS+(N1-8*h);
const u32*W2=JW+(N1-4*h),*W2s=JWS+(N1-4*h);
const u32*W1=JW+(N1-2*h),*W1s=JWS+(N1-2*h);
const __m256i p2v=_mm256_set1_epi32((int)P2),p4v=_mm256_set1_epi32((int)P4);
const __m256i p8v=_mm256_set1_epi32((int)P8);
const __m256i p16v=_mm256_set1_epi32((int)(16u*P));
const u32 shift = (u32)((8u - (u32)(((size_t)dst >> 2) & 7u)) & 7u);
const u32 npass = shift ? 3u : 1u;
const u32 cend0 = shift + (((N2-shift)>>3)<<3);
for(u32 s=0;s<N1;s+=8*h){
u32*rb=a+(size_t)s*ST;
for(u32 j=0;j<h;j++){
u32*r=rb+(size_t)j*ST;
__m256i t1v[4],t1sv[4],t20,t20s,t21,t21s,t3v,t3sv;
for(int m=0;m<4;m++){ t1v[m]=_mm256_set1_epi32((int)W4[j+m*h]); t1sv[m]=_mm256_set1_epi32((int)W4s[j+m*h]); }
t20=_mm256_set1_epi32((int)W2[j]); t20s=_mm256_set1_epi32((int)W2s[j]);
t21=_mm256_set1_epi32((int)W2[j+h]); t21s=_mm256_set1_epi32((int)W2s[j+h]);
t3v=_mm256_set1_epi32((int)W1[j]); t3sv=_mm256_set1_epi32((int)W1s[j]);
for(u32 pass=0; pass<npass; ++pass){
u32 cstart = (pass==0)? shift : (pass==1? 0u : (N2-8u));
u32 cend = (pass==0)? cend0 : (cstart+8u);
const int AL = (pass==0);
for(u32 c=cstart;c<cend;c+=8){
u32* dbase=OPQ_P(dst+(size_t)(s+j)*N2);
__m256i x[8];
for(int m=0;m<8;m++) _mm_prefetch((const char*)(r+(size_t)m*h*ST+c+20),_MM_HINT_T0);
for(int m=0;m<8;m++) x[m]=LD256(r+(size_t)m*h*ST+c);
for(int m=0;m<8;m+=2){ // stage h first (DIT order)
__m256i t=vshoup(x[m+1],t3v,t3sv,pv);
__asm__("":"+x"(t));
__m256i u=x[m];
x[m]=_mm256_add_epi32(u,t);
x[m+1]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2);
}
for(int q=0;q<4;q++){ // stage 2h
int lo=4*(q>>1)+(q&1), hi=lo+2;
__m256i tw=(lo&1)?t21:t20, tws=(lo&1)?t21s:t20s;
__m256i t=vshoup(x[hi],tw,tws,pv);
__asm__("":"+x"(t));
__m256i u=x[lo];
x[lo]=_mm256_add_epi32(u,t);
x[hi]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2);
}
for(int m=0;m<4;m++){ // stage 4h last
__m256i t=vshoup(x[m+4],t1v[m],t1sv[m],pv);
__asm__("":"+x"(t));
__m256i u=x[m];
x[m]=_mm256_add_epi32(u,t);
x[m+4]=_mm256_add_epi32(_mm256_sub_epi32(u,t),p2);
}
for(int m=0;m<8;m++){
u32 rr=s+j+m*h;
if(rr<full){
__m256i v=x[m];
{ __m256i q=_mm256_srli_epi32(v,27);
v=_mm256_sub_epi32(v,_mm256_mullo_epi32(q,pv)); }
v=_mm256_min_epu32(v,_mm256_sub_epi32(v,pv));
if(AL) _mm256_stream_si256((__m256i*)(dbase+c+(size_t)m*h*N2),v);
else ST256(dbase+c+(size_t)m*h*N2,v);
} else if(rr==full){ __m256i xc=_mm256_min_epu32(x[m],_mm256_sub_epi32(x[m],p16v));
xc=_mm256_min_epu32(xc,_mm256_sub_epi32(xc,p8v)); ST256(g_tail+c,xc); }
}
}
}
}
}
}
// fused radix-4 pass covering the two DIF stages 2h and h (block 4h). The second level's
// twiddles are w_{2h}^j, which for the tail call h=1 are all 1 -> that level is pure add/sub.
template<u32 ST,int ZA=0> TGT static void col4_dif(u32*a,u32 h,u32 nrows=N1){
COL8SETUP;
const u32*W2=CW+(N1-4*h),*W2s=CWS+(N1-4*h);
const u32*W1=CW+(N1-2*h),*W1s=CWS+(N1-2*h);
u32 t20=W2[0],t20s=W2s[0],t21=W2[h],t21s=W2s[h],t3=W1[0],t3s=W1s[0];
__m256i v20=_mm256_set1_epi32((int)t20),v20s=_mm256_set1_epi32((int)t20s);
__m256i v21=_mm256_set1_epi32((int)t21),v21s=_mm256_set1_epi32((int)t21s);
__m256i v3=_mm256_set1_epi32((int)t3),v3s=_mm256_set1_epi32((int)t3s);
for(u32 s=0;s<nrows;s+=4*h){
u32*rb=a+(size_t)s*ST;
for(u32 j=0;j<h;j++){
u32*r=rb+(size_t)j*ST;
for(u32 c=0;c<N2;c+=16){
__m256i x0=_mm256_load_si256((const __m256i*)(r+c));
__m256i x1=_mm256_load_si256((const __m256i*)(r+(size_t)h*ST+c));
__m256i x2=_mm256_load_si256((const __m256i*)(r+(size_t)2*h*ST+c));
__m256i x3=_mm256_load_si256((const __m256i*)(r+(size_t)3*h*ST+c));
__m256i y0=_mm256_load_si256((const __m256i*)(r+c+8));
__m256i y1=_mm256_load_si256((const __m256i*)(r+(size_t)h*ST+c+8));
__m256i y2=_mm256_load_si256((const __m256i*)(r+(size_t)2*h*ST+c+8));
__m256i y3=_mm256_load_si256((const __m256i*)(r+(size_t)3*h*ST+c+8));
// ZA: all four operands are in the <2P class => x+y < 4P => the 4P min is
// the identity (Shoup). Output bound is <4P either way, stage B unmoved.
__m256i A0=ZA?_mm256_add_epi32(x0,x2):vadd4(x0,x2,p4,p4m1), B0=vsub4x(x0,x2,p4);
__m256i A1=ZA?_mm256_add_epi32(x1,x3):vadd4(x1,x3,p4,p4m1), B1=vshoup(vdiff4(x1,x3,p4),v21,v21s,pv);
__m256i C0=ZA?_mm256_add_epi32(y0,y2):vadd4(y0,y2,p4,p4m1), D0=vsub4x(y0,y2,p4);
__m256i C1=ZA?_mm256_add_epi32(y1,y3):vadd4(y1,y3,p4,p4m1), D1=vshoup(vdiff4(y1,y3,p4),v21,v21s,pv);
x0=A0; x2=B0; x1=A1; x3=B1;
y0=C0; y2=D0; y1=C1; y3=D1;
A0=vadd4(x0,x1,p4,p4m1); B0=vsub4x(x0,x1,p4);
A1=vadd4(x2,x3,p4,p4m1); B1=vsub4x(x2,x3,p4);
C0=vadd4(y0,y1,p4,p4m1); D0=vsub4x(y0,y1,p4);
C1=vadd4(y2,y3,p4,p4m1); D1=vsub4x(y2,y3,p4);
_mm256_store_si256((__m256i*)(r+c),A0);
_mm256_store_si256((__m256i*)(r+(size_t)h*ST+c),B0);
_mm256_store_si256((__m256i*)(r+(size_t)2*h*ST+c),A1);
_mm256_store_si256((__m256i*)(r+(size_t)3*h*ST+c),B1);
_mm256_store_si256((__m256i*)(r+c+8),C0);
_mm256_store_si256((__m256i*)(r+(size_t)h*ST+c+8),D0);
_mm256_store_si256((__m256i*)(r+(size_t)2*h*ST+c+8),C1);
_mm256_store_si256((__m256i*)(r+(size_t)3*h*ST+c+8),D1);
}
}
}
}
template<u32 ST> TGT static void col4_dit(u32*a,u32 h,u32 nrows=N1){
COL8SETUP;
const u32*W2=JW+(N1-4*h),*W2s=JWS+(N1-4*h);
__m256i v20=_mm256_set1_epi32((int)W2[0]),v20s=_mm256_set1_epi32((int)W2s[0]);
__m256i v21=_mm256_set1_epi32((int)W2[h]),v21s=_mm256_set1_epi32((int)W2s[h]);
for(u32 s=0;s<nrows;s+=4*h){
u32*rb=a+(size_t)s*ST;
for(u32 j=0;j<h;j++){
u32*r=rb+(size_t)j*ST;
for(u32 c=0;c<N2;c+=16){
__m256i x0=_mm256_load_si256((const __m256i*)(r+c));
__m256i y0=_mm256_load_si256((const __m256i*)(r+(size_t)h*ST+c));
__m256i x1=_mm256_load_si256((const __m256i*)(r+(size_t)2*h*ST+c));
__m256i y1=_mm256_load_si256((const __m256i*)(r+(size_t)3*h*ST+c));
__m256i x2=_mm256_load_si256((const __m256i*)(r+c+8));
__m256i y2=_mm256_load_si256((const __m256i*)(r+(size_t)h*ST+c+8));
__m256i x3=_mm256_load_si256((const __m256i*)(r+(size_t)2*h*ST+c+8));
__m256i y3=_mm256_load_si256((const __m256i*)(r+(size_t)3*h*ST+c+8));
__m256i a0=_mm256_add_epi32(x0,y0); __m256i b0=_mm256_add_epi32(_mm256_sub_epi32(x0,y0),p4);
__m256i a1=_mm256_add_epi32(x1,y1); __m256i b1=_mm256_add_epi32(_mm256_sub_epi32(x1,y1),p4);
__m256i a2=_mm256_add_epi32(x2,y2); __m256i b2=_mm256_add_epi32(_mm256_sub_epi32(x2,y2),p4);
__m256i a3=_mm256_add_epi32(x3,y3); __m256i b3=_mm256_add_epi32(_mm256_sub_epi32(x3,y3),p4);
__m256i t2=a1;
__m256i o0=_mm256_add_epi32(a0,t2); __m256i o2=_mm256_add_epi32(_mm256_sub_epi32(a0,t2),p4);
__m256i t3=vshoup(b1,v21,v21s,pv);
__m256i o1=_mm256_add_epi32(b0,t3); __m256i o3=_mm256_add_epi32(_mm256_sub_epi32(b0,t3),p4);
__m256i u2=a3;
__m256i p0=_mm256_add_epi32(a2,u2); __m256i p2=_mm256_add_epi32(_mm256_sub_epi32(a2,u2),p4);
__m256i u3=vshoup(b3,v21,v21s,pv);
__m256i p1=_mm256_add_epi32(b2,u3); __m256i p3=_mm256_add_epi32(_mm256_sub_epi32(b2,u3),p4);
_mm256_store_si256((__m256i*)(r+c),o0);
_mm256_store_si256((__m256i*)(r+(size_t)h*ST+c),o1);
_mm256_store_si256((__m256i*)(r+(size_t)2*h*ST+c),o2);
_mm256_store_si256((__m256i*)(r+(size_t)3*h*ST+c),o3);
_mm256_store_si256((__m256i*)(r+c+8),p0);
_mm256_store_si256((__m256i*)(r+(size_t)h*ST+c+8),p1);
_mm256_store_si256((__m256i*)(r+(size_t)2*h*ST+c+8),p2);
_mm256_store_si256((__m256i*)(r+(size_t)3*h*ST+c+8),p3);
}
}
}
}
// ---------------- column transform: length N1 with stride N2, blocked by 8 columns ----------------
template<u32 ST> TGT static void col_tr(u32*a,u32 dir){
const __m256i p4=_mm256_set1_epi32((int)P4),p4m1=_mm256_set1_epi32((int)P4M1),pv=_mm256_set1_epi32((int)g_pvs);
const u32*W = dir? JW : CW; const u32*WS = dir? JWS : CWS;
if(!dir){
col8_dif1<ST>(a,256); col8_dif<ST>(a,32); col8_dif<ST>(a,4); // stages 1024..4
col4_dif<ST>(a,1); // stages 2 and 1
for(u32 h=0;h<0;h>>=1){
const u32*w=W+(N1-2*h),*ws=WS+(N1-2*h);
for(u32 s=0;s<N1;s+=2*h){ u32*r0=a+(size_t)s*ST,*r1=r0+(size_t)h*ST;
for(u32 j=0;j<h;j++,r0+=N2,r1+=N2){
__m256i wv=_mm256_set1_epi32((int)w[j]),wsv=_mm256_set1_epi32((int)ws[j]);
for(u32 c=0;c<N2;c+=32){
__m256i x0=ldu256((r0+c)),y0=ldu256((r1+c));
__m256i x1=ldu256((r0+c+8)),y1=ldu256((r1+c+8));
__m256i x2=ldu256((r0+c+16)),y2=ldu256((r1+c+16));
__m256i x3=ldu256((r0+c+24)),y3=ldu256((r1+c+24));
stu256((r0+c), vadd4(x0,y0,p4,p4m1));
stu256((r0+c+8), vadd4(x1,y1,p4,p4m1));
stu256((r0+c+16), vadd4(x2,y2,p4,p4m1));
stu256((r0+c+24), vadd4(x3,y3,p4,p4m1));
stu256((r1+c), vshoup(vdiff4(x0,y0,p4),wv,wsv,pv));
stu256((r1+c+8), vshoup(vdiff4(x1,y1,p4),wv,wsv,pv));
stu256((r1+c+16), vshoup(vdiff4(x2,y2,p4),wv,wsv,pv));
stu256((r1+c+24), vshoup(vdiff4(x3,y3,p4),wv,wsv,pv));
}
}}
}
} else {
for(u32 h=3;h<=2;h<<=1){
const u32*w=W+(N1-2*h),*ws=WS+(N1-2*h);
for(u32 s=0;s<N1;s+=2*h){ u32*r0=a+(size_t)s*ST,*r1=r0+(size_t)h*ST;
for(u32 j=0;j<h;j++,r0+=N2,r1+=N2){
__m256i wv=_mm256_set1_epi32((int)w[j]),wsv=_mm256_set1_epi32((int)ws[j]);
for(u32 c=0;c<N2;c+=32){
__m256i x0=ldu256((r0+c)),y0=ldu256((r1+c));
__m256i x1=ldu256((r0+c+8)),y1=ldu256((r1+c+8));
__m256i x2=ldu256((r0+c+16)),y2=ldu256((r1+c+16));
__m256i x3=ldu256((r0+c+24)),y3=ldu256((r1+c+24));
__m256i t0=vshoup(y0,wv,wsv,pv),t1=vshoup(y1,wv,wsv,pv),t2=vshoup(y2,wv,wsv,pv),t3=vshoup(y3,wv,wsv,pv);
__m256i s0=vadd4(x0,t0,p4,p4m1),s1=vadd4(x1,t1,p4,p4m1),s2=vadd4(x2,t2,p4,p4m1),s3=vadd4(x3,t3,p4,p4m1);
__m256i d0=vsubr(x0,t0,p4),d1=vsubr(x1,t1,p4),d2=vsubr(x2,t2,p4),d3=vsubr(x3,t3,p4);
stu256((r0+c), s0);stu256((r0+c+8), s1);
stu256((r0+c+16), s2);stu256((r0+c+24), s3);
stu256((r1+c), d0);stu256((r1+c+8), d1);
stu256((r1+c+16), d2);stu256((r1+c+24), d3);
}
}}
}
col4_dit<ST>(a,1); col8_dit<ST>(a,4); col8_dit<ST>(a,32); col8_dit<ST>(a,256);
}
}
// ---------------- diagonal (geometric sequence), Montgomery chain (from n4.h) ----------------
// ================= lane p1002_chirp_wsa: vmont -> vshoup =================
// ws(K) = floor(K*2^32/P) mod 2^32 for K in [0,2P) (including the *unreduced* Shoup outputs of
// the constant chain) is 32K + mulhi(K,Q64), exact or exactly ONE LOW -- never anything else.
// Q64 = floor((2^32 mod P)*2^32/P) = floor(Rm*2^32/P) = 2181569633
// Rm = 2^32 mod P = 2^26-32 ; 2^32 = 32P + Rm
// Verified exhaustively on [0,2P) (264241154 values, 0 exceptions, 0 high-side errors) and on
// all 2097152 constants the chirp actually builds (99.5821% exact, 0.4179% one-low, 0 other).
// A single mulhi can NEVER be exact (Granlund-Montgomery condition violated) -- see FINDINGS.
static const u32 Q64M = 2181569633u;
TGT static inline __m256i wsfv(__m256i K,__m256i q64v){
return _mm256_add_epi32(_mm256_slli_epi32(K,5),vmulhi32(K,q64v));
}
TGT static inline __m256i vredP(__m256i v,__m256i pv){ // v<2P -> [0,P)
return _mm256_min_epu32(v,_mm256_sub_epi32(v,pv));
}
// ===== MODE 2: 32-element step. SAME instruction count as MODE 1 (one chain step + one ws
// + one reduce per constant vector, 2 data multiplies per constant), but only 4 constant
// vectors are live instead of 8 and the d-vectors disappear entirely (the chain advances by
// r32 instead of r64), so the register pressure is halved and the spills go away.
template<int MODE> TGT static void row_scale2b(u32*ra,u32*rb,u32 start,u32 ratio){
(void)MODE;
const __m256i pv=_mm256_set1_epi32((int)g_pvs);
const __m256i q64v=_mm256_set1_epi32((int)Q64M);
u32 v[8];
v[0]=(u32)(start%P);
for(int t=1;t<8;t++) v[t]=(u32)((u64)v[t-1]*ratio%P);
u32 r2=(u32)((u64)ratio*ratio%P);
u32 r4=(u32)((u64)r2*r2%P);
u32 r8=(u32)((u64)r4*r4%P);
u32 r16=(u32)((u64)r8*r8%P);
u32 r32b=(u32)((u64)r16*r16%P);
__m256i c0=ldu256(v);
const __m256i r8v=_mm256_set1_epi32((int)r8), r8s=_mm256_set1_epi32((int)ws_of(r8));
__m256i c1=vshoup(c0,r8v,r8s,pv);
__m256i c2=vshoup(c1,r8v,r8s,pv);
__m256i c3=vshoup(c2,r8v,r8s,pv);
c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv);
const __m256i stpq=_mm256_set1_epi32((int)r32b), stpws=_mm256_set1_epi32((int)ws_of(r32b));
for(u32 j=0;j<N2;j+=32){
{ __m256i w=wsfv(c0,q64v);
stu256((ra+j), vshoup(ldu256((ra+j)), c0,w,pv));
stu256((rb+j), vshoup(ldu256((rb+j)), c0,w,pv)); }
{ __m256i w=wsfv(c1,q64v);
stu256((ra+j+8), vshoup(ldu256((ra+j+8)),c1,w,pv));
stu256((rb+j+8), vshoup(ldu256((rb+j+8)),c1,w,pv)); }
{ __m256i w=wsfv(c2,q64v);
stu256((ra+j+16), vshoup(ldu256((ra+j+16)),c2,w,pv));
stu256((rb+j+16), vshoup(ldu256((rb+j+16)),c2,w,pv)); }
{ __m256i w=wsfv(c3,q64v);
stu256((ra+j+24), vshoup(ldu256((ra+j+24)),c3,w,pv));
stu256((rb+j+24), vshoup(ldu256((rb+j+24)),c3,w,pv)); }
c0=vshoup(c0,stpq,stpws,pv); c1=vshoup(c1,stpq,stpws,pv);
c2=vshoup(c2,stpq,stpws,pv); c3=vshoup(c3,stpq,stpws,pv);
c0=vredP(c0,pv); c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv);
}
}
TGT static void row_scale32(u32*ra,u32*rb,u32 start,u32 ratio){ row_scale2b<2>(ra,rb,start,ratio); }
template<int MODE> TGT static void row_scale2m(u32*ra,u32*rb,u32 start,u32 ratio){
const __m256i pv=_mm256_set1_epi32((int)g_pvs);
__m256i pinv=_mm256_set1_epi32((int)g_pinv);
u32 Rm=(u32)((1ull<<32)%P);
const __m256i q64v=_mm256_set1_epi32((int)Q64M);
u32 v[8];
v[0]=MODE?(u32)(start%P):(u32)((u64)(start%P)*Rm%P);
for(int t=1;t<8;t++) v[t]=(u32)((u64)v[t-1]*ratio%P);
u32 r2=(u32)((u64)ratio*ratio%P);
u32 r4=(u32)((u64)r2*r2%P);
u32 r8=(u32)((u64)r4*r4%P);
u32 r16=(u32)((u64)r8*r8%P);
u32 r32=(u32)((u64)r16*r16%P);
u32 r64=(u32)((u64)r32*r32%P);
__m256i c0=ldu256(v);
const __m256i r8v=_mm256_set1_epi32((int)r8), r8s=_mm256_set1_epi32((int)ws_of(r8));
__m256i c1=vshoup(c0,r8v,r8s,pv);
__m256i c2=vshoup(c1,r8v,r8s,pv);
__m256i c3=vshoup(c2,r8v,r8s,pv);
if(MODE){ c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv); }
const __m256i stpq=_mm256_set1_epi32((int)r32), stpws=_mm256_set1_epi32((int)ws_of(r32));
const __m256i stpq2=_mm256_set1_epi32((int)r64), stpws2=_mm256_set1_epi32((int)ws_of(r64));
for(u32 j=0;j<N2;j+=64){
const __m256i q64v2=q64v;
__m256i d0=vshoup(c0,stpq,stpws,pv), d1=vshoup(c1,stpq,stpws,pv);
__m256i d2=vshoup(c2,stpq,stpws,pv), d3=vshoup(c3,stpq,stpws,pv);
if(MODE){ d0=vredP(d0,pv); d1=vredP(d1,pv); d2=vredP(d2,pv); d3=vredP(d3,pv); }
{ const __m256i*R=(const __m256i*)(ra+j); __m256i x0=ldu256(R+0),x1=ldu256(R+1),
x2=ldu256(R+2),x3=ldu256(R+3),x4=ldu256(R+4),
x5=ldu256(R+5),x6=ldu256(R+6),x7=ldu256(R+7);
if(MODE){
__m256i w0=wsfv(c0,q64v2),w1=wsfv(c1,q64v2),w2=wsfv(c2,q64v2),w3=wsfv(c3,q64v2),
w4=wsfv(d0,q64v2),w5=wsfv(d1,q64v2),w6=wsfv(d2,q64v2),w7=wsfv(d3,q64v2);
stu256((ra+j), vshoup(x0,c0,w0,pv));
stu256((ra+j+8), vshoup(x1,c1,w1,pv));
stu256((ra+j+16), vshoup(x2,c2,w2,pv));
stu256((ra+j+24), vshoup(x3,c3,w3,pv));
stu256((ra+j+32), vshoup(x4,d0,w4,pv));
stu256((ra+j+40), vshoup(x5,d1,w5,pv));
stu256((ra+j+48), vshoup(x6,d2,w6,pv));
stu256((ra+j+56), vshoup(x7,d3,w7,pv));
} else {
stu256((ra+j), vmont(x0,c0,pv,pinv));
stu256((ra+j+8), vmont(x1,c1,pv,pinv));
stu256((ra+j+16), vmont(x2,c2,pv,pinv));
stu256((ra+j+24), vmont(x3,c3,pv,pinv));
stu256((ra+j+32), vmont(x4,d0,pv,pinv));
stu256((ra+j+40), vmont(x5,d1,pv,pinv));
stu256((ra+j+48), vmont(x6,d2,pv,pinv));
stu256((ra+j+56), vmont(x7,d3,pv,pinv));
} }
{ const __m256i*R=(const __m256i*)(rb+j); __m256i x0=ldu256(R+0),x1=ldu256(R+1),
x2=ldu256(R+2),x3=ldu256(R+3),x4=ldu256(R+4),
x5=ldu256(R+5),x6=ldu256(R+6),x7=ldu256(R+7);
if(MODE){
__m256i w0=wsfv(c0,q64v2),w1=wsfv(c1,q64v2),w2=wsfv(c2,q64v2),w3=wsfv(c3,q64v2),
w4=wsfv(d0,q64v2),w5=wsfv(d1,q64v2),w6=wsfv(d2,q64v2),w7=wsfv(d3,q64v2);
stu256((rb+j), vshoup(x0,c0,w0,pv));
stu256((rb+j+8), vshoup(x1,c1,w1,pv));
stu256((rb+j+16), vshoup(x2,c2,w2,pv));
stu256((rb+j+24), vshoup(x3,c3,w3,pv));
stu256((rb+j+32), vshoup(x4,d0,w4,pv));
stu256((rb+j+40), vshoup(x5,d1,w5,pv));
stu256((rb+j+48), vshoup(x6,d2,w6,pv));
stu256((rb+j+56), vshoup(x7,d3,w7,pv));
} else {
stu256((rb+j), vmont(x0,c0,pv,pinv));
stu256((rb+j+8), vmont(x1,c1,pv,pinv));
stu256((rb+j+16), vmont(x2,c2,pv,pinv));
stu256((rb+j+24), vmont(x3,c3,pv,pinv));
stu256((rb+j+32), vmont(x4,d0,pv,pinv));
stu256((rb+j+40), vmont(x5,d1,pv,pinv));
stu256((rb+j+48), vmont(x6,d2,pv,pinv));
stu256((rb+j+56), vmont(x7,d3,pv,pinv));
} }
c0=vshoup(c0,stpq2,stpws2,pv); c1=vshoup(c1,stpq2,stpws2,pv);
c2=vshoup(c2,stpq2,stpws2,pv); c3=vshoup(c3,stpq2,stpws2,pv);
if(MODE){ c0=vredP(c0,pv); c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv); }
}
}
template<int MODE> TGT static void row_scalem(u32*row,u32 start,u32 ratio){
const __m256i pv=_mm256_set1_epi32((int)g_pvs);
__m256i pinv=_mm256_set1_epi32((int)g_pinv);
u32 Rm=(u32)((1ull<<32)%P);
const __m256i q64v=_mm256_set1_epi32((int)Q64M);
u32 v[8] __attribute__((aligned(32)));
v[0]=MODE?(u32)(start%P):(u32)((u64)(start%P)*Rm%P);
for(int t=1;t<8;t++) v[t]=(u32)((u64)v[t-1]*ratio%P);
u32 r2=(u32)((u64)ratio*ratio%P);
u32 r4=(u32)((u64)r2*r2%P);
u32 r8=(u32)((u64)r4*r4%P);
u32 r16=(u32)((u64)r8*r8%P);
u32 r32=(u32)((u64)r16*r16%P);
u32 r64=(u32)((u64)r32*r32%P);
__m256i c0=_mm256_load_si256((const __m256i*)v);
const __m256i r8v=_mm256_set1_epi32((int)r8), r8s=_mm256_set1_epi32((int)ws_of(r8));
__m256i c1=vshoup(c0,r8v,r8s,pv);
__m256i c2=vshoup(c1,r8v,r8s,pv);
__m256i c3=vshoup(c2,r8v,r8s,pv);
if(MODE){ c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv); }
const __m256i stpq=_mm256_set1_epi32((int)r32), stpws=_mm256_set1_epi32((int)ws_of(r32));
const __m256i stpq2=_mm256_set1_epi32((int)r64), stpws2=_mm256_set1_epi32((int)ws_of(r64));
if(MODE){
for(u32 j=0;j<N2;j+=32){
{ __m256i w0=wsfv(c0,q64v),w1=wsfv(c1,q64v);
_mm256_store_si256((__m256i*)(row+j), vshoup(_mm256_load_si256((const __m256i*)(row+j)), c0,w0,pv));
_mm256_store_si256((__m256i*)(row+j+8), vshoup(_mm256_load_si256((const __m256i*)(row+j+8)), c1,w1,pv));}
{ __m256i w2=wsfv(c2,q64v),w3=wsfv(c3,q64v);
_mm256_store_si256((__m256i*)(row+j+16),vshoup(_mm256_load_si256((const __m256i*)(row+j+16)),c2,w2,pv));
_mm256_store_si256((__m256i*)(row+j+24),vshoup(_mm256_load_si256((const __m256i*)(row+j+24)),c3,w3,pv));}
c0=vshoup(c0,stpq,stpws,pv); c1=vshoup(c1,stpq,stpws,pv);
c2=vshoup(c2,stpq,stpws,pv); c3=vshoup(c3,stpq,stpws,pv);
c0=vredP(c0,pv); c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv);
}
} else {
for(u32 j=0;j<N2;j+=64){
__m256i x0=_mm256_load_si256((const __m256i*)(row+j));
__m256i x1=_mm256_load_si256((const __m256i*)(row+j+8));
__m256i x2=_mm256_load_si256((const __m256i*)(row+j+16));
__m256i x3=_mm256_load_si256((const __m256i*)(row+j+24));
__m256i d0=vshoup(c0,stpq,stpws,pv), d1=vshoup(c1,stpq,stpws,pv);
__m256i d2=vshoup(c2,stpq,stpws,pv), d3=vshoup(c3,stpq,stpws,pv);
if(MODE){ d0=vredP(d0,pv); d1=vredP(d1,pv); d2=vredP(d2,pv); d3=vredP(d3,pv); }
__m256i x4=_mm256_load_si256((const __m256i*)(row+j+32));
__m256i x5=_mm256_load_si256((const __m256i*)(row+j+40));
__m256i x6=_mm256_load_si256((const __m256i*)(row+j+48));
__m256i x7=_mm256_load_si256((const __m256i*)(row+j+56));
if(MODE){
__m256i w0=wsfv(c0,q64v),w1=wsfv(c1,q64v),w2=wsfv(c2,q64v),w3=wsfv(c3,q64v),
w4=wsfv(d0,q64v),w5=wsfv(d1,q64v),w6=wsfv(d2,q64v),w7=wsfv(d3,q64v);
_mm256_store_si256((__m256i*)(row+j), vshoup(x0,c0,w0,pv));
_mm256_store_si256((__m256i*)(row+j+8), vshoup(x1,c1,w1,pv));
_mm256_store_si256((__m256i*)(row+j+16),vshoup(x2,c2,w2,pv));
_mm256_store_si256((__m256i*)(row+j+24),vshoup(x3,c3,w3,pv));
_mm256_store_si256((__m256i*)(row+j+32),vshoup(x4,d0,w4,pv));
_mm256_store_si256((__m256i*)(row+j+40),vshoup(x5,d1,w5,pv));
_mm256_store_si256((__m256i*)(row+j+48),vshoup(x6,d2,w6,pv));
_mm256_store_si256((__m256i*)(row+j+56),vshoup(x7,d3,w7,pv));
} else {
_mm256_store_si256((__m256i*)(row+j), vmont(x0,c0,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+8), vmont(x1,c1,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+16),vmont(x2,c2,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+24),vmont(x3,c3,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+32),vmont(x4,d0,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+40),vmont(x5,d1,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+48),vmont(x6,d2,pv,pinv));
_mm256_store_si256((__m256i*)(row+j+56),vmont(x7,d3,pv,pinv)); }
c0=vshoup(c0,stpq2,stpws2,pv); c1=vshoup(c1,stpq2,stpws2,pv);
c2=vshoup(c2,stpq2,stpws2,pv); c3=vshoup(c3,stpq2,stpws2,pv);
if(MODE){ c0=vredP(c0,pv); c1=vredP(c1,pv); c2=vredP(c2,pv); c3=vredP(c3,pv); }
}
}
}
// ---------------- function interface ----------------
static u32 g_pw[N1], g_pwi[N1];
TGT static void pointwise(u32*a,u32*b){
const __m256i pv=_mm256_set1_epi32((int)g_pvs);
__m256i pinv=_mm256_set1_epi32((int)g_pinv);
for(u32 i=0;i<N1*STDFT;i+=8){
__m256i x=_mm256_load_si256((const __m256i*)(a+i));
__m256i y=_mm256_load_si256((const __m256i*)(b+i));
_mm256_store_si256((__m256i*)(a+i),vmont(x,y,pv,pinv));
}
}
static int g_first=1;
static u32 g_scrA[N2], g_scrB[N2];
template<u32 ST,int PS> TGT static void fwd_pad_src(u32*dst,const u32*src,const u32*scr,u32 nrows){
col8_dif1s<ST,PS>(dst,src,scr,nrows,256);
col8_dif<ST>(dst,32); col8_dif<ST>(dst,4); col4_dif<ST>(dst,1);
}
template<u32 ST> TGT static void fwd_pad(u32*a,u32 srclen){
u32 rows=(srclen+N2-1)/N2;
if(!g_first) for(u32 i=rows;i<N1;i++) memset(a+(size_t)i*ST,0,(size_t)N2*4);
col_tr<ST>(a,0);
for(u32 i=0;i<N1;i++){ u32 k1=REV1[i]; row_scalem<0>(a+(size_t)i*ST,1,g_pw[k1]); row_dif_pw_z(a+(size_t)i*ST,0); }
}
template<u32 ST,int DOROW=0,int IBLK=32> TGT static void inv_pad(u32*dst,u32 dstlen,u32*a){
if(DOROW==1) for(u32 s=0;s<N1;s+=IBLK){
for(u32 i=s;i<s+IBLK;i++){ u32 k1=REV1[i]; row_dif_pw(a+(size_t)i*ST, gBp[i], 1); row_scalem<0>(a+(size_t)i*ST,g_ninv,g_pwi[k1]); }
for(u32 t=0;t<IBLK;t+=4) col4_dit<ST>(a+(size_t)(s+t)*ST,1,4);
col8_dit<ST>(a+(size_t)s*ST,4,IBLK); }
{ if(DOROW==0){ col4_dit<ST>(a,1); col8_dit<ST>(a,4); }
u32 full = dstlen/N2;
col8_dit_dst<ST>(a,256,dst,full);
const __m256i pv2=_mm256_set1_epi32((int)g_pvs),p22=_mm256_set1_epi32((int)P2),p42=_mm256_set1_epi32((int)P4);
for(u32 rr=full; rr<N1; rr++){
u32 off=rr*N2; if(off>=dstlen) break;
u32 lim=(dstlen-off>=N2)?N2:(dstlen-off); u32 i=0; u32* src=g_tail;
for(;i+8<=lim;i+=8){ __m256i v=LD256(src+i);
v=_mm256_min_epu32(v,_mm256_sub_epi32(v,p42));
v=_mm256_min_epu32(v,_mm256_sub_epi32(v,p22));
v=_mm256_min_epu32(v,_mm256_sub_epi32(v,pv2));
ST256(dst+off+i,v); }
for(;i<lim;i++){ u32 v=src[i]; if(v>=P4)v-=P4; if(v>=P2)v-=P2; if(v>=P)v-=P; dst[off+i]=v; }
}
}
(void)dstlen;
}
TGT static void prep(){
init_small();
{ u32 a=1,b=1; u32 wi=powmod32(g_w,P-2);
for(u32 i=0;i<N1;i++){ g_pw[i]=a; g_pwi[i]=b; a=(u32)((u64)a*g_w%P); b=(u32)((u64)b*wi%P); } }
}
TGT static u64 g_maxin=0,g_maxout=0,g_maxrs=0;
template<int MODE> TGT void polyQ(unsigned *a, int n, unsigned *b, int m, unsigned *c){
g_pvs = g_vsrc;
g_nt = (((size_t)c & 31u)==0);
prep();
u32 la=(u32)n+1, lb=(u32)m+1;
u32 rowsA=(la+N2-1)/N2, rowsB=(lb+N2-1)/N2;
{ u32 off=(rowsA-1)*N2, lm=la-off;
memcpy(g_scrA,a+off,(size_t)lm*4); memset(g_scrA+lm,0,(size_t)(N2-lm)*4); }
{ u32 off=(rowsB-1)*N2, lm=lb-off;
memcpy(g_scrB,b+off,(size_t)lm*4); memset(g_scrB+lm,0,(size_t)(N2-lm)*4); }
u32 tot=la+lb-1;
{ /* c may be only 4B aligned: step up to the next 32B boundary INSIDE c. */
u32* cb=(u32*)(((uintptr_t)c + 31u) & ~(uintptr_t)31u);
size_t skip=(size_t)(cb-(u32*)c);
int use = ((size_t)tot >= (size_t)BSPLIT*STDFT + skip);
g_bhi = use ? g_BX : 0;
g_blo = use ? cb : g_B;
if(use){ for(u32 i=0;i<BSPLIT;i++) gBp[i]=cb+(size_t)i*STDFT;
for(u32 i=BSPLIT;i<N1;i++) gBp[i]=g_BX+(size_t)(i-BSPLIT)*STDFT; }
else { for(u32 i=0;i<N1;i++) gBp[i]=g_B+(size_t)i*STDFT; } }
col8_dif1s<STDFT,0>(g_A,a,g_scrA,rowsA,256);
col8_dif1s<STDFT,1>(g_blo,b,g_scrB,rowsB,256);
// h=32: 8h=256-row blocks; the pad left rows (row%512)<256 at <4P and the other
// half at <2P, so stage 4h is the identity in blocks 1,3,5,7 (D4=1).
for(u32 b0=N1-256;b0<N1;b0-=256){
if((b0>>8)&1u){ col8_dif<STDFT,128,1>(g_A+(size_t)b0*STDFT,32,256);
col8_dif<STDFT,128,1>(gBp[b0],32,256); }
else { col8_dif<STDFT,128,0>(g_A+(size_t)b0*STDFT,32,256);
col8_dif<STDFT,128,0>(gBp[b0],32,256); }
for(u32 k16=0;k16<8;k16++){ u32 s=b0+224u-32u*k16;
if(s&32){ col8_dif<STDFT,24,1>(g_A+(size_t)s*STDFT,4,32);
col8_dif<STDFT,24,1>(gBp[s],4,32); }
else { col8_dif<STDFT,24,0>(g_A+(size_t)s*STDFT,4,32);
col8_dif<STDFT,24,0>(gBp[s],4,32); }
for(u32 t=0;t<32;t+=4){
if(t&4) col4_dif<STDFT,1>(g_A+(size_t)(s+t)*STDFT,1,4);
else col4_dif<STDFT,0>(g_A+(size_t)(s+t)*STDFT,1,4);
if(t&4) col4_dif<STDFT,1>(gBp[s+t],1,4);
else col4_dif<STDFT,0>(gBp[s+t],1,4);
for(u32 i=s+t;i<s+t+4;i++){ u32 k1=REV1[i];
if(TRK) for(u32 q=0;q<N2;q++){ u32 vv=g_A[(size_t)i*STDFT+q]; if(vv>g_maxin) g_maxin=vv; }
if(MODE==2) row_scale32(g_A+(size_t)i*STDFT,gBp[i],1,g_pw[k1]);
else row_scale2m<MODE>(g_A+(size_t)i*STDFT,gBp[i],1,g_pw[k1]);
if(TRK) for(u32 q=0;q<N2;q++){ u32 vv=g_A[(size_t)i*STDFT+q]; if(vv>g_maxrs) g_maxrs=vv; }
row_dif_pw_z(g_A+(size_t)i*STDFT,0); row_dif_pw_z(gBp[i],0);
row_dif_pw(g_A+(size_t)i*STDFT,gBp[i],1);
row_scalem<MODE>(g_A+(size_t)i*STDFT,g_ninv,g_pwi[k1]); }
col4_dit<STDFT>(g_A+(size_t)(s+t)*STDFT,1,4);
}
col8_dit<STDFT,24>(g_A+(size_t)s*STDFT,4,32);
}
col8_dit<STDFT>(g_A+(size_t)b0*STDFT,32,256);
}
inv_pad<STDFT,2,32>(c,tot,g_A);
g_first=0;
if(g_nt) _mm_sfence();
}
// A2D a4-0005-eeeeeeee
// ISSUE_mulldB 5903e824
// q5_nop2 site=tail8_dif4p h2_vpermd_deleted
// cx_fin_v3: v1 + col8_dit (3 stages) + col8_dit_dst (3 stages)
TGT void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c){
polyQ<2>(a,n,b,m,c);
}
// ---- lwc1002e8_lane/1002: B's transform array relocated into the OUTPUT buffer c.
// c must hold n+m+1 = 2000001 u32; rows 0..1791 of B never need more than 1792*STDFT
// = 1906688 elements, so the write stays strictly inside c. Rows 1792..2047 (a whole
// 256-row column block) stay in g_BX. The old g_B region is never touched, so its
// 2128 written pages disappear from mem_kb (measured 855 ticks per first-touched page
// on this row). Falls back to the standing layout when c is small or not 32B aligned.
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 18.803 ms | 17 MB + 104 KB | Accepted | Score: 100 | 显示更多 |