代码

// DUCK_PUBLIC_SURROGATE_BENCHMARK
#define BENCH_P 2013265921u
#define BENCH_GOOD 1
#include <bits/stdc++.h>
#include <x86intrin.h>
// Coprime prime-factor decomposition, I. J. Good (1958),
// https://academic.oup.com/jrsssb/article/20/2/361/7027226
// Algebra and wrapper code here are independent; fast_field notices retained.
// Derived from public Duck submission 141230, attributed notices retained below.
// Pure field extraction: no 0..9 digit wrapper, no production inputs or seeds.
// Generic padded outer stages and exact generic Shoup arithmetic are original.
#include <immintrin.h>
#include <cstdint>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <vector>
#include <algorithm>
#pragma GCC push_options
namespace fast_e7_field {
static uint32_t g_mod=2013265921u,g_generator=31u,g_shmagic=572662301u,g_shint=2u;
/* seat14_fresh_20261031 agg4 2026-10-31 0submit: base cand_twdel_agg4.cpp (sha256 8b8c7813, 101919cp, sid141141) + ONE change: dh4's fused S=4 level, four horizontal final_dif4_vec calls -> transposed vertical radix-4 (two 4x4 lane transposes, one vS, same mod-p math); inherited notices below unchanged. Full 80000004B output identical: mode0 FNV fb206c93723a3583 C0=10, all-9 77fffe38ccc593f7 C0=81, unaligned-c fallback identical. */
/* agg4 ADOPT of saffah_codex_6a_agg3 #140076 <https://duck.ac/submission/140076> (283,842.545us, 2026-10-23 00:51:39):
   body adopted in full with attribution (their header below intact; their chain = our #139932 <https://duck.ac/submission/139932>
   + their #139747/#140059 lazy inverse boundaries). Nothing of ours removed.
   Gates: full 80,000,004-byte output byte-identical to #139932, random and all-9 inputs.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140059: pointer-hoisted engine with three lazy middle inverse boundaries, all prior attribution below.
Idea: Extend the range proof across the padded16384-coefficient slab boundary. In the known even,L4 path, only child slabs1/2/3 feed Montgomery in merge4_dit_s3, so their H4096 outputs can remain below2P. Child0 stays canonical and every other geometry retains its original path. No padding element becomes data. Shorten private names solely for source size.
Purpose: Delete canonicalization before the next Montgomery consumer. Verify padded groups, modular ranges and complete20M outputs at three alignments.
*/
/* References: saffah_v41f1017_agg4 https://duck.ac/submission/139932: pointer-hoisted engine below; saffah_codex_6a_agg3 https://duck.ac/submission/139747: our three proven lazy inverse boundaries.
Idea: Combine those independent improvements. Keep canonical unmultiplied children, but omit normalization for fused64 and middle H64/H256/H1024 outputs consumed only by Montgomery. Bounds a<2P,b<P and unchanged canonical consumers still apply. Shorten three private helper names solely for source size.
Purpose: Validate stage residues/bounds and every full output on the new pointer layout.
*/
/* 1002e7: ac4 FUSE-IM (inv->merge; ours #139464<https://duck.ac/submission/139464>) + agg4 H1FUSE (h1->h4 DIT). Sources per RULES sec.3 (no license; acct+URL+content): saffah_codex_6a_agg3 #139373 <https://duck.ac/submission/139373> base; #139367 <https://duck.ac/submission/139367> merge5; #139358 <https://duck.ac/submission/139358>; #139353 <https://duck.ac/submission/139353>; #133546 <https://duck.ac/submission/133546> slabs; #128549 <https://duck.ac/submission/128549> changes; #125868 <https://duck.ac/submission/125868>; #121456 <https://duck.ac/submission/121456>; #121443 <https://duck.ac/submission/121443>; #121420 <https://duck.ac/submission/121420>; #120684 <https://duck.ac/submission/120684>; #120621 <https://duck.ac/submission/120621>; saffah_cc_v41_agg1 #130191 <https://duck.ac/submission/130191> base; #127195 <https://duck.ac/submission/127195>; #121257 <https://duck.ac/submission/121257>; #122356 <https://duck.ac/submission/122356>; #122333 <https://duck.ac/submission/122333>; #122392 <https://duck.ac/submission/122392>; #126906 <https://duck.ac/submission/126906>; #111674 <https://duck.ac/submission/111674>.
   agg4 chain: H1FUSE = ours #139571 <https://duck.ac/submission/139571> (inverse h1->h4 DIT ring fusion).
   agg4 additions on that base (all value-preserving; full 80,000,004-byte outputs byte-identical):
   (1) S64FUSE: inverse S=16 -> S=64 ring levels fused in registers (removes one 64KB-slab store/load round trip);
   (2) DISPHOIST: loop-invariant switch(S) dispatch hoisted out of the ring block loops in fwd_sb_ring and ntt_inv_tail;
   (3) PADGATE: PADJUMP/PADJUMP_S four boundary compares replaced by one nearest-boundary gate (min counter), sync at crossings;
   (4) PMSPLIT: dh4 g_pmdelta branch split into norm/pm loop variants with delta folded into the base pointer;
   (5) S64 kernel always_inline into ntt_inv_tail.
   (6) PTRHOIST: ds2/s3 and merge4_dit_s2/s3 tails, stage16_dif_top_batch, and the H16 ring level use raw incremented data/twiddle pointers instead of gcc9 per-iteration pointer-chain reloads; byte-exact, Ir -38.90M vs cand_disp.
   (7) PTRHOIST2: single-base table/data addressing from build_tabs layout (MP=H+TPAD) for stage4_dif_h, merge4_dit_h, ds2/s3, merge4_dit_s2/s3, fused inverse l-tables and bodyC; H4 pointer loops; H16 noinline single-base loop with its four prefetches removed. Ir -66.45M vs cand_disp (additional -27.56M vs PTRHOIST).
   Gates rand sha1 4bc5ad88828cf2aefc06810ad4032ec5c2ef1677, all-9 9276aacf2a6b58674a5d65128fd2908318e308f0.
 */
#define NTT_MAXN (1u<<25)
#pragma GCC optimize("O3,rename-registers")
#include <immintrin.h>
#include <cstdint>
#include <cstring>
#include <cstdlib>
__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));
}
typedef uint32_t u32;
typedef uint64_t u64;

#define PF_DIST 192


#define NTP g_mod
#define TPAD 704u
#define TOPPAD 128u
static inline u32 padpos(u32 i){
return i+TOPPAD*((i>>20)+(i>>18)+(i>>16)+(i>>14));
}

#define PADJUMP()                   if (j >= nxt) { if (j >= nxt0) { bp += TOPPAD; nxt0 += (1u << 20); } \
if (j >= nxt1) { bp += TOPPAD; nxt1 += (1u << 18); } \
if (j >= nxt2) { bp += TOPPAD; nxt2 += (1u << 16); } \
if (j >= nxt3) { bp += TOPPAD; nxt3 += (1u << 14); } \
nxt = nxt0; if (nxt1 < nxt) nxt = nxt1; if (nxt2 < nxt) nxt = nxt2; if (nxt3 < nxt) nxt = nxt3; }
#define PADJUMP_S()                   if (s >= nxt) { if (s >= nxt0) { bp += TOPPAD; nxt0 += (1u << 20); } \
if (s >= nxt1) { bp += TOPPAD; nxt1 += (1u << 18); } \
if (s >= nxt2) { bp += TOPPAD; nxt2 += (1u << 16); } \
if (s >= nxt3) { bp += TOPPAD; nxt3 += (1u << 14); } \
nxt = nxt0; if (nxt1 < nxt) nxt = nxt1; if (nxt2 < nxt) nxt = nxt2; if (nxt3 < nxt) nxt = nxt3; }
#define TGT __attribute__((target("avx2")))
static u32 g_scale_digit[10];
static u32 g_scale_lo[8]__attribute__((aligned(32)));
static u32 g_scale_hi[8]__attribute__((aligned(32)));


#ifndef NTT_MAXN
#define NTT_MAXN (1u << 25)
#endif
static u32 g_n,g_Sb;
static u32 g_pinv;
static int g_log2n,g_L,g_odd,g_nb,g_BT=13;
static u32 g_ciF,g_cisF,g_ciI,g_cisI;
static u32 *g_tabF;

#define Q_T5OFF (5 * NTT_MAXN / 2 + 4096 - 4 * (1u << 22))
#define Q_SCF (g_tabF + 2*g_n + 131072u)
#define Q_SCI (g_tabI + 2*g_n + 131072u)
static u32 *g_tabI;
struct Tabs{
const u32*l1[40],*l2[40],*l3[40],*m1[40],*m2[40],*m3[40];
const u32*b2,*b2s;
const u32*h1[40],*h2[40],*h3[40],*q1[40],*q2[40],*q3[40];
}
;
static Tabs g_TF,g_TI;
static inline u32 mulmod(u32 a,u32 b){
return(u32)((u64)a*b%NTP);
}
static u32 powmod(u32 a,u32 e){
u32 r=1;
while(e){
if(e&1)r=mulmod(r,a);
a=mulmod(a,a);
e>>=1;
}
return r;
}

#define RRW ((u32)((1ull<<32)%NTP))
static u32 g_wf,g_wif;


static inline u32 shoup_c(u32 w){
return(u32)(((u64)w<<32)/NTP);
}
TGT static inline __m256i vmulhi(__m256i a,__m256i b){
__m256i e=_mm256_mul_epu32(a,b);
__m256i o=_mm256_mul_epu32(_mm256_bsrli_epi128(a,4),_mm256_bsrli_epi128(b,4));
e=_mm256_bsrli_epi128(e,4);
return _mm256_blend_epi32(e,o,0xAA);
}
TGT static inline __m256i vmont(__m256i a,__m256i b,__m256i pv){
const __m256i pinv=_mm256_set1_epi32((int)g_pinv);
__m256i tlo=_mm256_mullo_epi32(a,b);
__m256i thi=vmulhi(a,b);
__m256i m=_mm256_mullo_epi32(tlo,pinv);
__m256i mhi=vmulhi(m,pv);
__m256i nz=_mm256_andnot_si256(_mm256_cmpeq_epi32(tlo,_mm256_setzero_si256()),_mm256_set1_epi32(1));
__m256i u=_mm256_add_epi32(_mm256_add_epi32(thi,mhi),nz);
return _mm256_min_epu32(u,_mm256_sub_epi32(u,pv));
}
TGT static inline __m256i mpa(__m256i a,__m256i b,__m256i pv){
const __m256i np=_mm256_set1_epi64x((long long)g_pinv);
const __m256i p64=_mm256_set1_epi64x((long long)NTP);
__m256i ae=_mm256_mul_epu32(a,b);
__m256i ao=_mm256_mul_epu32(_mm256_bsrli_epi128(a,4),_mm256_bsrli_epi128(b,4));
__m256i me=_mm256_mul_epu32(ae,np),mo=_mm256_mul_epu32(ao,np);
ae=_mm256_add_epi64(ae,_mm256_mul_epu32(me,p64));
ao=_mm256_add_epi64(ao,_mm256_mul_epu32(mo,p64));
__m256i value=_mm256_or_si256(_mm256_bsrli_epi128(ae,4),ao);
return _mm256_min_epu32(value,_mm256_sub_epi32(value,pv));
}
TGT static inline __m256i mpb(__m256i a,__m256i b,__m256i bpi,__m256i pv){
const __m256i p64=_mm256_set1_epi64x((long long)NTP);
__m256i ae=_mm256_mul_epu32(a,b);
__m256i ao=_mm256_mul_epu32(_mm256_bsrli_epi128(a,4),_mm256_bsrli_epi128(b,4));
__m256i me=_mm256_mul_epu32(a,bpi),mo=_mm256_mul_epu32(_mm256_bsrli_epi128(a,4),_mm256_bsrli_epi128(bpi,4));
ae=_mm256_add_epi64(ae,_mm256_mul_epu32(me,p64));
ao=_mm256_add_epi64(ao,_mm256_mul_epu32(mo,p64));
__m256i value=_mm256_or_si256(_mm256_bsrli_epi128(ae,4),ao);
return _mm256_min_epu32(value,_mm256_sub_epi32(value,pv));
}
TGT static inline __m128i mpc(__m128i a,__m128i b,__m128i bpi,__m128i pv){
const __m128i p64=_mm_set1_epi64x((long long)NTP);
__m128i ae=_mm_mul_epu32(a,b);
__m128i ao=_mm_mul_epu32(_mm_bsrli_si128(a,4),_mm_bsrli_si128(b,4));
__m128i me=_mm_mul_epu32(a,bpi),mo=_mm_mul_epu32(_mm_bsrli_si128(a,4),_mm_bsrli_si128(bpi,4));
ae=_mm_add_epi64(ae,_mm_mul_epu32(me,p64));
ao=_mm_add_epi64(ao,_mm_mul_epu32(mo,p64));
__m128i value=_mm_or_si128(_mm_bsrli_si128(ae,4),ao);
return _mm_min_epu32(value,_mm_sub_epi32(value,pv));
}


#define MLO_V g_shmagic
TGT static inline __m256i der_ws(__m256i w){
 __m256i q=_mm256_add_epi32(_mm256_mullo_epi32(w,_mm256_set1_epi32(g_shint)),vmulhi(w,_mm256_set1_epi32((int)g_shmagic)));
 __m256i q1=_mm256_add_epi32(q,_mm256_set1_epi32(1)),p=_mm256_set1_epi32((int)NTP);
 __m256i high=vmulhi(q1,p),low=_mm256_mullo_epi32(q1,p),zero=_mm256_setzero_si256();
 __m256i correction=_mm256_or_si256(_mm256_xor_si256(_mm256_cmpeq_epi32(high,w),_mm256_set1_epi32(-1)),_mm256_cmpeq_epi32(low,zero));
 return _mm256_add_epi32(q,_mm256_and_si256(correction,_mm256_set1_epi32(1)));
}
TGT static inline __m256i vS(__m256i a,__m256i w,__m256i ws,__m256i pv){
__m256i t=_mm256_mullo_epi32(a,w);
__m256i q=vmulhi(a,ws);
__m256i r=_mm256_sub_epi32(t,_mm256_mullo_epi32(q,pv));
return _mm256_min_epu32(r,_mm256_sub_epi32(r,pv));
}
TGT static inline __m256i vaddm(__m256i a,__m256i b,__m256i pv){
__m256i s=_mm256_add_epi32(a,b);
return _mm256_min_epu32(s,_mm256_sub_epi32(s,pv));
}
TGT static inline __m256i vsubm(__m256i a,__m256i b,__m256i pv){
__m256i d=_mm256_sub_epi32(a,b);
return _mm256_min_epu32(d,_mm256_add_epi32(d,pv));
}
TGT static inline __m128i vmulhi4(__m128i a,__m128i b){
__m128i e=_mm_mul_epu32(a,b);
__m128i o=_mm_mul_epu32(_mm_bsrli_si128(a,4),_mm_bsrli_si128(b,4));
e=_mm_bsrli_si128(e,4);
return _mm_blend_epi32(e,o,0xAA);
}
TGT static inline __m128i sh4(__m128i a,__m128i w,__m128i ws,__m128i pv){
__m128i t=_mm_mullo_epi32(a,w);
__m128i q=vmulhi4(a,ws);
__m128i r=_mm_sub_epi32(t,_mm_mullo_epi32(q,pv));
return _mm_min_epu32(r,_mm_sub_epi32(r,pv));
}
TGT static inline __m128i add4(__m128i a,__m128i b,__m128i pv){
__m128i s=_mm_add_epi32(a,b);
return _mm_min_epu32(s,_mm_sub_epi32(s,pv));
}
TGT static inline __m128i sub4(__m128i a,__m128i b,__m128i pv){
__m128i d=_mm_sub_epi32(a,b);
return _mm_min_epu32(d,_mm_add_epi32(d,pv));
}
TGT static inline __m256i ld2(const __m128i*lo,const __m128i*hi){
return _mm256_inserti128_si256(_mm256_castsi128_si256(_mm_loadu_si128(lo)),_mm_loadu_si128(hi),1);
}
TGT static inline void st2(__m128i*lo,__m128i*hi,__m256i v){
_mm_storeu_si128(lo,_mm256_castsi256_si128(v));
_mm_storeu_si128(hi,_mm256_extracti128_si256(v,1));
}
TGT static inline __m256i bc128(__m128i v){
return _mm256_inserti128_si256(_mm256_castsi128_si256(v),v,1);
}
TGT static inline __m256i vL(__m256i a,__m256i b,__m256i pv){
return _mm256_add_epi32(_mm256_sub_epi32(a,b),pv);
}
TGT static inline __m128i sub_lazy4(__m128i a,__m128i b,__m128i pv){
return _mm_add_epi32(_mm_sub_epi32(a,b),pv);
}






TGT static void stage4_dif_t(u32*a,u32 blk,u32 h,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
for(u32 i=0;
i<blk;
i+=4*h){
u32*p=a+i;
u32 j=0;
for(;
j+16<=h;
j+=16){
__m256i x0=ldu256((p+j));
__m256i x1=ldu256((p+j+h));
__m256i x2=ldu256((p+j+2*h));
__m256i x3=ldu256((p+j+3*h));
__m256i y0=ldu256((p+j+8));
__m256i y1=ldu256((p+j+h+8));
__m256i y2=ldu256((p+j+2*h+8));
__m256i y3=ldu256((p+j+3*h+8));
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i A2=vaddm(y0,y2,pv),B2=vsubm(y0,y2,pv);
__m256i C2=vaddm(y1,y3,pv),D2=vL(y1,y3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i v0=vaddm(A2,C2,pv),v2=vL(A2,C2,pv);
__m256i E=vS(D,cI,cIs,pv);
__m256i F=vS(D2,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
__m256i v1=_mm256_add_epi32(B2,F),v3=vL(B2,F,pv);
{
u1=mpa(u1,ldu256((t1+j)),pv);
}
{
u2=mpa(u2,ldu256((t2+j)),pv);
}
{
u3=mpa(u3,ldu256((t3+j)),pv);
}
{
v1=mpa(v1,ldu256((t1+j+8)),pv);
}
{
v2=mpa(v2,ldu256((t2+j+8)),pv);
}
{
v3=mpa(v3,ldu256((t3+j+8)),pv);
}
stu256((p+j),u0);
stu256((p+j+h),u1);
stu256((p+j+2*h),u2);
stu256((p+j+3*h),u3);
stu256((p+j+8),v0);
stu256((p+j+h+8),v1);
stu256((p+j+2*h+8),v2);
stu256((p+j+3*h+8),v3);
}
for(;
j+8<=h;
j+=8){
__builtin_prefetch(p+j+192);
__builtin_prefetch(p+j+h+192);
__builtin_prefetch(p+j+2*h+192);
__builtin_prefetch(p+j+3*h+192);
__m256i x0=ldu256((p+j));
__m256i x1=ldu256((p+j+h));
__m256i x2=ldu256((p+j+2*h));
__m256i x3=ldu256((p+j+3*h));
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i E=vS(D,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
{
u1=mpa(u1,ldu256((t1+j)),pv);
}
{
u2=mpa(u2,ldu256((t2+j)),pv);
}
{
u3=mpa(u3,ldu256((t3+j)),pv);
}
stu256((p+j),u0);
stu256((p+j+h),u1);
stu256((p+j+2*h),u2);
stu256((p+j+3*h),u3);
}
}
}
template<int H>TGT static void stage4_dif_h(u32*a,u32 blk,u32 _h,const u32*tb,u32 ci,u32 cis){
const u32 MP=H+TPAD;
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
for(u32 i=0;
i<blk;
i+=4*H){
u32*q0=a+i;
const u32*w=tb;
u32 j=0;
for(;
j+16<=H;
j+=16,q0+=16,w+=16){
__builtin_prefetch(q0+192);
__builtin_prefetch(q0+H+192);
__builtin_prefetch(q0+2*H+192);
__builtin_prefetch(q0+3*H+192);
__m256i x0=ldu256(q0),x1=ldu256(q0+H),x2=ldu256(q0+2*H),x3=ldu256(q0+3*H);
__m256i y0=ldu256(q0+8),y1=ldu256(q0+H+8),y2=ldu256(q0+2*H+8),y3=ldu256(q0+3*H+8);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i A2=vaddm(y0,y2,pv),B2=vsubm(y0,y2,pv);
__m256i C2=vaddm(y1,y3,pv),D2=vL(y1,y3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i v0=vaddm(A2,C2,pv),v2=vL(A2,C2,pv);
__m256i E=vS(D,cI,cIs,pv),F=vS(D2,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
__m256i v1=_mm256_add_epi32(B2,F),v3=vL(B2,F,pv);
u1=mpb(u1,ldu256(w),ldu256(w+3*MP),pv);
u2=mpb(u2,ldu256(w+MP),ldu256(w+4*MP),pv);
u3=mpb(u3,ldu256(w+2*MP),ldu256(w+5*MP),pv);
v1=mpb(v1,ldu256(w+8),ldu256(w+3*MP+8),pv);
v2=mpb(v2,ldu256(w+MP+8),ldu256(w+4*MP+8),pv);
v3=mpb(v3,ldu256(w+2*MP+8),ldu256(w+5*MP+8),pv);
stu256(q0,u0);
stu256(q0+H,u1);
stu256(q0+2*H,u2);
stu256(q0+3*H,u3);
stu256(q0+8,v0);
stu256(q0+H+8,v1);
stu256(q0+2*H+8,v2);
stu256(q0+3*H+8,v3);
}
for(;
j+8<=H;
j+=8,q0+=8,w+=8){
__builtin_prefetch(q0+192);
__builtin_prefetch(q0+H+192);
__builtin_prefetch(q0+2*H+192);
__builtin_prefetch(q0+3*H+192);
__m256i x0=ldu256(q0),x1=ldu256(q0+H),x2=ldu256(q0+2*H),x3=ldu256(q0+3*H);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i E=vS(D,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
u1=mpb(u1,ldu256(w),ldu256(w+3*MP),pv);
u2=mpb(u2,ldu256(w+MP),ldu256(w+4*MP),pv);
u3=mpb(u3,ldu256(w+2*MP),ldu256(w+5*MP),pv);
stu256(q0,u0);
stu256(q0+H,u1);
stu256(q0+2*H,u2);
stu256(q0+3*H,u3);
}
}
}

TGT static inline void transpose4(__m256i&a,__m256i&b,__m256i&c,__m256i&d){
__m256i t0=_mm256_unpacklo_epi32(a,b),t1=_mm256_unpackhi_epi32(a,b),t2=_mm256_unpacklo_epi32(c,d),t3=_mm256_unpackhi_epi32(c,d);
a=_mm256_unpacklo_epi64(t0,t2);
b=_mm256_unpackhi_epi64(t0,t2);
c=_mm256_unpacklo_epi64(t1,t3);
d=_mm256_unpackhi_epi64(t1,t3);
}
TGT static void dh1(u32*a,u32 blk,u32 ci,u32 cis);
static long g_pmdelta=0;

#define E7Z3_PM_ON() (g_pmdelta != 0)



TGT static inline __m256i e7z3_minn(__m256i u,__m256i pv){
return _mm256_min_epu32(u,_mm256_sub_epi32(u,pv));
}
TGT static inline __m256i e7z3_pmn(__m256i u,const u32*p,__m256i pv){
return mpa(ld2((const __m128i*)(p),(const __m128i*)(p+16)),u,pv);
}
TGT __attribute__((always_inline))static inline void dh4(u32*a,u32 nsub,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
if(nsub<32){
const __m128i pv4=_mm256_castsi256_si128(pv),cI4=_mm256_castsi256_si128(cI),cIs4=_mm256_castsi256_si128(cIs);
__m128i W1=_mm_loadu_si128((const __m128i*)t1);
__m128i WS1=_mm_loadu_si128((const __m128i*)s1);
__m128i W2=_mm_loadu_si128((const __m128i*)t2);
__m128i WS2=_mm_loadu_si128((const __m128i*)s2);
__m128i W3=_mm_loadu_si128((const __m128i*)t3);
__m128i WS3=_mm_loadu_si128((const __m128i*)s3);
for(u32 i=0;
i<nsub;
i+=16){
u32*p=a+i;
__m128i x0=_mm_loadu_si128((const __m128i*)(p));
__m128i x1=_mm_loadu_si128((const __m128i*)(p+4));
__m128i x2=_mm_loadu_si128((const __m128i*)(p+8));
__m128i x3=_mm_loadu_si128((const __m128i*)(p+12));
__m128i A=add4(x0,x2,pv4),B=sub4(x0,x2,pv4);
__m128i C=add4(x1,x3,pv4),D=sub_lazy4(x1,x3,pv4);
__m128i u0=add4(A,C,pv4),u2=sub_lazy4(A,C,pv4);
__m128i E=sh4(D,cI4,cIs4,pv4);
__m128i u1=_mm_add_epi32(B,E),u3=sub_lazy4(B,E,pv4);
u1=mpc(u1,W1,WS1,pv4);
u2=mpc(u2,W2,WS2,pv4);
u3=mpc(u3,W3,WS3,pv4);
_mm_storeu_si128((__m128i*)(p),u0);
_mm_storeu_si128((__m128i*)(p+4),u1);
_mm_storeu_si128((__m128i*)(p+8),u2);
_mm_storeu_si128((__m128i*)(p+12),u3);
}
dh1(a,nsub,ci,cis);
return;
}
__m256i W1=bc128(_mm_loadu_si128((const __m128i*)t1));
__m256i WS1=bc128(_mm_loadu_si128((const __m128i*)s1));
__m256i W2=bc128(_mm_loadu_si128((const __m128i*)t2));
__m256i WS2=bc128(_mm_loadu_si128((const __m128i*)s2));
__m256i W3=bc128(_mm_loadu_si128((const __m128i*)t3));
__m256i WS3=bc128(_mm_loadu_si128((const __m128i*)s3));
if(g_pmdelta==0){

#pragma GCC unroll 2
for(u32 i=0,*p=a;
i<nsub;
i+=32,p+=32){
__m256i x0=ld2((const __m128i*)(p),(const __m128i*)(p+16));
__m256i x1=ld2((const __m128i*)(p+4),(const __m128i*)(p+20));
__m256i x2=ld2((const __m128i*)(p+8),(const __m128i*)(p+24));
__m256i x3=ld2((const __m128i*)(p+12),(const __m128i*)(p+28));
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i E=vS(D,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
u1=mpb(u1,W1,WS1,pv);
u2=mpb(u2,W2,WS2,pv);
u3=mpb(u3,W3,WS3,pv);
transpose4(u0,u1,u2,u3);
{
__m256i A4=vaddm(u0,u2,pv),B4=vsubm(u0,u2,pv),C4=vaddm(u1,u3,pv),D4=vL(u1,u3,pv),E4=vS(D4,cI,cIs,pv);
u0=vaddm(A4,C4,pv);
u2=vL(A4,C4,pv);
u1=_mm256_add_epi32(B4,E4);
u3=vL(B4,E4,pv);
transpose4(u0,u1,u2,u3);
}
u0=e7z3_minn(u0,pv);
u1=e7z3_minn(u1,pv);
u2=e7z3_minn(u2,pv);
u3=e7z3_minn(u3,pv);
st2((__m128i*)(p),(__m128i*)(p+16),u0);
st2((__m128i*)(p+4),(__m128i*)(p+20),u1);
st2((__m128i*)(p+8),(__m128i*)(p+24),u2);
st2((__m128i*)(p+12),(__m128i*)(p+28),u3);
}
}
else{

#pragma GCC unroll 2
for(u32 i=0,*p=a,*pd=(u32*)((char*)a+g_pmdelta);
i<nsub;
i+=32,p+=32,pd+=32){
__m256i x0=ld2((const __m128i*)(p),(const __m128i*)(p+16));
__m256i x1=ld2((const __m128i*)(p+4),(const __m128i*)(p+20));
__m256i x2=ld2((const __m128i*)(p+8),(const __m128i*)(p+24));
__m256i x3=ld2((const __m128i*)(p+12),(const __m128i*)(p+28));
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i E=vS(D,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
u1=mpb(u1,W1,WS1,pv);
u2=mpb(u2,W2,WS2,pv);
u3=mpb(u3,W3,WS3,pv);
transpose4(u0,u1,u2,u3);
{
__m256i A4=vaddm(u0,u2,pv),B4=vsubm(u0,u2,pv),C4=vaddm(u1,u3,pv),D4=vL(u1,u3,pv),E4=vS(D4,cI,cIs,pv);
u0=vaddm(A4,C4,pv);
u2=vL(A4,C4,pv);
u1=_mm256_add_epi32(B4,E4);
u3=vL(B4,E4,pv);
transpose4(u0,u1,u2,u3);
}
u0=e7z3_pmn(u0,pd,pv);
u1=e7z3_pmn(u1,pd+4,pv);
u2=e7z3_pmn(u2,pd+8,pv);
u3=e7z3_pmn(u3,pd+12,pv);
st2((__m128i*)(pd),(__m128i*)(pd+16),u0);
st2((__m128i*)(pd+4),(__m128i*)(pd+20),u1);
st2((__m128i*)(pd+8),(__m128i*)(pd+24),u2);
st2((__m128i*)(pd+12),(__m128i*)(pd+28),u3);
}
}
}
TGT static void dh1(u32*a,u32 blk,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
for(u32 i=0;
i<blk;
i+=8){
u32*p=a+i;
__m256i x=ldu256(p);
__m256i t=_mm256_shuffle_epi32(x,_MM_SHUFFLE(1,0,3,2));
__m256i s=vaddm(x,t,pv);
__m256i d=vsubm(x,t,pv);
__m256i Z=_mm256_shuffle_epi32(d,0x00);
__m256i P=_mm256_blend_epi32(s,Z,0xAA);
__m256i E=vS(d,cI,cIs,pv);
__m256i E2=_mm256_shuffle_epi32(E,0x55);
__m256i W=_mm256_shuffle_epi32(s,0x55);
__m256i Q=_mm256_blend_epi32(W,E2,0xAA);
__m256i sum=vaddm(P,Q,pv);
__m256i dif=vsubm(P,Q,pv);
stu256(p,_mm256_blend_epi32(sum,dif,0xCC));
}
}
TGT static void stage2_dif(u32*a,u32 blk,const u32*t,const u32*s){
const __m256i pv=_mm256_set1_epi32((int)NTP);
u32 h=blk>>1;
for(u32 j=0;
j<h;
j+=8){
__m256i x0=ldu256((a+j));
__m256i x1=ldu256((a+j+h));
__m256i u0=vaddm(x0,x1,pv);
__m256i u1=mpa(vsubm(x0,x1,pv),ldu256((t+j)),pv);
stu256((a+j),u0);
stu256((a+j+h),u1);
}
}



TGT static void merge4_dit_t(u32*a,u32 blk,u32 h,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
for(u32 i=0;
i<blk;
i+=4*h){
u32*p=a+i;
for(u32 j=0;
j<h;
j+=8){
__m256i x0=ldu256((p+j));
__m256i x1=ldu256((p+j+h));
__m256i x2=ldu256((p+j+2*h));
__m256i x3=ldu256((p+j+3*h));
{
x1=mpa(x1,ldu256((t1+j)),pv);
}
{
x2=mpa(x2,ldu256((t2+j)),pv);
}
{
x3=mpa(x3,ldu256((t3+j)),pv);
}
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i E=vS(D,cI,cIs,pv);
stu256((p+j),vaddm(A,C,pv));
stu256((p+j+h),vaddm(B,E,pv));
stu256((p+j+2*h),vsubm(A,C,pv));
stu256((p+j+3*h),vsubm(B,E,pv));
}
}
}
template<int H,bool LazyOutput=false>TGT static void merge4_dit_h(u32*a,u32 blk,u32 _h,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis){
const u32 MP=H+TPAD;
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
(void)t2;
(void)t3;
(void)s1;
(void)s2;
(void)s3;
for(u32 i=0;
i<blk;
i+=4*H){
u32*q0=a+i;
const u32*w=t1;
u32 j=0;
for(;
j<H;
j+=16,q0+=16,w+=16){
__builtin_prefetch(q0+192);
__builtin_prefetch(q0+H+192);
__builtin_prefetch(q0+2*H+192);
__builtin_prefetch(q0+3*H+192);
__m256i x0=ldu256(q0),y0=ldu256(q0+8);
__m256i x1=ldu256(q0+H),y1=ldu256(q0+H+8);
__m256i x2=ldu256(q0+2*H),y2=ldu256(q0+2*H+8);
__m256i x3=ldu256(q0+3*H),y3=ldu256(q0+3*H+8);
x1=mpb(x1,ldu256(w),ldu256(w+3*MP),pv);
y1=mpb(y1,ldu256(w+8),ldu256(w+3*MP+8),pv);
x2=mpb(x2,ldu256(w+MP),ldu256(w+4*MP),pv);
y2=mpb(y2,ldu256(w+MP+8),ldu256(w+4*MP+8),pv);
x3=mpb(x3,ldu256(w+2*MP),ldu256(w+5*MP),pv);
y3=mpb(y3,ldu256(w+2*MP+8),ldu256(w+5*MP+8),pv);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i A2=vaddm(y0,y2,pv),B2=vsubm(y0,y2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i C2=vaddm(y1,y3,pv),D2=vL(y1,y3,pv);
__m256i E=vS(D,cI,cIs,pv),F=vS(D2,cI,cIs,pv);
stu256(q0,(LazyOutput?_mm256_add_epi32(A,C):vaddm(A,C,pv)));
stu256(q0+8,(LazyOutput?_mm256_add_epi32(A2,C2):vaddm(A2,C2,pv)));
stu256(q0+H,(LazyOutput?_mm256_add_epi32(B,E):vaddm(B,E,pv)));
stu256(q0+H+8,(LazyOutput?_mm256_add_epi32(B2,F):vaddm(B2,F,pv)));
stu256(q0+2*H,(LazyOutput?vL(A,C,pv):vsubm(A,C,pv)));
stu256(q0+2*H+8,(LazyOutput?vL(A2,C2,pv):vsubm(A2,C2,pv)));
stu256(q0+3*H,(LazyOutput?vL(B,E,pv):vsubm(B,E,pv)));
stu256(q0+3*H+8,(LazyOutput?vL(B2,F,pv):vsubm(B2,F,pv)));
}
for(;
j<H;
j+=8,q0+=8,w+=8){
__m256i x0=ldu256(q0),x1=ldu256(q0+H),x2=ldu256(q0+2*H),x3=ldu256(q0+3*H);
x1=mpb(x1,ldu256(w),ldu256(w+3*MP),pv);
x2=mpb(x2,ldu256(w+MP),ldu256(w+4*MP),pv);
x3=mpb(x3,ldu256(w+2*MP),ldu256(w+5*MP),pv);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i E=vS(D,cI,cIs,pv);
stu256(q0,(LazyOutput?_mm256_add_epi32(A,C):vaddm(A,C,pv)));
stu256(q0+H,(LazyOutput?_mm256_add_epi32(B,E):vaddm(B,E,pv)));
stu256(q0+2*H,(LazyOutput?vL(A,C,pv):vsubm(A,C,pv)));
stu256(q0+3*H,(LazyOutput?vL(B,E,pv):vsubm(B,E,pv)));
}
}
}

TGT static void merge4_dit_h4(u32*a,u32 nsub,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
if(nsub<32){
const __m128i pv4=_mm256_castsi256_si128(pv),cI4=_mm256_castsi256_si128(cI),cIs4=_mm256_castsi256_si128(cIs);
__m128i W1=_mm_loadu_si128((const __m128i*)t1);
__m128i WS1=_mm_loadu_si128((const __m128i*)s1);
__m128i W2=_mm_loadu_si128((const __m128i*)t2);
__m128i WS2=_mm_loadu_si128((const __m128i*)s2);
__m128i W3=_mm_loadu_si128((const __m128i*)t3);
__m128i WS3=_mm_loadu_si128((const __m128i*)s3);
for(u32 i=0;
i<nsub;
i+=16){
u32*p=a+i;
__m128i x0=_mm_loadu_si128((const __m128i*)(p));
__m128i x1=mpc(_mm_loadu_si128((const __m128i*)(p+4)),W1,WS1,pv4);
__m128i x2=mpc(_mm_loadu_si128((const __m128i*)(p+8)),W2,WS2,pv4);
__m128i x3=mpc(_mm_loadu_si128((const __m128i*)(p+12)),W3,WS3,pv4);
__m128i A=add4(x0,x2,pv4),B=sub4(x0,x2,pv4);
__m128i C=add4(x1,x3,pv4),D=sub4(x1,x3,pv4);
__m128i E=sh4(D,cI4,cIs4,pv4);
_mm_storeu_si128((__m128i*)(p),add4(A,C,pv4));
_mm_storeu_si128((__m128i*)(p+4),add4(B,E,pv4));
_mm_storeu_si128((__m128i*)(p+8),sub4(A,C,pv4));
_mm_storeu_si128((__m128i*)(p+12),sub4(B,E,pv4));
}
return;
}
__m256i W1=bc128(_mm_loadu_si128((const __m128i*)t1));
__m256i WS1=bc128(_mm_loadu_si128((const __m128i*)s1));
__m256i W2=bc128(_mm_loadu_si128((const __m128i*)t2));
__m256i WS2=bc128(_mm_loadu_si128((const __m128i*)s2));
__m256i W3=bc128(_mm_loadu_si128((const __m128i*)t3));
__m256i WS3=bc128(_mm_loadu_si128((const __m128i*)s3));
for(u32 i=0;
i<nsub;
i+=32){
u32*p=a+i;
__m256i x0=ld2((const __m128i*)(p),(const __m128i*)(p+16));
__m256i x1=ld2((const __m128i*)(p+4),(const __m128i*)(p+20));
__m256i x2=ld2((const __m128i*)(p+8),(const __m128i*)(p+24));
__m256i x3=ld2((const __m128i*)(p+12),(const __m128i*)(p+28));
x1=mpb(x1,W1,WS1,pv);
x2=mpb(x2,W2,WS2,pv);
x3=mpb(x3,W3,WS3,pv);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i E=vS(D,cI,cIs,pv);
st2((__m128i*)(p),(__m128i*)(p+16),vaddm(A,C,pv));
st2((__m128i*)(p+4),(__m128i*)(p+20),vaddm(B,E,pv));
st2((__m128i*)(p+8),(__m128i*)(p+24),vsubm(A,C,pv));
st2((__m128i*)(p+12),(__m128i*)(p+28),vsubm(B,E,pv));
}
}
TGT static inline __m256i merge4_h1v(__m256i x,__m256i pv,__m256i cI,__m256i cIs){
__m256i t=_mm256_shuffle_epi32(x,_MM_SHUFFLE(1,0,3,2)),sv=vaddm(x,t,pv),d=vsubm(x,t,pv),Z=_mm256_shuffle_epi32(d,0x00),P=_mm256_blend_epi32(sv,Z,0xAA),E=vS(d,cI,cIs,pv),E2=_mm256_shuffle_epi32(E,0x55),W=_mm256_shuffle_epi32(sv,0x55),Q=_mm256_blend_epi32(W,E2,0xAA),sum=vaddm(P,Q,pv),dif=vsubm(P,Q,pv);
return _mm256_blend_epi32(sum,dif,0xCC);
}
TGT static inline __m256i merge4_h1v_lazy(__m256i x,__m256i pv,__m256i cI,__m256i cIs){
__m256i t=_mm256_shuffle_epi32(x,_MM_SHUFFLE(1,0,3,2)),sv=vaddm(x,t,pv),d=vsubm(x,t,pv),Z=_mm256_shuffle_epi32(d,0x00),P=_mm256_blend_epi32(sv,Z,0xAA),E=vS(d,cI,cIs,pv),E2=_mm256_shuffle_epi32(E,0x55),W=_mm256_shuffle_epi32(sv,0x55),Q=_mm256_blend_epi32(W,E2,0xAA),sum=_mm256_add_epi32(P,Q),dif=vL(P,Q,pv);
return _mm256_blend_epi32(sum,dif,0xCC);
}
TGT static void mhh(u32*a,u32 nsub,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP),cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
__m256i W1=bc128(_mm_loadu_si128((const __m128i*)t1)),WS1=bc128(_mm_loadu_si128((const __m128i*)s1)),W2=bc128(_mm_loadu_si128((const __m128i*)t2)),WS2=bc128(_mm_loadu_si128((const __m128i*)s2)),W3=bc128(_mm_loadu_si128((const __m128i*)t3)),WS3=bc128(_mm_loadu_si128((const __m128i*)s3));
for(u32 i=0;
i+32<=nsub;
i+=32){
u32*p=a+i;
__m256i g0=merge4_h1v(ldu256((const void*)(p)),pv,cI,cIs),g1=merge4_h1v(ldu256((const void*)(p+8)),pv,cI,cIs),g2=merge4_h1v(ldu256((const void*)(p+16)),pv,cI,cIs),g3=merge4_h1v(ldu256((const void*)(p+24)),pv,cI,cIs);
__m256i x0=_mm256_permute2x128_si256(g0,g2,0x20),x1=_mm256_permute2x128_si256(g0,g2,0x31),x2=_mm256_permute2x128_si256(g1,g3,0x20),x3=_mm256_permute2x128_si256(g1,g3,0x31);
x1=mpb(x1,W1,WS1,pv);
x2=mpb(x2,W2,WS2,pv);
x3=mpb(x3,W3,WS3,pv);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv),C=vaddm(x1,x3,pv),D=vL(x1,x3,pv),E=vS(D,cI,cIs,pv);
st2((__m128i*)(p),(__m128i*)(p+16),vaddm(A,C,pv));
st2((__m128i*)(p+4),(__m128i*)(p+20),vaddm(B,E,pv));
st2((__m128i*)(p+8),(__m128i*)(p+24),vsubm(A,C,pv));
st2((__m128i*)(p+12),(__m128i*)(p+28),vsubm(B,E,pv));
}
for(u32 i=(nsub&~31u);
i<nsub;
i+=8){
u32*p=a+i;
stu256((void*)p,merge4_h1v(ldu256((const void*)p),pv,cI,cIs));
}
}

#define H4B32(Q,R1,R2,R3,R4,LAZY) {const u32*q=(Q);__m256i g0=merge4_h1v(ldu256((const void*)(q)),pv,cI,cIs),g1=merge4_h1v_lazy(ldu256((const void*)(q+8)),pv,cI,cIs),g2=merge4_h1v(ldu256((const void*)(q+16)),pv,cI,cIs),g3=merge4_h1v_lazy(ldu256((const void*)(q+24)),pv,cI,cIs);__m256i e0=_mm256_permute2x128_si256(g0,g2,0x20),e1=_mm256_permute2x128_si256(g0,g2,0x31),e2=_mm256_permute2x128_si256(g1,g3,0x20),e3=_mm256_permute2x128_si256(g1,g3,0x31);e1=mpb(e1,W1,WS1,pv);e2=mpb(e2,W2,WS2,pv);e3=mpb(e3,W3,WS3,pv);__m256i a1=vaddm(e0,e2,pv),b1=vsubm(e0,e2,pv),c1=vaddm(e1,e3,pv),d1=vL(e1,e3,pv),f1=vS(d1,cI,cIs,pv);if constexpr(LAZY){R1=_mm256_add_epi32(a1,c1);R2=_mm256_add_epi32(b1,f1);R3=vL(a1,c1,pv);R4=vL(b1,f1,pv);}else{R1=vaddm(a1,c1,pv);R2=vaddm(b1,f1,pv);R3=vsubm(a1,c1,pv);R4=vsubm(b1,f1,pv);}}
TGT __attribute__((always_inline))static inline void mhh64(u32*a,u32 nsub,const u32*t1,const u32*t2,const u32*t3,const u32*s1,const u32*s2,const u32*s3,u32 ci,u32 cis,const u32*u1,const u32*u2,const u32*u3,const u32*v1,const u32*v2,const u32*v3){
const __m256i pv=_mm256_set1_epi32((int)NTP),cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
const __m256i W1=bc128(_mm_loadu_si128((const __m128i*)t1)),WS1=bc128(_mm_loadu_si128((const __m128i*)s1)),W2=bc128(_mm_loadu_si128((const __m128i*)t2)),WS2=bc128(_mm_loadu_si128((const __m128i*)s2)),W3=bc128(_mm_loadu_si128((const __m128i*)t3)),WS3=bc128(_mm_loadu_si128((const __m128i*)s3));
for(u32 i=0,*p=a;
i+64<=nsub;
i+=64,p+=64){
__m256i R1,R2,R3,R4,R1b,R2b,R3b,R4b;
H4B32(p,R1,R2,R3,R4,false)H4B32(p+32,R1b,R2b,R3b,R4b,true)__m256i x0=_mm256_permute2x128_si256(R1,R2,0x20),y0=_mm256_permute2x128_si256(R3,R4,0x20),x1=_mm256_permute2x128_si256(R1,R2,0x31),y1=_mm256_permute2x128_si256(R3,R4,0x31),x2=_mm256_permute2x128_si256(R1b,R2b,0x20),y2=_mm256_permute2x128_si256(R3b,R4b,0x20),x3=_mm256_permute2x128_si256(R1b,R2b,0x31),y3=_mm256_permute2x128_si256(R3b,R4b,0x31);
x1=mpb(x1,ldu256((const void*)(u1)),ldu256((const void*)(v1)),pv);
y1=mpb(y1,ldu256((const void*)(u1+8)),ldu256((const void*)(v1+8)),pv);
x2=mpb(x2,ldu256((const void*)(u2)),ldu256((const void*)(v2)),pv);
y2=mpb(y2,ldu256((const void*)(u2+8)),ldu256((const void*)(v2+8)),pv);
x3=mpb(x3,ldu256((const void*)(u3)),ldu256((const void*)(v3)),pv);
y3=mpb(y3,ldu256((const void*)(u3+8)),ldu256((const void*)(v3+8)),pv);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv),A2=vaddm(y0,y2,pv),B2=vsubm(y0,y2,pv),C=vaddm(x1,x3,pv),D=vL(x1,x3,pv),C2=vaddm(y1,y3,pv),D2=vL(y1,y3,pv),E=vS(D,cI,cIs,pv),F=vS(D2,cI,cIs,pv);
if(nsub>=256u && (nsub&255u)==0u && (i&192u)!=0u){
stu256((void*)(p),_mm256_add_epi32(A,C));
stu256((void*)(p+8),_mm256_add_epi32(A2,C2));
stu256((void*)(p+16),_mm256_add_epi32(B,E));
stu256((void*)(p+24),_mm256_add_epi32(B2,F));
stu256((void*)(p+32),vL(A,C,pv));
stu256((void*)(p+40),vL(A2,C2,pv));
stu256((void*)(p+48),vL(B,E,pv));
stu256((void*)(p+56),vL(B2,F,pv));
}
else{
stu256((void*)(p),vaddm(A,C,pv));
stu256((void*)(p+8),vaddm(A2,C2,pv));
stu256((void*)(p+16),vaddm(B,E,pv));
stu256((void*)(p+24),vaddm(B2,F,pv));
stu256((void*)(p+32),vsubm(A,C,pv));
stu256((void*)(p+40),vsubm(A2,C2,pv));
stu256((void*)(p+48),vsubm(B,E,pv));
stu256((void*)(p+56),vsubm(B2,F,pv));
}
}
for(u32 i=(nsub&~63u);
i<nsub;
i+=8){
u32*p=a+i;
stu256((void*)p,merge4_h1v(ldu256((const void*)p),pv,cI,cIs));
}
}
TGT static void mh1(u32*a,u32 blk,u32 ci,u32 cis){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
for(u32 i=0;
i<blk;
i+=8){
u32*p=a+i;
__m256i x=ldu256(p);
__m256i t=_mm256_shuffle_epi32(x,_MM_SHUFFLE(1,0,3,2));
__m256i s=vaddm(x,t,pv);
__m256i d=vsubm(x,t,pv);
__m256i Z=_mm256_shuffle_epi32(d,0x00);
__m256i P=_mm256_blend_epi32(s,Z,0xAA);
__m256i E=vS(d,cI,cIs,pv);
__m256i E2=_mm256_shuffle_epi32(E,0x55);
__m256i W=_mm256_shuffle_epi32(s,0x55);
__m256i Q=_mm256_blend_epi32(W,E2,0xAA);
__m256i sum=vaddm(P,Q,pv);
__m256i dif=vsubm(P,Q,pv);
stu256(p,_mm256_blend_epi32(sum,dif,0xCC));
}
}
TGT static void merge2_dit(u32*a,u32 blk,const u32*t,const u32*s){
const __m256i pv=_mm256_set1_epi32((int)NTP);
u32 h=blk>>1;
for(u32 j=0;
j<h;
j+=8){
__m256i x0=ldu256((a+j));
__m256i x1=mpa(ldu256((a+j+h)),ldu256((t+j)),pv);
stu256((a+j),vaddm(x0,x1,pv));
stu256((a+j+h),vsubm(x0,x1,pv));
}
}
TGT static inline void ntst256(void*p,__m256i v){
if(((uintptr_t)p&31u)==0)_mm256_stream_si256((__m256i*)p,v);
else stu256(p,v);
}
TGT static void gen_geo_val(u32*dst,u32 cnt,u32 v0,u32 r);
TGT static void gen_geo(u32*dst,u32*sc,u32 cnt,u32 v0,u32 r){
if(!sc){
gen_geo_val(dst,cnt,v0,r);
return;
}
if(!cnt)return;
if(cnt<32){
u32 x=v0%NTP;
for(u32 i=0;
i<cnt;
i++){
dst[i]=x;
sc[i]=shoup_c(x);
x=mulmod(x,r);
}
return;
}
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i mlo=_mm256_set1_epi32((int)g_shmagic);
u32 seq[32],cur=v0%NTP;
for(int i=0;
i<32;
i++){
seq[i]=cur;
cur=mulmod(cur,r);
}
if(v0!=1u)cur=powmod(r,32);
u32 rs=shoup_c(cur);
const __m256i rv=_mm256_set1_epi32((int)cur);
const __m256i rsv=_mm256_set1_epi32((int)rs);
__m256i a=ldu256(seq),b=ldu256(seq+8),c=ldu256(seq+16),d=ldu256(seq+24);
u32 i=0;
for(;
i+32<=cnt;
i+=32){
ntst256(dst+i,a);
ntst256(dst+i+8,b);
ntst256(dst+i+16,c);
ntst256(dst+i+24,d);
ntst256(sc+i,der_ws(a));
ntst256(sc+i+8,der_ws(b));
ntst256(sc+i+16,der_ws(c));
ntst256(sc+i+24,der_ws(d));
a=vS(a,rv,rsv,pv);
b=vS(b,rv,rsv,pv);
c=vS(c,rv,rsv,pv);
d=vS(d,rv,rsv,pv);
}
if(i<cnt){
u32 x=mulmod(v0,powmod(r,i));
for(;
i<cnt;
i++){
dst[i]=x;
sc[i]=shoup_c(x);
x=mulmod(x,r);
}
}
}
TGT static void gen_geo_val(u32*dst,u32 cnt,u32 v0,u32 r){
if(!cnt)return;
if(cnt<32){
u32 x=v0%NTP;
for(u32 i=0;
i<cnt;
i++){
dst[i]=x;
x=mulmod(x,r);
}
return;
}
const __m256i pv=_mm256_set1_epi32((int)NTP);
u32 seq[32],cur=v0%NTP;
for(int i=0;
i<32;
i++){
seq[i]=cur;
cur=mulmod(cur,r);
}
if(v0!=1u)cur=powmod(r,32);
if(v0!=1u)cur=powmod(r,32);
u32 rs=shoup_c(cur);
const __m256i rv=_mm256_set1_epi32((int)cur);
const __m256i rsv=_mm256_set1_epi32((int)rs);
__m256i a=ldu256(seq),b=ldu256(seq+8),c=ldu256(seq+16),d=ldu256(seq+24);
u32 i=0;
for(;
i+32<=cnt;
i+=32){
ntst256(dst+i,a);
ntst256(dst+i+8,b);
ntst256(dst+i+16,c);
ntst256(dst+i+24,d);
a=vS(a,rv,rsv,pv);
b=vS(b,rv,rsv,pv);
c=vS(c,rv,rsv,pv);
d=vS(d,rv,rsv,pv);
}
if(i<cnt){
u32 x=mulmod(v0,powmod(r,i));
for(;
i<cnt;
i++){
dst[i]=x;
x=mulmod(x,r);
}
}
}
static void build_tabs(Tabs&T,u32*tab,u32 n,u32 Sb,int L,int odd,u32 root,u32*scbase,u32 skipsc){
u32*p=tab;
u32*sp=scbase;
u32 S=n;
for(int l=0;
l<L;
l++){
u32 M=S>>2;
const u32 MPr=M+TPAD;
T .l1[l]=p;
T .l2[l]=p+MPr;
T .l3[l]=p+2*MPr;
T .m1[l]=sp;
T .m2[l]=sp+MPr;
T .m3[l]=sp+2*MPr;
u32 base=powmod(root,n/S);
u32 b1=base,b2=mulmod(base,base),b3=mulmod(b2,base);
{
const bool nosc=((skipsc>>l)&1u)!=0;
const u32 unit=(u32)((1ull<<32)%NTP);
gen_geo(p,nosc?0:sp,M,unit,b1);
gen_geo(p+MPr,nosc?0:sp+MPr,M,unit,b2);
gen_geo(p+2*MPr,nosc?0:sp+2*MPr,M,unit,b3);
}
p+=3*MPr;
sp+=3*MPr;
S>>=2;
}
T .b2=0;
T .b2s=0;
if(odd){
u32 M=Sb>>1;
T .b2=p;
T .b2s=p+M;
gen_geo(p,0,M,(u32)((1ull<<32)%NTP),powmod(root,n/Sb));
p+=2*M;
}
int s=0;
u32 S2v=4;
while(s<g_nb){
u32 M=S2v>>2;
if(M>=32)p=(u32*)(((uintptr_t)p+31u)&~(uintptr_t)31u);
const u32 MP=M+TPAD;
T .h1[s]=p;
T .h2[s]=p+MP;
T .h3[s]=p+2*MP;
T .q1[s]=p+3*MP;
T .q2[s]=p+4*MP;
T .q3[s]=p+5*MP;
u32 base=powmod(root,n/S2v);
u32 b1=base,b2=mulmod(base,base),b3=mulmod(b2,base);
gen_geo(p,0,M,(u32)((1ull<<32)%NTP),b1);
gen_geo(p+MP,0,M,(u32)((1ull<<32)%NTP),b2);
gen_geo(p+2*MP,0,M,(u32)((1ull<<32)%NTP),b3);
for(u32 i2=0;
i2<M;
i2++){
p[3*MP+i2]=(u32)((u64)p[i2]*(u64)g_pinv);
p[4*MP+i2]=(u32)((u64)p[MP+i2]*(u64)g_pinv);
p[5*MP+i2]=(u32)((u64)p[2*MP+i2]*(u64)g_pinv);
}
p+=6*MP;
S2v<<=2;
s++;
}
}
static void ntt_init(u32 n){
{
u32 iv=1;
for(int q=0;
q<5;
q++)iv*=2u-NTP*iv;
g_pinv=0u-iv;
}
g_n=n;
int lg=0;
while((1u<<lg)<n)lg++;
g_log2n=lg;
const int BT=g_BT;
int L=0;
while(lg-2*(L+1)>=BT)L++;
g_L=L;
g_Sb=n>>(2*L);
g_odd=((lg-2*L)&1)?1:0;
u32 S2=g_Sb;
if(g_odd)S2>>=1;
g_nb=0;
while(S2>=4){
S2>>=2;
g_nb++;
}
u32 w=powmod(g_generator,(NTP-1)/n);
u32 wi=powmod(w,NTP-2);
u32 iF=powmod(w,n>>2),iI=powmod(wi,n>>2);
g_ciF=iF;
g_cisF=shoup_c(iF);
g_ciI=iI;
g_cisI=shoup_c(iI);
build_tabs(g_TF,g_tabF,n,g_Sb,L,g_odd,w,Q_SCF,~0u);
build_tabs(g_TI,g_tabI,n,g_Sb,L,g_odd,wi,Q_SCI,~0u);
}
int g_fuse16=1;

#define E7_TILE 512u
static u32 g_scrA[16*E7_TILE]__attribute__((aligned(64)));


#ifndef E7Z_SBFUSE
#define E7Z_SBFUSE 1
#endif
#ifndef E7Z_IVFUSE
#define E7Z_IVFUSE 1
#endif

TGT __attribute__((always_inline))static inline void stage4_h16_loop(u32*a,u32 total,const u32*tb,u32 ci,u32 cis){
const u32 MP=16+TPAD;
const __m256i pv=_mm256_set1_epi32((int)NTP);
const __m256i cI=_mm256_set1_epi32((int)ci),cIs=_mm256_set1_epi32((int)cis);
for(u32 i=0;
i<total;
i+=64){
u32*q0=a+i;
const u32*w=tb;
u32 j=0;
for(;
j+16<=16;
j+=16,q0+=16,w+=16){
__m256i x0=ldu256(q0),x1=ldu256(q0+16),x2=ldu256(q0+32),x3=ldu256(q0+48);
__m256i y0=ldu256(q0+8),y1=ldu256(q0+24),y2=ldu256(q0+40),y3=ldu256(q0+56);
__m256i A=vaddm(x0,x2,pv),B=vsubm(x0,x2,pv);
__m256i C=vaddm(x1,x3,pv),D=vL(x1,x3,pv);
__m256i A2=vaddm(y0,y2,pv),B2=vsubm(y0,y2,pv);
__m256i C2=vaddm(y1,y3,pv),D2=vL(y1,y3,pv);
__m256i u0=vaddm(A,C,pv),u2=vL(A,C,pv);
__m256i v0=vaddm(A2,C2,pv),v2=vL(A2,C2,pv);
__m256i E=vS(D,cI,cIs,pv),F=vS(D2,cI,cIs,pv);
__m256i u1=_mm256_add_epi32(B,E),u3=vL(B,E,pv);
__m256i v1=_mm256_add_epi32(B2,F),v3=vL(B2,F,pv);
u1=mpb(u1,ldu256(w),ldu256(w+3*MP),pv);
u2=mpb(u2,ldu256(w+MP),ldu256(w+4*MP),pv);
u3=mpb(u3,ldu256(w+2*MP),ldu256(w+5*MP),pv);
v1=mpb(v1,ldu256(w+8),ldu256(w+3*MP+8),pv);
v2=mpb(v2,ldu256(w+MP+8),ldu256(w+4*MP+8),pv);
v3=mpb(v3,ldu256(w+2*MP+8),ldu256(w+5*MP+8),pv);
stu256(q0,u0);
stu256(q0+16,u1);
stu256(q0+32,u2);
stu256(q0+48,u3);
stu256(q0+8,v0);
stu256(q0+24,v1);
stu256(q0+40,v2);
stu256(q0+56,v3);
}
}
}
TGT static void fwd_sb_ring(u32*a,u32 b){
u32*p=a+padpos(b);
if(g_odd)stage2_dif(p,g_Sb,g_TF .b2,g_TF .b2s);
for(int s=g_nb-1;
s>=0;
s--){
u32 S2=4u<<(2*s);
if(S2>=32){
if(S2==16384){
for(u32 i=0;
i<g_Sb;
i+=S2)stage4_dif_h<4096>(p+i,16384,4096,g_TF .h1[s],g_ciF,g_cisF);
}
else if(S2==4096){
for(u32 i=0;
i<g_Sb;
i+=S2)stage4_dif_h<1024>(p+i,4096,1024,g_TF .h1[s],g_ciF,g_cisF);
}
else if(S2==1024){
for(u32 i=0;
i<g_Sb;
i+=S2)stage4_dif_h<256>(p+i,1024,256,g_TF .h1[s],g_ciF,g_cisF);
}
else if(S2==256){
for(u32 i=0;
i<g_Sb;
i+=S2)stage4_dif_h<64>(p+i,256,64,g_TF .h1[s],g_ciF,g_cisF);
}
else if(S2==64){
stage4_h16_loop(p,g_Sb,g_TF .h1[s],g_ciF,g_cisF);
}
else{
for(u32 i=0;
i<g_Sb;
i+=S2)stage4_dif_t(p+i,S2,S2>>2,g_TF .h1[s],g_TF .h2[s],g_TF .h3[s],g_TF .q1[s],g_TF .q2[s],g_TF .q3[s],g_ciF,g_cisF);
}
}
else if(S2==16){
dh4(p,g_Sb,g_TF .h1[s],g_TF .h2[s],g_TF .h3[s],g_TF .q1[s],g_TF .q2[s],g_TF .q3[s],g_ciF,g_cisF);
}
else{
}
}
}

TGT static void outer_dif(u32*a,u32 h,const u32*t1,const u32*t2,const u32*t3,u32 ci,u32 cis){
 const __m256i p=_mm256_set1_epi32((int)NTP),imag=_mm256_set1_epi32(ci),imags=_mm256_set1_epi32(cis);
 for(u32 j=0;j<h;j+=8){u32 off=padpos(j),s=padpos(h);__m256i x0=ldu256(a+off),x1=ldu256(a+s+off),x2=ldu256(a+padpos(2*h)+off),x3=ldu256(a+padpos(3*h)+off);
  __m256i A=vaddm(x0,x2,p),B=vsubm(x0,x2,p),C=vaddm(x1,x3,p),D=vL(x1,x3,p),E=vS(D,imag,imags,p);
  stu256(a+off,vaddm(A,C,p));stu256(a+s+off,mpa(_mm256_add_epi32(B,E),ldu256(t1+j),p));stu256(a+padpos(2*h)+off,mpa(vL(A,C,p),ldu256(t2+j),p));stu256(a+padpos(3*h)+off,mpa(vL(B,E,p),ldu256(t3+j),p));
 }
}
TGT static void outer_dit(u32*a,u32 h,const u32*t1,const u32*t2,const u32*t3,u32 ci,u32 cis){
 const __m256i p=_mm256_set1_epi32((int)NTP),imag=_mm256_set1_epi32(ci),imags=_mm256_set1_epi32(cis);
 for(u32 j=0;j<h;j+=8){u32 off=padpos(j),s=padpos(h);__m256i x0=ldu256(a+off),x1=ldu256(a+s+off),x2=ldu256(a+padpos(2*h)+off),x3=ldu256(a+padpos(3*h)+off);
  x0=_mm256_min_epu32(x0,_mm256_sub_epi32(x0,p));x1=_mm256_min_epu32(x1,_mm256_sub_epi32(x1,p));x2=_mm256_min_epu32(x2,_mm256_sub_epi32(x2,p));x3=_mm256_min_epu32(x3,_mm256_sub_epi32(x3,p));
  x1=mpa(x1,ldu256(t1+j),p);x2=mpa(x2,ldu256(t2+j),p);x3=mpa(x3,ldu256(t3+j),p);
  __m256i A=vaddm(x0,x2,p),B=vsubm(x0,x2,p),C=vaddm(x1,x3,p),D=vL(x1,x3,p),E=vS(D,imag,imags,p);
  stu256(a+off,vaddm(A,C,p));stu256(a+s+off,vaddm(B,E,p));stu256(a+padpos(2*h)+off,vsubm(A,C,p));stu256(a+padpos(3*h)+off,vsubm(B,E,p));
 }
}
TGT static void ntt_fwd_tail(u32*a,u32 S,int lstart,int l0);
TGT static void ntt_fwd(u32*a,const u32*src=0,u32 len=0){
 (void)src;(void)len;for(int l=0;l<g_L;++l){u32 S=g_n>>(2*l),h=S/4;for(u32 b=0;b<g_n;b+=S)outer_dif(a+padpos(b),h,g_TF.l1[l],g_TF.l2[l],g_TF.l3[l],g_ciF,g_cisF);}
 for(u32 b=0;b<g_n;b+=g_Sb)fwd_sb_ring(a,b);
}
TGT static void ntt_fwd_pair_batch(u32**aa,u32**bb,int nb,const long*deltas){(void)deltas;g_pmdelta=0;for(int i=0;i<nb;++i){ntt_fwd(aa[i]);ntt_fwd(bb[i]);}}


#define I44 512u

TGT static void ntt_inv_tail(u32*a,const u32*Bmul){
for(u32 b=0;
b<g_n;
b+=g_Sb){
u32*p=a+padpos(b);
if(Bmul&&!E7Z3_PM_ON()){
const __m256i pv=_mm256_set1_epi32((int)NTP);
const u32*q=Bmul+padpos(b);
for(u32 i=0;
i<g_Sb;
i+=8){
__m256i x=ldu256((p+i));
__m256i y=ldu256((q+i));
stu256((p+i),vS(x,y,der_ws(y),pv));
}
}
const int h1fused=(g_nb>=2&&g_Sb>=32u&&(g_Sb%32u)==0u);
const int s64fused=(h1fused&&g_nb>=3&&g_Sb>=64u&&(g_Sb%64u)==0u);
for(int s=0;
s<g_nb;
s++){
u32 S=4u<<(2*s);
if(S>=32){
if(S==16384){
for(u32 i=0;
i<g_Sb;
i+=S){
if(g_Sb==16384&&!g_odd&&g_L==4&&(b&(3*g_Sb)))merge4_dit_h<4096,true>(p+i,16384,4096,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
else merge4_dit_h<4096>(p+i,16384,4096,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
}
else if(S==4096){
for(u32 i=0;
i<g_Sb;
i+=S){
if(g_Sb>=4*S&&(i&(3*S))!=0)merge4_dit_h<1024,true>(p+i,4096,1024,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
else merge4_dit_h<1024>(p+i,4096,1024,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
}
else if(S==1024){
for(u32 i=0;
i<g_Sb;
i+=S){
if(g_Sb>=4*S&&(i&(3*S))!=0)merge4_dit_h<256,true>(p+i,1024,256,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
else merge4_dit_h<256>(p+i,1024,256,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
}
else if(S==256){
for(u32 i=0;
i<g_Sb;
i+=S){
if(g_Sb>=4*S&&(i&(3*S))!=0)merge4_dit_h<64,true>(p+i,256,64,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
else merge4_dit_h<64>(p+i,256,64,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
}
else if(S==64){
if(!s64fused){
for(u32 i=0;
i<g_Sb;
i+=S)merge4_dit_h<16>(p+i,64,16,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
}
else{
for(u32 i=0;
i<g_Sb;
i+=S)merge4_dit_t(p+i,S,S>>2,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
}
else if(S==16){
if(s64fused)mhh64(p,g_Sb,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI,g_TI .h1[s+1],g_TI .h2[s+1],g_TI .h3[s+1],g_TI .q1[s+1],g_TI .q2[s+1],g_TI .q3[s+1]);
else if(h1fused)mhh(p,g_Sb,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
else merge4_dit_h4(p,g_Sb,g_TI .h1[s],g_TI .h2[s],g_TI .h3[s],g_TI .q1[s],g_TI .q2[s],g_TI .q3[s],g_ciI,g_cisI);
}
else{
if(!h1fused)mh1(p,g_Sb,g_ciI,g_cisI);
}
}
if(g_odd)merge2_dit(p,g_Sb,g_TI .b2,g_TI .b2s);

}

for(int l=g_L-1;l>=0;--l){u32 S=g_n>>(2*l),h=S/4;for(u32 b=0;b<g_n;b+=S)outer_dit(a+padpos(b),h,g_TI.l1[l],g_TI.l2[l],g_TI.l3[l],g_ciI,g_cisI);}
}

TGT static void ntt_inv_batch(u32**blocks,int count){g_pmdelta=0;for(int i=0;i<count;++i)ntt_inv_tail(blocks[i],0);}

struct AlignedBuffer {u32*p=nullptr;size_t capacity=0;~AlignedBuffer(){free(p);}void reserve(size_t words){if(words<=capacity)return;free(p);p=nullptr;if(posix_memalign((void**)&p,64,words*4))std::abort();capacity=words;}};
static AlignedBuffer tableF,tableI,paddedA,paddedB;
static u32 active_n=0,active_p=0;
static u32 primitive_root(u32 p){u32 x=p-1;std::vector<u32>f;for(u32 d=2;u64(d)*d<=x;++d)if(x%d==0){f.push_back(d);while(x%d==0)x/=d;}if(x>1)f.push_back(x);for(u32 g=2;;++g){bool ok=true;for(u32 q:f)if(powmod(g,(p-1)/q)==1){ok=false;break;}if(ok)return g;}}
static void initialize(u32 n,u32 p){assert(n>=64&&(n&(n-1))==0&&p<(1u<<31)&&(p-1)%n==0);g_mod=p;if(active_n==n&&active_p==p)return;g_shmagic=u32(~0ull/u64(p));g_shint=u32((~0ull/u64(p))>>32);if(active_p!=p)g_generator=primitive_root(p);tableF.reserve(3ull*n+262144);tableI.reserve(3ull*n+262144);g_tabF=tableF.p;g_tabI=tableI.p;ntt_init(n);g_pmdelta=0;_mm_sfence();active_n=n;active_p=p;}
static void load_padded(u32*dst,const u32*src,u32 n){for(u32 i=0;i<n;i+=16384u)std::memcpy(dst+padpos(i),src+i,size_t(std::min(16384u,n-i))*4);}
TGT static void store_padded(u32*dst,const u32*src,u32 n,u32 scale){const __m256i p=_mm256_set1_epi32((int)NTP),s=_mm256_set1_epi32(scale),ss=_mm256_set1_epi32(shoup_c(scale));for(u32 base=0;base<n;base+=16384u){const u32*x=src+padpos(base);u32 count=std::min(16384u,n-base);for(u32 i=0;i<count;i+=8)stu256(dst+base+i,vS(ldu256(x+i),s,ss,p));}}
// In-place cyclic product of Montgomery-represented coefficients. Output is
// Montgomery-represented too, compatible with the batch canonical NTT backend.
TGT static bool cyclic_impl(u32*a,u32*b,u32 n,u32 p,bool mont){
 if(!n||(n&(n-1))||p>=(1u<<31)||(p-1)%n)return false;
 if(n<64){std::vector<u32>c(n);u32 inv=1;for(int i=0;i<5;++i)inv*=2-p*inv;u32 small_pinv=0-inv;for(u32 i=0;i<n;++i)for(u32 j=0;j<n;++j){u64 z=u64(a[i])*b[j],m=u32(z)*u64(small_pinv);u32 v=mont?u32((z+u64(u32(m))*p)>>32):u32(z%p);if(v>=p)v-=p;u32&vout=c[(i+j)&(n-1)];vout+=v;if(vout>=p)vout-=p;}std::copy(c.begin(),c.end(),a);return true;}
 initialize(n,p);paddedA.reserve(padpos(n)+64);paddedB.reserve(padpos(n)+64);load_padded(paddedA.p,a,n);load_padded(paddedB.p,b,n);u32*aa[1]={paddedA.p},*bb[1]={paddedB.p};ntt_fwd_pair_batch(aa,bb,1,nullptr);
 const __m256i pv=_mm256_set1_epi32(p);for(u32 base=0;base<n;base+=16384u){u32*x=paddedA.p+padpos(base),*y=paddedB.p+padpos(base);u32 count=std::min(16384u,n-base);for(u32 i=0;i<count;i+=8){__m256i ax=ldu256(x+i),by=ldu256(y+i);ax=_mm256_min_epu32(ax,_mm256_sub_epi32(ax,pv));by=_mm256_min_epu32(by,_mm256_sub_epi32(by,pv));stu256(x+i,vmont(ax,by,pv));}}
 ntt_inv_batch(aa,1);u32 scale=powmod(n,p-2);if(!mont)scale=mulmod(scale,RRW);store_padded(a,paddedA.p,n,scale);return true;
}

TGT static bool cyclic_mont(u32*a,u32*b,u32 n,u32 p){return cyclic_impl(a,b,n,p,true);}
TGT static bool cyclic_normal(u32*a,u32*b,u32 n,u32 p){return cyclic_impl(a,b,n,p,false);}
static u32 pair_mask(){u32 mask=0;for(int k=0;k<g_nb;++k)mask|=1u<<(2*k);for(int k=0;k<g_L;++k)mask|=1u<<(g_log2n-2-2*k);return mask;}
static u32 reorder_index(u32 i,u32 lo){u32 hi=lo<<1;return ((i&lo)<<1)|((i&hi)>>1)|(i&~(lo|hi));}
TGT static void binary_layout(u32*dst,const u32*src,u32 n,bool to_padded){
 const __m256i p=_mm256_set1_epi32(NTP),perm=_mm256_setr_epi32(0,2,1,3,4,6,5,7);u32 mask=pair_mask();
 for(u32 i=0;i<n;i+=16){u32 offset=padpos(reorder_index(i,mask));const u32*in=to_padded?src+i:src+offset;__m256i a=ldu256(in),b=ldu256(in+8);a=_mm256_min_epu32(a,_mm256_sub_epi32(a,p));b=_mm256_min_epu32(b,_mm256_sub_epi32(b,p));a=_mm256_permutevar8x32_epi32(a,perm);b=_mm256_permutevar8x32_epi32(b,perm);u32*out=to_padded?dst+offset:dst+i;stu256(out,_mm256_permute2x128_si256(a,b,0x20));stu256(out+8,_mm256_permute2x128_si256(a,b,0x31));}
}
static u32 bitreverse(u32 x,unsigned log){u32 y=0;for(unsigned i=0;i<log;++i){y=y*2+(x&1);x>>=1;}return y;}
static void tiny_transform(u32*a,u32 n,u32 p,bool inverse,bool normal){u32 oldmod=g_mod;g_mod=p;u32 root=powmod(primitive_root(p),(p-1)/n);if(inverse)root=powmod(root,p-2);std::vector<u32>out(n);unsigned log=__builtin_ctz(n);for(u32 k=0;k<n;++k){u32 step=powmod(root,k),w=1,sum=0;for(u32 j=0;j<n;++j){u32 value=a[inverse?bitreverse(j,log):j];sum=(sum+u64(value)*w)%p;w=mulmod(w,step);}if(inverse&&normal)sum=mulmod(sum,powmod(n,p-2));out[inverse?k:bitreverse(k,log)]=sum;}std::copy(out.begin(),out.end(),a);g_mod=oldmod;}
// Public residue-array transforms use BINARY bit-reversed frequency slots,
// exactly matching the canonical Nyaan interface, irrespective of raw radix4
// engine layout and the position of its odd radix2 stage.
TGT static bool ntt_forward(u32*a,u32 n,u32 p){if(!n||(n&(n-1))||p>=1u<<31||(p-1)%n)return false;if(n<64){tiny_transform(a,n,p,false,false);return true;}initialize(n,p);paddedA.reserve(padpos(n)+64);load_padded(paddedA.p,a,n);ntt_fwd(paddedA.p);binary_layout(a,paddedA.p,n,false);return true;}
TGT static bool ntt_inverse(u32*a,u32 n,u32 p,bool normalize=true){if(!n||(n&(n-1))||p>=1u<<31||(p-1)%n)return false;if(n<64){tiny_transform(a,n,p,true,normalize);return true;}initialize(n,p);paddedA.reserve(padpos(n)+64);binary_layout(paddedA.p,a,n,true);u32*blocks[1]={paddedA.p};ntt_inv_batch(blocks,1);store_padded(a,paddedA.p,n,normalize?powmod(n,p-2):1);return true;}

// The pruned inverse needs a different base length. Keep its complete table
// context separately, preserving the large forward cache across node calls.
struct SavedContext {u32 mod=0,generator=0,shmagic=0,shint=0,n=0,Sb=0,pinv=0,ciF=0,cisF=0,ciI=0,cisI=0,wf=0,wif=0,active_n=0,active_p=0;int log2n=0,L=0,odd=0,nb=0,BT=13;long pmdelta=0;Tabs TF{},TI{};AlignedBuffer F,I;};
static SavedContext prefix_context;
static void swap_context(SavedContext&s){
 std::swap(s.mod,g_mod);std::swap(s.generator,g_generator);std::swap(s.shmagic,g_shmagic);std::swap(s.shint,g_shint);std::swap(s.n,g_n);std::swap(s.Sb,g_Sb);std::swap(s.pinv,g_pinv);std::swap(s.ciF,g_ciF);std::swap(s.cisF,g_cisF);std::swap(s.ciI,g_ciI);std::swap(s.cisI,g_cisI);std::swap(s.wf,g_wf);std::swap(s.wif,g_wif);std::swap(s.active_n,active_n);std::swap(s.active_p,active_p);std::swap(s.log2n,g_log2n);std::swap(s.L,g_L);std::swap(s.odd,g_odd);std::swap(s.nb,g_nb);std::swap(s.BT,g_BT);std::swap(s.pmdelta,g_pmdelta);std::swap(s.TF,g_TF);std::swap(s.TI,g_TI);std::swap(s.F.p,tableF.p);std::swap(s.F.capacity,tableF.capacity);std::swap(s.I.p,tableI.p);std::swap(s.I.capacity,tableI.capacity);g_tabF=tableF.p;g_tabI=tableI.p;
}
// Compute only the first wanted time coefficients from binary-bit-reversed
// frequency input. Small subtransforms use the optimized complete inverse;
// upper DIT stages retain only the requested prefixes. Input is destroyed.
TGT static bool ntt_inverse_prefix(u32*a,u32 n,u32 wanted,u32 p,bool normalize=true){
 if(!n||(n&(n-1))||!wanted||wanted>n||p>=1u<<31||(p-1)%n)return false;
 u32 base=1;while(base<wanted)base*=2;base=std::min(n,std::max(base,64u));
 if(base==n)return ntt_inverse(a,n,p,normalize);
 swap_context(prefix_context);initialize(base,p);paddedA.reserve(padpos(base)+64);u32*blocks[1]={paddedA.p};
 for(u32 offset=0;offset<n;offset+=base){binary_layout(paddedA.p,a+offset,base,true);ntt_inv_batch(blocks,1);store_padded(a+offset,paddedA.p,base,1);}
 u32 count=(wanted+7)&~7u;std::vector<u32>tw(count),ts(count);const __m256i pv=_mm256_set1_epi32(p);
 for(u32 step=base;step<n;step*=2){u32 root=powmod(g_generator,(p-1)/(step*2));root=powmod(root,p-2);u32 v=1;for(u32 i=0;i<count;++i){tw[i]=v;ts[i]=shoup_c(v);v=mulmod(v,root);}
  for(u32 offset=0;offset<n;offset+=step*2){u32*x=a+offset,*y=x+step;for(u32 i=0;i<count;i+=8){__m256i left=ldu256(x+i),right=vS(ldu256(y+i),ldu256(tw.data()+i),ldu256(ts.data()+i),pv);__m256i sum=_mm256_add_epi32(left,right);stu256(x+i,_mm256_min_epu32(sum,_mm256_sub_epi32(sum,pv)));}}
 }
 if(normalize){u32 scale=powmod(n,p-2);const __m256i sv=_mm256_set1_epi32(scale),ss=_mm256_set1_epi32(shoup_c(scale));for(u32 i=0;i<count;i+=8)stu256(a+i,vS(ldu256(a+i),sv,ss,pv));}
 swap_context(prefix_context);return true;
}
}
#pragma GCC pop_options

#undef NTT_MAXN
#undef PF_DIST
#undef NTP
#undef TPAD
#undef TOPPAD
#undef PADJUMP
#undef PADJUMP_S
#undef TGT
#undef Q_T5OFF
#undef Q_SCF
#undef Q_SCI
#undef RRW
#undef MLO_V
#undef E7Z3_PM_ON
#undef H4B32
#undef E7_TILE
#undef E7Z_SBFUSE
#undef E7Z_IVFUSE
#undef I44

#include <memory>
#include <array>
#pragma GCC push_options
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
namespace mixed3_field {
using u32=uint32_t;using u64=uint64_t;
inline u32 mul(u32 a,u32 b,u32 p){return u64(a)*b%p;}
inline u32 power(u32 a,u32 e,u32 p){u32 z=1;for(;e;e>>=1,a=mul(a,a,p))if(e&1)z=mul(z,a,p);return z;}
inline u32 shoup(u32 a,u32 p){return (u64(a)<<32)/p;}
inline u32 smul(u32 a,u32 b,u32 bs,u32 p){u32 z=u32(u64(a)*b)-u32((u64(a)*bs)>>32)*p;return z>=p?z-p:z;}
inline u32 add(u32 a,u32 b,u32 p){u32 z=a+b;return z>=p?z-p:z;}
inline u32 sub(u32 a,u32 b,u32 p){return a>=b?a-b:a+p-b;}
inline __m256i av(__m256i a,__m256i b,__m256i p){__m256i z=_mm256_add_epi32(a,b);return _mm256_min_epu32(z,_mm256_sub_epi32(z,p));}
inline __m256i sv(__m256i a,__m256i b,__m256i p){__m256i z=_mm256_sub_epi32(a,b);return _mm256_min_epu32(z,_mm256_add_epi32(z,p));}
inline u32 primitive_root(u32 p){u32 x=p-1;std::vector<u32>f;for(u32 d=2;u64(d)*d<=x;++d)if(x%d==0){f.push_back(d);while(x%d==0)x/=d;}if(x>1)f.push_back(x);for(u32 g=2;;++g){bool ok=true;for(u32 q:f)if(power(g,(p-1)/q,p)==1){ok=false;break;}if(ok)return g;}}
inline __m256i derive_shoup(__m256i w,u32 mod){u64 reciprocal=~0ull/u64(mod);__m256i q=_mm256_add_epi32(_mm256_mullo_epi32(w,_mm256_set1_epi32(u32(reciprocal>>32))),fast_e7_field::vmulhi(w,_mm256_set1_epi32(u32(reciprocal))));__m256i q1=_mm256_add_epi32(q,_mm256_set1_epi32(1)),pv=_mm256_set1_epi32(mod),high=fast_e7_field::vmulhi(q1,pv),low=_mm256_mullo_epi32(q1,pv),zero=_mm256_setzero_si256();__m256i correction=_mm256_or_si256(_mm256_xor_si256(_mm256_cmpeq_epi32(high,w),_mm256_set1_epi32(-1)),_mm256_cmpeq_epi32(low,zero));return _mm256_add_epi32(q,_mm256_and_si256(correction,_mm256_set1_epi32(1)));}
inline u32 mont_pinv(u32 p){u32 x=1;for(int i=0;i<5;++i)x*=2-p*x;return 0-x;}
inline __m256i mont_product(__m256i a,__m256i b,u32 mod,u32 pinv){__m256i np=_mm256_set1_epi64x(pinv),p64=_mm256_set1_epi64x(mod),ae=_mm256_mul_epu32(a,b),ao=_mm256_mul_epu32(_mm256_bsrli_epi128(a,4),_mm256_bsrli_epi128(b,4)),me=_mm256_mul_epu32(ae,np),mo=_mm256_mul_epu32(ao,np);ae=_mm256_add_epi64(ae,_mm256_mul_epu32(me,p64));ao=_mm256_add_epi64(ao,_mm256_mul_epu32(mo,p64));auto value=_mm256_or_si256(_mm256_bsrli_epi128(ae,4),ao),pv=_mm256_set1_epi32(mod);return _mm256_min_epu32(value,_mm256_sub_epi32(value,pv));}
struct Twiddles{std::vector<u32>w1,w1s,w2,w2s;};
class Plan {
 u32 n,m,p,gen,zeta,zetai,pinv;std::vector<u32>temp;Twiddles fwd,inv;
 void select(){}
 void build(Twiddles&t,u32 root){t.w1.resize(m);t.w2.resize(m);t.w1s.resize(m);t.w2s.resize(m);t.w1[0]=t.w2[0]=1;const __m256i pp=_mm256_set1_epi32(p);
  for(u32 block=1;block<m;block*=2){u32 factor=power(root,m/(2*block),p),factor2=mul(factor,factor,p),s1=shoup(factor,p),s2=shoup(factor2,p),i=0;auto v1=_mm256_set1_epi32(factor),v2=_mm256_set1_epi32(factor2),vs1=_mm256_set1_epi32(s1),vs2=_mm256_set1_epi32(s2);
   for(;i+8<=block;i+=8){_mm256_storeu_si256((__m256i*)(t.w1.data()+block+i),fast_e7_field::vS(_mm256_loadu_si256((const __m256i*)(t.w1.data()+i)),v1,vs1,pp));_mm256_storeu_si256((__m256i*)(t.w2.data()+block+i),fast_e7_field::vS(_mm256_loadu_si256((const __m256i*)(t.w2.data()+i)),v2,vs2,pp));}
   for(;i<block;++i){t.w1[block+i]=smul(t.w1[i],factor,s1,p);t.w2[block+i]=smul(t.w2[i],factor2,s2,p);}
  }
  u32 i=0;for(;i+8<=m;i+=8){_mm256_storeu_si256((__m256i*)(t.w1s.data()+i),derive_shoup(_mm256_loadu_si256((const __m256i*)(t.w1.data()+i)),p));_mm256_storeu_si256((__m256i*)(t.w2s.data()+i),derive_shoup(_mm256_loadu_si256((const __m256i*)(t.w2.data()+i)),p));}for(;i<m;++i){t.w1s[i]=shoup(t.w1[i],p);t.w2s[i]=shoup(t.w2[i],p);}
 }
 void combine(const u32*a,u32*dst,bool inverse){const auto&t=inverse?inv:fwd;const __m256i pp=_mm256_set1_epi32(p),zz=_mm256_set1_epi32(inverse?zetai:zeta),zs=_mm256_set1_epi32(shoup(inverse?zetai:zeta,p));u32 i=0;
  for(;i+8<=m;i+=8){__m256i x0=_mm256_loadu_si256((const __m256i*)(a+i)),x1=_mm256_loadu_si256((const __m256i*)(a+m+i)),x2=_mm256_loadu_si256((const __m256i*)(a+2*m+i));
   if(!inverse){x1=fast_e7_field::vS(x1,_mm256_loadu_si256((const __m256i*)(t.w1.data()+i)),_mm256_loadu_si256((const __m256i*)(t.w1s.data()+i)),pp);x2=fast_e7_field::vS(x2,_mm256_loadu_si256((const __m256i*)(t.w2.data()+i)),_mm256_loadu_si256((const __m256i*)(t.w2s.data()+i)),pp);}
   auto delta=fast_e7_field::vS(sv(x1,x2,pp),zz,zs,pp),y0=av(x0,av(x1,x2,pp),pp),y1=av(sv(x0,x2,pp),delta,pp),y2=sv(sv(x0,x1,pp),delta,pp);
   if(inverse){y1=fast_e7_field::vS(y1,_mm256_loadu_si256((const __m256i*)(t.w1.data()+i)),_mm256_loadu_si256((const __m256i*)(t.w1s.data()+i)),pp);y2=fast_e7_field::vS(y2,_mm256_loadu_si256((const __m256i*)(t.w2.data()+i)),_mm256_loadu_si256((const __m256i*)(t.w2s.data()+i)),pp);}
   _mm256_storeu_si256((__m256i*)(dst+i),y0);_mm256_storeu_si256((__m256i*)(dst+m+i),y1);_mm256_storeu_si256((__m256i*)(dst+2*m+i),y2);
  }
  for(;i<m;++i){u32 x0=a[i],x1=a[m+i],x2=a[2*m+i];if(!inverse){x1=smul(x1,t.w1[i],t.w1s[i],p);x2=smul(x2,t.w2[i],t.w2s[i],p);}u32 zz0=inverse?zetai:zeta,delta=smul(sub(x1,x2,p),zz0,shoup(zz0,p),p),y0=add(x0,add(x1,x2,p),p),y1=add(sub(x0,x2,p),delta,p),y2=sub(sub(x0,x1,p),delta,p);if(inverse){y1=smul(y1,t.w1[i],t.w1s[i],p);y2=smul(y2,t.w2[i],t.w2s[i],p);}dst[i]=y0;dst[m+i]=y1;dst[2*m+i]=y2;}
 }
public:
 Plan(u32 length,u32 mod):n(length),m(length/3),p(mod),gen(0),temp(length){assert(n>=3&&n%3==0&&!(m&(m-1))&&p<(1u<<31)&&(p-1)%n==0);select();gen=primitive_root(p);pinv=mont_pinv(p);u32 root=power(gen,(p-1)/n,p);zeta=power(root,m,p);zetai=mul(zeta,zeta,p);build(fwd,root);build(inv,power(root,p-2,p));}
 u32 length()const{return n;}u32 modulus()const{return p;}
 bool forward(u32*a){select();for(u32 i=0;i<m;++i){temp[i]=a[3*i];temp[m+i]=a[3*i+1];temp[2*m+i]=a[3*i+2];}for(unsigned r=0;r<3;++r)if(!fast_e7_field::ntt_forward(temp.data()+r*m,m,p))return false;combine(temp.data(),a,false);return true;}
 bool inverse_prefix(u32*a,u32 wanted,bool normalize=true,u32 extra_scale=1){if(!wanted||wanted>n)return false;select();combine(a,temp.data(),true);for(u32 r=0;r<3&&r<wanted;++r){u32 count=(wanted+2-r)/3;if(!(count>=m/128?fast_e7_field::ntt_inverse(temp.data()+r*m,m,p,normalize):fast_e7_field::ntt_inverse_prefix(temp.data()+r*m,m,count,p,normalize)))return false;}u32 scale=normalize?power(3,p-2,p):1;scale=mul(scale,extra_scale,p);u32 ss=shoup(scale,p);for(u32 i=0;i<wanted;++i)a[i]=smul(temp[(i%3)*m+i/3],scale,ss,p);return true;}
 void pointwise(u32*a,const u32*b){select();const __m256i pp=_mm256_set1_epi32(p);u32 i=0;for(;i+8<=n;i+=8)_mm256_storeu_si256((__m256i*)(a+i),mont_product(_mm256_loadu_si256((const __m256i*)(a+i)),_mm256_loadu_si256((const __m256i*)(b+i)),p,pinv));for(;i<n;++i){u64 z=u64(a[i])*b[i];u32 q=u32(z)*pinv;u32 v=(z+u64(q)*p)>>32;a[i]=v>=p?v-p:v;}}
 bool cyclic(u32*a,u32*b,bool mont=true){if(!forward(a)||(a!=b&&!forward(b)))return false;pointwise(a,b);return inverse_prefix(a,n,true,mont?1:u32((1ull<<32)%p));}
};
static std::unique_ptr<Plan>active;
inline Plan&select(u32 n,u32 p){if(!active||active->length()!=n||active->modulus()!=p)active=std::make_unique<Plan>(n,p);return *active;}
inline bool mixed3_forward(u32*a,u32 n,u32 p){return select(n,p).forward(a);}
inline bool mixed3_inverse_prefix(u32*a,u32 n,u32 wanted,u32 p,bool normalize=true){return select(n,p).inverse_prefix(a,wanted,normalize);}
inline bool cyclic_mont(u32*a,u32*b,u32 n,u32 p){return select(n,p).cyclic(a,b,true);}
inline bool cyclic_normal(u32*a,u32*b,u32 n,u32 p){return select(n,p).cyclic(a,b,false);}
}
#pragma GCC pop_options

#pragma GCC push_options
#pragma GCC optimize("O3")
#pragma GCC target("avx2")
namespace good3_field {
using mixed3_field::u32;using mixed3_field::u64;using mixed3_field::power;using mixed3_field::mul;using mixed3_field::shoup;using mixed3_field::smul;using mixed3_field::add;using mixed3_field::sub;using mixed3_field::av;using mixed3_field::sv;
class Plan {
 u32 n,m,p,gen,zeta,zetai,pinv;std::vector<u32>temp;
 void select(){}
 void combine(const u32*a,u32*dst,bool inverse){const __m256i pp=_mm256_set1_epi32(p),zz=_mm256_set1_epi32(inverse?zetai:zeta),zs=_mm256_set1_epi32(shoup(inverse?zetai:zeta,p));u32 i=0;for(;i+8<=m;i+=8){auto x0=_mm256_loadu_si256((const __m256i*)(a+i)),x1=_mm256_loadu_si256((const __m256i*)(a+m+i)),x2=_mm256_loadu_si256((const __m256i*)(a+2*m+i)),delta=fast_e7_field::vS(sv(x1,x2,pp),zz,zs,pp);_mm256_storeu_si256((__m256i*)(dst+i),av(x0,av(x1,x2,pp),pp));_mm256_storeu_si256((__m256i*)(dst+m+i),av(sv(x0,x2,pp),delta,pp));_mm256_storeu_si256((__m256i*)(dst+2*m+i),sv(sv(x0,x1,pp),delta,pp));}for(;i<m;++i){u32 x0=a[i],x1=a[m+i],x2=a[2*m+i],z=inverse?zetai:zeta,d=smul(sub(x1,x2,p),z,shoup(z,p),p);dst[i]=add(x0,add(x1,x2,p),p);dst[m+i]=add(sub(x0,x2,p),d,p);dst[2*m+i]=sub(sub(x0,x1,p),d,p);}}
public:
 Plan(u32 length,u32 mod):n(length),m(length/3),p(mod),gen(0),temp(length){assert(n>=3&&n%3==0&&!(m&(m-1))&&p<(1u<<31)&&(p-1)%n==0);select();gen=mixed3_field::primitive_root(p);pinv=mixed3_field::mont_pinv(p);zeta=power(gen,(p-1)/3,p);zetai=mul(zeta,zeta,p);}
 u32 length()const{return n;}u32 modulus()const{return p;}
 bool forward(u32*a){select();unsigned invm=(m%3)==1?1:2;u32 q=0;unsigned which[3][3];for(unsigned r=0;r<3;++r)for(unsigned s=0;s<3;++s)which[r][s]=((r+3-s)*invm)%3;
  for(;q+3<=m;q+=3)for(unsigned r=0;r<3;++r){temp[r*m+q]=a[which[r][0]*m+q];temp[r*m+q+1]=a[which[r][1]*m+q+1];temp[r*m+q+2]=a[which[r][2]*m+q+2];}
  for(;q<m;++q)for(unsigned r=0;r<3;++r)temp[r*m+q]=a[which[r][q%3]*m+q];
  for(unsigned r=0;r<3;++r)if(!fast_e7_field::ntt_forward(temp.data()+r*m,m,p))return false;combine(temp.data(),a,false);return true;}
 bool inverse_prefix(u32*a,u32 wanted,bool normalize=true,u32 extra_scale=1){if(!wanted||wanted>n)return false;select();combine(a,temp.data(),true);u32 count=std::min(wanted,m);for(unsigned r=0;r<3;++r)if(!(count>=m/128?fast_e7_field::ntt_inverse(temp.data()+r*m,m,p,normalize):fast_e7_field::ntt_inverse_prefix(temp.data()+r*m,m,count,p,normalize)))return false;u32 scale=mul(normalize?power(3,p-2,p):1,extra_scale,p),ss=shoup(scale,p),mask=m-1,i=0;for(;i+3<=wanted;i+=3){a[i]=smul(temp[i&mask],scale,ss,p);a[i+1]=smul(temp[m+((i+1)&mask)],scale,ss,p);a[i+2]=smul(temp[2*m+((i+2)&mask)],scale,ss,p);}for(;i<wanted;++i)a[i]=smul(temp[(i%3)*m+(i&mask)],scale,ss,p);return true;}
 void pointwise(u32*a,const u32*b){select();auto pp=_mm256_set1_epi32(p);u32 i=0;for(;i+8<=n;i+=8)_mm256_storeu_si256((__m256i*)(a+i),mixed3_field::mont_product(_mm256_loadu_si256((const __m256i*)(a+i)),_mm256_loadu_si256((const __m256i*)(b+i)),p,pinv));for(;i<n;++i){u64 z=u64(a[i])*b[i];u32 q=u32(z)*pinv;u32 v=(z+u64(q)*p)>>32;a[i]=v>=p?v-p:v;}}
 bool cyclic(u32*a,u32*b,bool mont=true){if(!forward(a)||(a!=b&&!forward(b)))return false;pointwise(a,b);return inverse_prefix(a,n,true,mont?1:u32((1ull<<32)%p));}
};
static std::unique_ptr<Plan>active;
inline Plan&select(u32 n,u32 p){if(!active||active->length()!=n||active->modulus()!=p)active=std::make_unique<Plan>(n,p);return *active;}
inline bool mixed3_forward(u32*a,u32 n,u32 p){return select(n,p).forward(a);}
inline bool mixed3_inverse_prefix(u32*a,u32 n,u32 wanted,u32 p,bool normalize=true){return select(n,p).inverse_prefix(a,wanted,normalize);}
inline bool cyclic_mont(u32*a,u32*b,u32 n,u32 p){return select(n,p).cyclic(a,b,true);}
inline bool cyclic_normal(u32*a,u32*b,u32 n,u32 p){return select(n,p).cyclic(a,b,false);}
}
#pragma GCC pop_options

#ifndef DUCK_SHA256_HPP
#define DUCK_SHA256_HPP
#include <stddef.h>
#include <stdint.h>
#include <string.h>
// Portable SHA-256; input is the exact byte sequence, with no native-word hashing.
namespace duck_sha256 {
static inline uint32_t rotate(uint32_t x,unsigned n) {return (x>>n)|(x<<(32-n));}
static void compress(uint32_t h[8],const unsigned char block[64]) {
    static const uint32_t k[64]={
      0x428a2f98u,0x71374491u,0xb5c0fbcfu,0xe9b5dba5u,0x3956c25bu,0x59f111f1u,0x923f82a4u,0xab1c5ed5u,
      0xd807aa98u,0x12835b01u,0x243185beu,0x550c7dc3u,0x72be5d74u,0x80deb1feu,0x9bdc06a7u,0xc19bf174u,
      0xe49b69c1u,0xefbe4786u,0x0fc19dc6u,0x240ca1ccu,0x2de92c6fu,0x4a7484aau,0x5cb0a9dcu,0x76f988dau,
      0x983e5152u,0xa831c66du,0xb00327c8u,0xbf597fc7u,0xc6e00bf3u,0xd5a79147u,0x06ca6351u,0x14292967u,
      0x27b70a85u,0x2e1b2138u,0x4d2c6dfcu,0x53380d13u,0x650a7354u,0x766a0abbu,0x81c2c92eu,0x92722c85u,
      0xa2bfe8a1u,0xa81a664bu,0xc24b8b70u,0xc76c51a3u,0xd192e819u,0xd6990624u,0xf40e3585u,0x106aa070u,
      0x19a4c116u,0x1e376c08u,0x2748774cu,0x34b0bcb5u,0x391c0cb3u,0x4ed8aa4au,0x5b9cca4fu,0x682e6ff3u,
      0x748f82eeu,0x78a5636fu,0x84c87814u,0x8cc70208u,0x90befffau,0xa4506cebu,0xbef9a3f7u,0xc67178f2u};
    uint32_t w[64];
    for(unsigned i=0;i<16;++i)w[i]=(uint32_t(block[4*i])<<24)|(uint32_t(block[4*i+1])<<16)|(uint32_t(block[4*i+2])<<8)|block[4*i+3];
    for(unsigned i=16;i<64;++i){uint32_t a=w[i-15],b=w[i-2];w[i]=w[i-16]+(rotate(a,7)^rotate(a,18)^(a>>3))+w[i-7]+(rotate(b,17)^rotate(b,19)^(b>>10));}
    uint32_t a=h[0],b=h[1],c=h[2],d=h[3],e=h[4],f=h[5],g=h[6],v=h[7];
    for(unsigned i=0;i<64;++i){
        uint32_t t1=v+(rotate(e,6)^rotate(e,11)^rotate(e,25))+((e&f)^(~e&g))+k[i]+w[i];
        uint32_t t2=(rotate(a,2)^rotate(a,13)^rotate(a,22))+((a&b)^(a&c)^(b&c));
        v=g;g=f;f=e;e=d+t1;d=c;c=b;b=a;a=t1+t2;
    }
    h[0]+=a;h[1]+=b;h[2]+=c;h[3]+=d;h[4]+=e;h[5]+=f;h[6]+=g;h[7]+=v;
}
inline void digest(const void*data,size_t length,unsigned char out[32]) {
    uint32_t h[8]={0x6a09e667u,0xbb67ae85u,0x3c6ef372u,0xa54ff53au,0x510e527fu,0x9b05688cu,0x1f83d9abu,0x5be0cd19u};
    const unsigned char*p=(const unsigned char*)data;
    size_t whole=length/64;
    for(size_t i=0;i<whole;++i)compress(h,p+64*i);
    unsigned char tail[128]={};size_t remain=length%64;
    if(remain)memcpy(tail,p+64*whole,remain);
    tail[remain]=0x80;size_t padded=remain<56?64:128;
    uint64_t bits=uint64_t(length)*8;
    for(unsigned i=0;i<8;++i)tail[padded-1-i]=(unsigned char)(bits>>(8*i));
    compress(h,tail);if(padded==128)compress(h,tail+64);
    for(unsigned i=0;i<8;++i)for(unsigned j=0;j<4;++j)out[4*i+j]=(unsigned char)(h[i]>>(24-8*j));
}
}
#endif


#if BENCH_GOOD
using FieldPlan=good3_field::Plan;
#else
using FieldPlan=mixed3_field::Plan;
#endif
struct Rng{uint64_t s;unsigned next(){uint64_t z=(s+=0x9e3779b97f4a7c15ULL);z=(z^(z>>30))*0xbf58476d1ce4e5b9ULL;z=(z^(z>>27))*0x94d049bb133111ebULL;return unsigned(z^(z>>31));}};
static uint64_t tick(){_mm_lfence();return __rdtsc();}
__attribute__((target("avx2"))) int main(){auto all=tick();constexpr unsigned p=BENCH_P,M=1u<<22,n=3*M,big=1u<<24,orig=10000000;const unsigned R=(1ull<<32)%p;std::vector<unsigned>a(n),b(n),x(big),y(big),weights(orig),as,ys;fast_e7_field::paddedA.reserve(fast_e7_field::padpos(big)+64);memset(fast_e7_field::paddedA.p,0,fast_e7_field::paddedA.capacity*4);Rng ar{0x1357214865430012ULL},br{0x2143658730024589ULL};for(unsigned i=0;i<orig;++i){unsigned value=uint64_t(ar.next()%p)*R%p;if(i==orig-1)a[0]=x[0]=value;else{a[n-orig+1+i]=value;x[big-orig+1+i]=value;}weights[i]=uint64_t(br.next()%p)*R%p;}auto t=tick();fast_e7_field::initialize(M,p);FieldPlan plan(n,p);auto setup3=tick()-t;t=tick();plan.forward(a.data());as.resize(n);for(unsigned i=0;i<n;i+=8){auto v=mixed3_field::mont_product(_mm256_loadu_si256((const __m256i*)(a.data()+i)),_mm256_set1_epi32(1),p,mixed3_field::mont_pinv(p));_mm256_storeu_si256((__m256i*)(a.data()+i),v);_mm256_storeu_si256((__m256i*)(as.data()+i),mixed3_field::derive_shoup(v,p));}auto af3=tick()-t;fast_e7_field::SavedContext other;fast_e7_field::swap_context(other);t=tick();fast_e7_field::initialize(big,p);auto setup2=tick()-t;t=tick();if(!fast_e7_field::ntt_forward(x.data(),big,p))__builtin_trap();ys.resize(big);for(unsigned i=0;i<big;i+=8){auto v=mixed3_field::mont_product(_mm256_loadu_si256((const __m256i*)(x.data()+i)),_mm256_set1_epi32(1),p,mixed3_field::mont_pinv(p));_mm256_storeu_si256((__m256i*)(x.data()+i),v);_mm256_storeu_si256((__m256i*)(ys.data()+i),mixed3_field::derive_shoup(v,p));}auto af2=tick()-t;fast_e7_field::swap_context(other);uint64_t tm[4][2],tp[4][2];unsigned wants[4]={16384,131072,1048576,2500000};
 for(unsigned wi=0;wi<4;++wi){unsigned wanted=wants[wi],r=2*wanted;for(unsigned rep=0;rep<2;++rep){t=tick();std::fill(b.begin(),b.end(),0);if(wanted<=65536){for(unsigned lo=wanted;lo<orig;lo+=32768){unsigned hi=std::min(lo+32768,orig);for(unsigned j=wanted;j<std::min(r,hi);++j){unsigned value=weights[j];for(unsigned d=((lo+j-1)/j)*j;d<hi;d+=j){unsigned&v=b[orig-1-d];v+=value;if(v>=p)v-=p;}}}}else{for(unsigned q=1;q<=(orig-1)/wanted;++q){unsigned je=std::min(r-1,(orig-1)/q);for(unsigned j=wanted;j<=je;++j){unsigned&v=b[orig-1-q*j];v+=weights[j];if(v>=p)v-=p;}}}if(!plan.forward(b.data()))__builtin_trap();for(unsigned i=0;i<n;i+=8)_mm256_storeu_si256((__m256i*)(b.data()+i),fast_e7_field::vS(_mm256_loadu_si256((const __m256i*)(b.data()+i)),_mm256_loadu_si256((const __m256i*)(a.data()+i)),_mm256_loadu_si256((const __m256i*)(as.data()+i)),_mm256_set1_epi32(p)));if(!plan.inverse_prefix(b.data(),wanted))__builtin_trap();tm[wi][rep]=tick()-t;}fast_e7_field::swap_context(other);for(unsigned rep=0;rep<2;++rep){t=tick();std::fill(y.begin(),y.end(),0);if(wanted<=65536){for(unsigned lo=wanted;lo<orig;lo+=32768){unsigned hi=std::min(lo+32768,orig);for(unsigned j=wanted;j<std::min(r,hi);++j){unsigned value=weights[j];for(unsigned d=((lo+j-1)/j)*j;d<hi;d+=j){unsigned&v=y[orig-1-d];v+=value;if(v>=p)v-=p;}}}}else{for(unsigned q=1;q<=(orig-1)/wanted;++q){unsigned je=std::min(r-1,(orig-1)/q);for(unsigned j=wanted;j<=je;++j){unsigned&v=y[orig-1-q*j];v+=weights[j];if(v>=p)v-=p;}}}if(!fast_e7_field::ntt_forward(y.data(),big,p))__builtin_trap();for(unsigned i=0;i<big;i+=8)_mm256_storeu_si256((__m256i*)(y.data()+i),fast_e7_field::vS(_mm256_loadu_si256((const __m256i*)(y.data()+i)),_mm256_loadu_si256((const __m256i*)(x.data()+i)),_mm256_loadu_si256((const __m256i*)(ys.data()+i)),_mm256_set1_epi32(p)));if(!(wanted>=big/128?fast_e7_field::ntt_inverse(y.data(),big,p,true):fast_e7_field::ntt_inverse_prefix(y.data(),big,wanted,p,true)))__builtin_trap();tp[wi][rep]=tick()-t;}if(!std::equal(b.begin(),b.begin()+wanted,y.begin()))__builtin_trap();fast_e7_field::swap_context(other);}
 unsigned char hash[32];duck_sha256::digest(b.data(),size_t(wants[3])*4,hash);auto total=tick()-all;printf("S3%llu S2%llu A3%llu A2%llu T%llu ",(unsigned long long)setup3,(unsigned long long)setup2,(unsigned long long)af3,(unsigned long long)af2,(unsigned long long)total);for(unsigned i=0;i<4;++i)printf("%u:%llu/%llu ",wants[i],(unsigned long long)tm[i][1],(unsigned long long)tp[i][1]);puts("ok1");}

评测结果