#pragma GCC optimize("O3")
// Blocked radix-4 DIF/DIT NTT mod 2013265921 using Shoup modular multiplication.
// Per stage: 6 tables (t1,t2,t3 = twiddle values; s1,s2,s3 = Shoup constants), contiguous.
// M180ARM:c_r2PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
// DRAM levels: block sizes n, n/4, ...; then in-cache base blocks of size ~2^15..2^19.
// DIF forward: natural in -> digit-reversed out. DIT inverse: digit-reversed in -> natural out.
#pragma once
#include <immintrin.h>
#include <cstdint>
#include <cstdint>
#include <cstring>
#include <cstdlib>
// AVX-256 split-load fix: g++-9 -O2 with NO -march lowers _mm256_loadu_si256 into
// vmovdqu128 + vinserti128 (2 loads + a shuffle). One asm instruction emits the real load.
__attribute__((target("avx2"))) static inline __m256i ldu256(const void *p) {
__m256i r; __asm__("vmovdqu %1, %0" : "=x"(r) : "m"(*(const __m256i *)p)); return r; }
// AVX-256 split-STORE fix, symmetric to ldu256 above: g++-9 -O2 with NO -march lowers
// stu256(2 uops, one of them on p5).
// One asm instruction emits the real 32-byte store, 1 uop, and -- unlike the aligned
// intrinsic _mm256_store_si256, which compiles to vmovdqa with NO alignment proof and
// therefore FAULTS at runtime on a misaligned pointer -- vmovdqu carries no precondition.
// The memory is declared as an OUTPUT operand ("=m"), so gcc knows the location is written
// and no full "memory" clobber (which would block optimisation) is needed.
__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 256
static inline void pf5(const void *p) { _mm_prefetch((const char *)p, _MM_HINT_T0); }
#define NTP 2013265921u
#define TGT __attribute__((target("avx2")))
#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;
// ARM TABN: measured need is 4 210 686 u32 for the radix-4 head (n=2^22) + M+margin for the
// single materialised r5 series (init_r5 aliases series 1..3 onto series 0). The old 5*NTT_MAXN/2
// reservation is 20x the used footprint and is what put .bss at 839 MB -- over custom_test's cap.
#define TABN (4300000u + (1u << 22) + 1024u)
static u32 g_tabF[TABN];
static u32 g_tabI[TABN];
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; }
static inline u32 shoup1s(u32 a, u32 w, u32 ws) {
u32 t = (u32)((u64)a * w);
u32 q = (u32)(((u64)a * ws) >> 32);
u32 r = t - (u32)((u64)q * NTP);
return r >= NTP ? r - NTP : r;
}
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_srli_epi64(a, 32), _mm256_srli_epi64(b, 32));
e = _mm256_srli_epi64(e, 32);
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));
}
// M = floor(2^64/p) = 2^33 + MLO ; ws ~= 2*w + mulhi(w, MLO) (error <= 1, valid for Shoup)
#define MLO_V 572662301u
TGT static inline __m256i der_ws(__m256i w) {
__m256i mlo = _mm256_set1_epi32((int)MLO_V);
return _mm256_add_epi32(_mm256_add_epi32(w, w), vmulhi(w, mlo));
}
TGT static inline __m256i vshoup(__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_srli_epi64(a, 32), _mm_srli_epi64(b, 32));
e = _mm_srli_epi64(e, 32);
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);
}
// fused level-0 DIF stage: reads coefficients directly from src ([0,len] valid, zero beyond).
// Requires len >= h and len < 2h+... (x2 = x3 = 0 assumed).
TGT static void stage0_dif_src(u32 *a, const u32 *src, u32 len, u32 h,
const u32 *t1, const u32 *t2, const u32 *t3, 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);
const __m256i zero = _mm256_setzero_si256();
u32 j = 0;
u32 lim1 = len - h;
for (; j + 8 <= lim1 + 1; j += 8) {
__builtin_prefetch(src + j + 256);
__builtin_prefetch(src + j + h + 256);
__builtin_prefetch(a + j + 256);
__builtin_prefetch(a + j + h + 256);
__builtin_prefetch(a + j + 2 * h + 256);
__builtin_prefetch(a + j + 3 * h + 256);
__builtin_prefetch(t1 + j + 256);
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + h));
__m256i u0 = vaddm(x0, x1, pv), u2 = vsubm(x0, x1, pv);
__m256i E = vshoup(x1, cI, cIs, pv);
__m256i u1 = vaddm(x0, E, pv), u3 = vsubm(x0, E, pv);
{ __m256i w = ldu256((t1 + j)); u1 = vshoup(u1, w, der_ws(w), pv); }
{ __m256i w = ldu256((t2 + j)); u2 = vshoup(u2, w, der_ws(w), pv); }
{ __m256i w = ldu256((t3 + j)); u3 = vshoup(u3, w, der_ws(w), pv); }
stu256((a + j), u0);
stu256((a + j + h), u1);
stu256((a + j + 2 * h), u2);
stu256((a + j + 3 * h), u3);
}
for (; j <= lim1 && j < h; j++) {
u32 x0 = (j <= len) ? src[j] : 0, x1 = (j + h <= len) ? src[j + h] : 0;
u32 u0 = (u32)(((u64)x0 + x1) % NTP);
u32 u2 = (u32)(((u64)x0 + NTP - x1) % NTP);
u32 E = mulmod(x1, ci);
u32 u1 = (u32)(((u64)x0 + E) % NTP), u3 = (u32)(((u64)x0 + NTP - E) % NTP);
u1 = shoup1s(u1, t1[j], shoup_c(t1[j]));
u2 = shoup1s(u2, t2[j], shoup_c(t2[j]));
u3 = shoup1s(u3, t3[j], shoup_c(t3[j]));
a[j] = u0; a[j + h] = u1; a[j + 2 * h] = u2; a[j + 3 * h] = u3;
}
for (; j + 8 <= h; j += 8) {
__m256i x0 = (j + 8 <= len + 1) ? ldu256((src + j)) : zero;
__m256i u0 = x0, u1 = x0, u2 = x0, u3 = x0;
{ __m256i w = ldu256((t1 + j)); u1 = vshoup(u1, w, der_ws(w), pv); }
{ __m256i w = ldu256((t2 + j)); u2 = vshoup(u2, w, der_ws(w), pv); }
{ __m256i w = ldu256((t3 + j)); u3 = vshoup(u3, w, der_ws(w), pv); }
stu256((a + j), u0);
stu256((a + j + h), u1);
stu256((a + j + 2 * h), u2);
stu256((a + j + 3 * h), u3);
}
for (; j < h; j++) {
u32 x0 = (j <= len) ? src[j] : 0;
u32 u1 = shoup1s(x0, t1[j], shoup_c(t1[j]));
u32 u2 = shoup1s(x0, t2[j], shoup_c(t2[j]));
u32 u3 = shoup1s(x0, t3[j], shoup_c(t3[j]));
a[j] = x0; a[j + h] = u1; a[j + 2 * h] = u2; a[j + 3 * h] = u3;
}
}
// ================= forward (DIF) =================
TGT static void stage4_dif(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) {
__builtin_prefetch(p + j + 256);
__builtin_prefetch(p + j + h + 256);
__builtin_prefetch(p + j + 2 * h + 256);
__builtin_prefetch(p + j + 3 * h + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__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 = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i A2 = vaddm(y0, y2, pv), B2 = vsubm(y0, y2, pv);
__m256i C2 = vaddm(y1, y3, pv), D2 = vsubm(y1, y3, pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i v0 = vaddm(A2, C2, pv), v2 = vsubm(A2, C2, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i F = vshoup(D2, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
__m256i v1 = vaddm(B2, F, pv), v3 = vsubm(B2, F, pv);
{ __m256i w = ldu256((t1 + j)); u1 = vshoup(u1, w, der_ws(w), pv); }
{ __m256i w = ldu256((t2 + j)); u2 = vshoup(u2, w, der_ws(w), pv); }
{ __m256i w = ldu256((t3 + j)); u3 = vshoup(u3, w, der_ws(w), pv); }
{ __m256i w = ldu256((t1 + j + 8)); v1 = vshoup(v1, w, der_ws(w), pv); }
{ __m256i w = ldu256((t2 + j + 8)); v2 = vshoup(v2, w, der_ws(w), pv); }
{ __m256i w = ldu256((t3 + j + 8)); v3 = vshoup(v3, w, der_ws(w), 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 + 256);
__builtin_prefetch(p + j + h + 256);
__builtin_prefetch(p + j + 2 * h + 256);
__builtin_prefetch(p + j + 3 * h + 256);
__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 = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
{ __m256i w = ldu256((t1 + j)); u1 = vshoup(u1, w, der_ws(w), pv); }
{ __m256i w = ldu256((t2 + j)); u2 = vshoup(u2, w, der_ws(w), pv); }
{ __m256i w = ldu256((t3 + j)); u3 = vshoup(u3, w, der_ws(w), pv); }
stu256((p + j), u0);
stu256((p + j + h), u1);
stu256((p + j + 2 * h), u2);
stu256((p + j + 3 * h), u3);
}
}
}
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 = vsubm(x1, x3, pv);
__m256i A2 = vaddm(y0, y2, pv), B2 = vsubm(y0, y2, pv);
__m256i C2 = vaddm(y1, y3, pv), D2 = vsubm(y1, y3, pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i v0 = vaddm(A2, C2, pv), v2 = vsubm(A2, C2, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i F = vshoup(D2, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
__m256i v1 = vaddm(B2, F, pv), v3 = vsubm(B2, F, pv);
{ __m256i w = ldu256((t1 + j)); __m256i ws = ldu256((s1 + j)); u1 = vshoup(u1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j)); __m256i ws = ldu256((s2 + j)); u2 = vshoup(u2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j)); __m256i ws = ldu256((s3 + j)); u3 = vshoup(u3, w, ws, pv); }
{ __m256i w = ldu256((t1 + j + 8)); __m256i ws = ldu256((s1 + j + 8)); v1 = vshoup(v1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j + 8)); __m256i ws = ldu256((s2 + j + 8)); v2 = vshoup(v2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j + 8)); __m256i ws = ldu256((s3 + j + 8)); v3 = vshoup(v3, w, ws, 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 + 256);
__builtin_prefetch(p + j + h + 256);
__builtin_prefetch(p + j + 2 * h + 256);
__builtin_prefetch(p + j + 3 * h + 256);
__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 = vsubm(x1, x3, pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
{ __m256i w = ldu256((t1 + j)); __m256i ws = ldu256((s1 + j)); u1 = vshoup(u1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j)); __m256i ws = ldu256((s2 + j)); u2 = vshoup(u2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j)); __m256i ws = ldu256((s3 + j)); u3 = vshoup(u3, w, ws, 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 *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 = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i A2 = vaddm(y0, y2, pv), B2 = vsubm(y0, y2, pv);
__m256i C2 = vaddm(y1, y3, pv), D2 = vsubm(y1, y3, pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i v0 = vaddm(A2, C2, pv), v2 = vsubm(A2, C2, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i F = vshoup(D2, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
__m256i v1 = vaddm(B2, F, pv), v3 = vsubm(B2, F, pv);
{ __m256i w = ldu256((t1 + j)); __m256i ws = ldu256((s1 + j)); u1 = vshoup(u1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j)); __m256i ws = ldu256((s2 + j)); u2 = vshoup(u2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j)); __m256i ws = ldu256((s3 + j)); u3 = vshoup(u3, w, ws, pv); }
{ __m256i w = ldu256((t1 + j + 8)); __m256i ws = ldu256((s1 + j + 8)); v1 = vshoup(v1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j + 8)); __m256i ws = ldu256((s2 + j + 8)); v2 = vshoup(v2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j + 8)); __m256i ws = ldu256((s3 + j + 8)); v3 = vshoup(v3, w, ws, 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 + 256);
__builtin_prefetch(p + j + H + 256);
__builtin_prefetch(p + j + 2 * H + 256);
__builtin_prefetch(p + j + 3 * H + 256);
__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 = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
{ __m256i w = ldu256((t1 + j)); __m256i ws = ldu256((s1 + j)); u1 = vshoup(u1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j)); __m256i ws = ldu256((s2 + j)); u2 = vshoup(u2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j)); __m256i ws = ldu256((s3 + j)); u3 = vshoup(u3, w, ws, pv); }
stu256((p + j), u0);
stu256((p + j + H), u1);
stu256((p + j + 2 * H), u2);
stu256((p + j + 3 * H), u3);
}
}
}
TGT static void stage4_dif_nt(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 = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i A2 = vaddm(y0, y2, pv), B2 = vsubm(y0, y2, pv);
__m256i C2 = vaddm(y1, y3, pv), D2 = vsubm(y1, y3, pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i v0 = vaddm(A2, C2, pv), v2 = vsubm(A2, C2, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i F = vshoup(D2, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
__m256i v1 = vaddm(B2, F, pv), v3 = vsubm(B2, F, pv);
u1 = vshoup(u1, ldu256((t1 + j)), der_ws(ldu256((t1 + j))), pv);
u2 = vshoup(u2, ldu256((t2 + j)), der_ws(ldu256((t2 + j))), pv);
u3 = vshoup(u3, ldu256((t3 + j)), der_ws(ldu256((t3 + j))), pv);
v1 = vshoup(v1, ldu256((t1 + j + 8)), der_ws(ldu256((t1 + j + 8))), pv);
v2 = vshoup(v2, ldu256((t2 + j + 8)), der_ws(ldu256((t2 + j + 8))), pv);
v3 = vshoup(v3, ldu256((t3 + j + 8)), der_ws(ldu256((t3 + j + 8))), pv);
_mm256_stream_si256((__m256i *)(p + j), u0);
_mm256_stream_si256((__m256i *)(p + j + h), u1);
_mm256_stream_si256((__m256i *)(p + j + 2 * h), u2);
_mm256_stream_si256((__m256i *)(p + j + 3 * h), u3);
_mm256_stream_si256((__m256i *)(p + j + 8), v0);
_mm256_stream_si256((__m256i *)(p + j + h + 8), v1);
_mm256_stream_si256((__m256i *)(p + j + 2 * h + 8), v2);
_mm256_stream_si256((__m256i *)(p + j + 3 * h + 8), v3);
}
for (; j + 8 <= 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));
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
u1 = vshoup(u1, ldu256((t1 + j)), der_ws(ldu256((t1 + j))), pv);
u2 = vshoup(u2, ldu256((t2 + j)), der_ws(ldu256((t2 + j))), pv);
u3 = vshoup(u3, ldu256((t3 + j)), der_ws(ldu256((t3 + j))), pv);
_mm256_stream_si256((__m256i *)(p + j), u0);
_mm256_stream_si256((__m256i *)(p + j + h), u1);
_mm256_stream_si256((__m256i *)(p + j + 2 * h), u2);
_mm256_stream_si256((__m256i *)(p + j + 3 * h), u3);
}
}
}
// ===== ARM R8: radix-8 first rung. Fuse the last TWO radix-4 DIF stages (S2=16 then S2=4)
// ===== into one pass: the S2=16 butterfly's four outputs are transformed IN REGISTERS by the
// ===== S2=4 butterfly (twiddle = 1) before they are stored, so the intermediate never round-trips.
TGT static inline __m256i h1_256(__m256i x, __m256i cI, __m256i cIs, __m256i pv) {
__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 = vshoup(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);
return _mm256_blend_epi32(sum, dif, 0xCC);
}
TGT static inline __m128i h1_128(__m128i x, __m128i cI, __m128i cIs, __m128i pv) {
__m128i t = _mm_shuffle_epi32(x, _MM_SHUFFLE(1, 0, 3, 2));
__m128i s = add4(x, t, pv);
__m128i d = sub4(x, t, pv);
__m128i Z = _mm_shuffle_epi32(d, 0x00);
__m128i P = _mm_blend_epi32(s, Z, 0xAA);
__m128i E = sh4(d, cI, cIs, pv);
__m128i E2 = _mm_shuffle_epi32(E, 0x55);
__m128i W = _mm_shuffle_epi32(s, 0x55);
__m128i Q = _mm_blend_epi32(W, E2, 0xAA);
__m128i sum = add4(P, Q, pv);
__m128i dif = sub4(P, Q, pv);
return _mm_blend_epi32(sum, dif, 0x0C);
}
TGT static void stage4_dif_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), X1 = _mm_loadu_si128((const __m128i *)s1);
__m128i W2 = _mm_loadu_si128((const __m128i *)t2), X2 = _mm_loadu_si128((const __m128i *)s2);
__m128i W3 = _mm_loadu_si128((const __m128i *)t3), X3 = _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 = _mm_add_epi32(sub4(x1, x3, pv4), pv4);
__m128i u0 = add4(A, C, pv4), u2 = sub4(A, C, pv4);
__m128i E = sh4(D, cI4, cIs4, pv4);
__m128i u1 = add4(B, E, pv4), u3 = sub4(B, E, pv4);
u1 = sh4(u1, W1, X1, pv4); u2 = sh4(u2, W2, X2, pv4); u3 = sh4(u3, W3, X3, pv4);
_mm_storeu_si128((__m128i *)(p), h1_128(u0, cI4, cIs4, pv4));
_mm_storeu_si128((__m128i *)(p + 4), h1_128(u1, cI4, cIs4, pv4));
_mm_storeu_si128((__m128i *)(p + 8), h1_128(u2, cI4, cIs4, pv4));
_mm_storeu_si128((__m128i *)(p + 12), h1_128(u3, cI4, cIs4, pv4));
}
return;
}
__m256i W1 = bc128(_mm_loadu_si128((const __m128i *)t1));
__m256i W2 = bc128(_mm_loadu_si128((const __m128i *)t2));
__m256i W3 = bc128(_mm_loadu_si128((const __m128i *)t3));
__m256i X1 = bc128(_mm_loadu_si128((const __m128i *)s1));
__m256i X2 = bc128(_mm_loadu_si128((const __m128i *)s2));
__m256i X3 = 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));
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i u0 = vaddm(A, C, pv), u2 = vsubm(A, C, pv);
__m256i E = vshoup(D, cI, cIs, pv);
__m256i u1 = vaddm(B, E, pv), u3 = vsubm(B, E, pv);
u1 = vshoup(u1, W1, X1, pv); u2 = vshoup(u2, W2, X2, pv); u3 = vshoup(u3, W3, X3, pv);
st2((__m128i *)(p), (__m128i *)(p + 16), h1_256(u0, cI, cIs, pv));
st2((__m128i *)(p + 4), (__m128i *)(p + 20), h1_256(u1, cI, cIs, pv));
st2((__m128i *)(p + 8), (__m128i *)(p + 24), h1_256(u2, cI, cIs, pv));
st2((__m128i *)(p + 12), (__m128i *)(p + 28), h1_256(u3, cI, cIs, pv));
}
}
TGT static void stage4_dif_h1(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 = vshoup(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 = vshoup(vsubm(x0, x1, pv), ldu256((t + j)),
ldu256((s + j)), pv);
stu256((a + j), u0);
stu256((a + j + h), u1);
}
}
// ================= inverse (DIT) =================
TGT static void merge4_dit(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) {
__builtin_prefetch(p + j + 256);
__builtin_prefetch(p + j + h + 256);
__builtin_prefetch(p + j + 2 * h + 256);
__builtin_prefetch(p + j + 3 * h + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__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 w = ldu256((t1 + j)); x1 = vshoup(x1, w, der_ws(w), pv); }
{ __m256i w = ldu256((t2 + j)); x2 = vshoup(x2, w, der_ws(w), pv); }
{ __m256i w = ldu256((t3 + j)); x3 = vshoup(x3, w, der_ws(w), pv); }
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i E = vshoup(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));
}
}
}
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));
{ __m256i w = ldu256((t1 + j)); __m256i ws = ldu256((s1 + j)); x1 = vshoup(x1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j)); __m256i ws = ldu256((s2 + j)); x2 = vshoup(x2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j)); __m256i ws = ldu256((s3 + j)); x3 = vshoup(x3, w, ws, pv); }
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = vsubm(x1, x3, pv);
__m256i E = vshoup(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>
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 __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));
{ __m256i w = ldu256((t1 + j)); __m256i ws = ldu256((s1 + j)); x1 = vshoup(x1, w, ws, pv); }
{ __m256i w = ldu256((t2 + j)); __m256i ws = ldu256((s2 + j)); x2 = vshoup(x2, w, ws, pv); }
{ __m256i w = ldu256((t3 + j)); __m256i ws = ldu256((s3 + j)); x3 = vshoup(x3, w, ws, pv); }
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i E = vshoup(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));
}
}
}
TGT static void merge4_dit_nt(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 = vshoup(x1, ldu256((t1 + j)), der_ws(ldu256((t1 + j))), pv);
x2 = vshoup(x2, ldu256((t2 + j)), der_ws(ldu256((t2 + j))), pv);
x3 = vshoup(x3, ldu256((t3 + j)), der_ws(ldu256((t3 + j))), pv);
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i E = vshoup(D, cI, cIs, pv);
_mm256_stream_si256((__m256i *)(p + j), vaddm(A, C, pv));
_mm256_stream_si256((__m256i *)(p + j + h), vaddm(B, E, pv));
_mm256_stream_si256((__m256i *)(p + j + 2 * h), vsubm(A, C, pv));
_mm256_stream_si256((__m256i *)(p + j + 3 * h), 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), X1 = _mm_loadu_si128((const __m128i *)s1);
__m128i W2 = _mm_loadu_si128((const __m128i *)t2), X2 = _mm_loadu_si128((const __m128i *)s2);
__m128i W3 = _mm_loadu_si128((const __m128i *)t3), X3 = _mm_loadu_si128((const __m128i *)s3);
for (u32 i = 0; i < nsub; i += 16) {
u32 *p = a + i;
__m128i x0 = h1_128(_mm_loadu_si128((const __m128i *)(p)), cI4, cIs4, pv4);
__m128i x1 = sh4(h1_128(_mm_loadu_si128((const __m128i *)(p + 4)), cI4, cIs4, pv4), W1, X1, pv4);
__m128i x2 = sh4(h1_128(_mm_loadu_si128((const __m128i *)(p + 8)), cI4, cIs4, pv4), W2, X2, pv4);
__m128i x3 = sh4(h1_128(_mm_loadu_si128((const __m128i *)(p + 12)), cI4, cIs4, pv4), W3, X3, pv4);
__m128i A = add4(x0, x2, pv4), B = sub4(x0, x2, pv4);
__m128i C = add4(x1, x3, pv4), D = _mm_add_epi32(sub4(x1, x3, pv4), 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 W2 = bc128(_mm_loadu_si128((const __m128i *)t2));
__m256i W3 = bc128(_mm_loadu_si128((const __m128i *)t3));
__m256i X1 = bc128(_mm_loadu_si128((const __m128i *)s1));
__m256i X2 = bc128(_mm_loadu_si128((const __m128i *)s2));
__m256i X3 = bc128(_mm_loadu_si128((const __m128i *)s3));
for (u32 i = 0; i < nsub; i += 32) {
u32 *p = a + i;
__m256i x0 = h1_256(ld2((const __m128i *)(p), (const __m128i *)(p + 16)), cI, cIs, pv);
__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 = vshoup(h1_256(x1, cI, cIs, pv), W1, X1, pv);
x2 = vshoup(h1_256(x2, cI, cIs, pv), W2, X2, pv);
x3 = vshoup(h1_256(x3, cI, cIs, pv), W3, X3, pv);
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i E = vshoup(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));
}
}
// ===== LANE sched_1002e7 PTFUSE: the pointwise pass FOLDED INTO ITS ONLY CONSUMER.
// ntt_inv ran `if (Bmul) for (i<g_Sb) A[i] = A[i]*B[i]` and then, at s=1, called
// merge4_dit_h4(p, g_Sb, ...) over EXACTLY the same range (nsub = g_Sb = 16384 at this
// row's shape). Folding the multiply into the kernel's loads is arithmetically
// identical: the butterfly is linear and the pointwise is elementwise, and every
// element is loaded exactly once per s=1 iteration before its butterfly. Net:
// -252 MB of traffic reduced to -84 MB (B still read, at the same offsets) and
// -~65 M instructions (the pass was ~3.5 instr/element; the fold adds 0.375).
TGT static void merge4_dit_h4pt(u32 *a, u32 nsub, const u32 *bq, 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) {
for (u32 i = 0; i < nsub; i++) a[i] = (u32)((u64)a[i] * bq[i] % NTP); /* tiny case: the pointwise here */
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), X1 = _mm_loadu_si128((const __m128i *)s1);
__m128i W2 = _mm_loadu_si128((const __m128i *)t2), X2 = _mm_loadu_si128((const __m128i *)s2);
__m128i W3 = _mm_loadu_si128((const __m128i *)t3), X3 = _mm_loadu_si128((const __m128i *)s3);
for (u32 i = 0; i < nsub; i += 16) {
u32 *p = a + i;
__m128i x0 = h1_128(_mm_loadu_si128((const __m128i *)(p)), cI4, cIs4, pv4);
__m128i x1 = sh4(h1_128(_mm_loadu_si128((const __m128i *)(p + 4)), cI4, cIs4, pv4), W1, X1, pv4);
__m128i x2 = sh4(h1_128(_mm_loadu_si128((const __m128i *)(p + 8)), cI4, cIs4, pv4), W2, X2, pv4);
__m128i x3 = sh4(h1_128(_mm_loadu_si128((const __m128i *)(p + 12)), cI4, cIs4, pv4), W3, X3, pv4);
__m128i A = add4(x0, x2, pv4), B = sub4(x0, x2, pv4);
__m128i C = add4(x1, x3, pv4), D = _mm_add_epi32(sub4(x1, x3, pv4), 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 W2 = bc128(_mm_loadu_si128((const __m128i *)t2));
__m256i W3 = bc128(_mm_loadu_si128((const __m128i *)t3));
__m256i X1 = bc128(_mm_loadu_si128((const __m128i *)s1));
__m256i X2 = bc128(_mm_loadu_si128((const __m128i *)s2));
__m256i X3 = bc128(_mm_loadu_si128((const __m128i *)s3));
for (u32 i = 0; i < nsub; i += 32) {
u32 *p = a + i;
const u32 *bp = bq + i;
__m256i x0 = h1_256(vmont(ld2((const __m128i *)(p), (const __m128i *)(p + 16)),
ld2((const __m128i *)(bp), (const __m128i *)(bp + 16)), pv), cI, cIs, pv);
__m256i x1 = vmont(ld2((const __m128i *)(p + 4), (const __m128i *)(p + 20)),
ld2((const __m128i *)(bp + 4), (const __m128i *)(bp + 20)), pv);
__m256i x2 = vmont(ld2((const __m128i *)(p + 8), (const __m128i *)(p + 24)),
ld2((const __m128i *)(bp + 8), (const __m128i *)(bp + 24)), pv);
__m256i x3 = vmont(ld2((const __m128i *)(p + 12), (const __m128i *)(p + 28)),
ld2((const __m128i *)(bp + 12), (const __m128i *)(bp + 28)), pv);
x1 = vshoup(h1_256(x1, cI, cIs, pv), W1, X1, pv);
x2 = vshoup(h1_256(x2, cI, cIs, pv), W2, X2, pv);
x3 = vshoup(h1_256(x3, cI, cIs, pv), W3, X3, pv);
__m256i A = vaddm(x0, x2, pv), B = vsubm(x0, x2, pv);
__m256i C = vaddm(x1, x3, pv), D = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), pv);
__m256i E = vshoup(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 void merge4_dit_h1(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 = vshoup(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 = vshoup(ldu256((a + j + h)),
ldu256((t + j)),
ldu256((s + j)), pv);
stu256((a + j), vaddm(x0, x1, pv));
stu256((a + j + h), vsubm(x0, x1, pv));
}
}
// ================= table generation =================
// dst[i] = v0*r^i ; sc[i] = floor(dst[i]*2^32/p) (within 1, sufficient for Shoup)
TGT static void gen_geo(u32 *dst, u32 *sc, u32 cnt, u32 v0, u32 r) {
if (!cnt) return;
const __m256i pv = _mm256_set1_epi32((int)NTP);
const u32 Mlo = 572662301u;
__m256i mlo = _mm256_set1_epi32((int)Mlo);
u32 seq[8], cur = v0 % NTP;
for (int i = 0; i < 8; i++) { seq[i] = cur; cur = mulmod(cur, r); }
u32 r8 = 1; for (int t = 0; t < 8; t++) r8 = mulmod(r8, r);
u32 rs8 = shoup_c(r8);
__m256i v = ldu256(seq);
__m256i r8v = _mm256_set1_epi32((int)r8), rs8v = _mm256_set1_epi32((int)rs8);
u32 i = 0;
for (; i + 8 <= cnt; i += 8) {
stu256((dst + i), v);
__m256i two = _mm256_add_epi32(v, v);
stu256((sc + i), _mm256_add_epi32(two, vmulhi(v, mlo)));
v = vshoup(v, r8v, rs8v, pv);
}
if (i < cnt) {
u32 x = v0 % NTP;
for (u32 t = 0; t < i; t++) x = mulmod(x, r);
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;
const __m256i pv = _mm256_set1_epi32((int)NTP);
u32 seq[8], cur = v0 % NTP;
for (int i = 0; i < 8; i++) { seq[i] = cur; cur = mulmod(cur, r); }
u32 r8 = 1; for (int t = 0; t < 8; t++) r8 = mulmod(r8, r);
u32 rs8 = shoup_c(r8);
__m256i v = ldu256(seq);
__m256i r8v = _mm256_set1_epi32((int)r8), rs8v = _mm256_set1_epi32((int)rs8);
u32 i = 0;
for (; i + 8 <= cnt; i += 8) { stu256((dst + i), v); v = vshoup(v, r8v, rs8v, pv); }
if (i < cnt) { u32 x = v0 % NTP; for (u32 t = 0; t < i; t++) x = mulmod(x, r); 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 *p = tab;
u32 S = n;
for (int l = 0; l < L; l++) {
u32 M = S >> 2;
T.l1[l] = p; T.l2[l] = p + M; T.l3[l] = p + 2 * M;
T.m1[l] = p; T.m2[l] = p + M; T.m3[l] = p + 2 * M;
u32 base = powmod(root, n / S);
u32 b1 = base, b2 = mulmod(base, base), b3 = mulmod(b2, base);
gen_geo_val(p, M, 1, b1);
gen_geo_val(p + M, M, 1, b2);
gen_geo_val(p + 2 * M, M, 1, b3);
p += 3 * M;
S >>= 2;
}
T.b2 = 0; T.b2s = 0;
if (odd) {
u32 M = Sb >> 1;
T.b2 = p; T.b2s = p + M;
gen_geo(p, p + M, M, 1, powmod(root, n / Sb));
p += 2 * M;
}
int s = 0;
u32 S2v = 4;
while (s < g_nb) {
u32 M = S2v >> 2;
T.h1[s] = p; T.h2[s] = p + M; T.h3[s] = p + 2 * M;
T.q1[s] = p + 3 * M; T.q2[s] = p + 4 * M; T.q3[s] = p + 5 * M;
u32 base = powmod(root, n / S2v);
u32 b1 = base, b2 = mulmod(base, base), b3 = mulmod(b2, base);
gen_geo(p, p + 3 * M, M, 1, b1);
gen_geo(p + M, p + 4 * M, M, 1, b2);
gen_geo(p + 2 * M, p + 5 * M, M, 1, b3);
p += 6 * M;
S2v <<= 2;
s++;
}
}
static u32 g_tabwords;
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(31u, (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);
{ u32 S=n; u32 tot=0; for (int l=0;l<L;l++){tot+=3*(S>>2);S>>=2;} u32 S2v=4; for(int s=0;s<g_nb;s++){tot+=6*(S2v>>2);S2v<<=2;} if(g_odd) tot+=2*(g_Sb>>1); g_tabwords=tot; }
build_tabs(g_TF, g_tabF, n, g_Sb, L, g_odd, w);
build_tabs(g_TI, g_tabI, n, g_Sb, L, g_odd, wi);
}
TGT static void ntt_fwd(u32 *a, const u32 *src = 0, u32 len = 0) {
u32 S = g_n;
int l0 = 0;
if (src && len >= (g_n >> 2) && len <= (g_n >> 1)) {
stage0_dif_src(a, src, len, g_n >> 2, g_TF.l1[0], g_TF.l2[0], g_TF.l3[0], g_ciF, g_cisF);
S = g_n >> 2;
l0 = 1;
}
for (int l = l0; l < g_L; l++) {
u32 h = S >> 2;
for (u32 b = 0; b < g_n; b += S)
stage4_dif(a + b, S, h, g_TF.l1[l], g_TF.l2[l], g_TF.l3[l], g_TF.m1[l], g_TF.m2[l], g_TF.m3[l], g_ciF, g_cisF);
_mm_sfence();
S = h;
}
for (u32 b = 0; b < g_n; b += g_Sb) {
u32 *p = a + 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) {
for (u32 i = 0; i < g_Sb; i += S2)
{ u32 *pp = p + i;
switch (S2) {
case 16384: stage4_dif_h<4096>(pp, 16384, 4096, 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); break;
case 4096: stage4_dif_h<1024>(pp, 4096, 1024, 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); break;
case 1024: stage4_dif_h<256>(pp, 1024, 256, 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); break;
case 256: stage4_dif_h<64>(pp, 256, 64, 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); break;
case 64: stage4_dif_h<16>(pp, 64, 16, 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); break;
default: stage4_dif_t(pp, 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) {
stage4_dif_h4(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);
} // S2=4 FUSED into the S2=16 pass above (ARM R8)
}
}
}
TGT static void ntt_inv(u32 *a, const u32 *Bmul = 0) {
for (u32 b = 0; b < g_n; b += g_Sb) {
u32 *p = a + b;
const u32 *qpt = Bmul ? (Bmul + b) : 0; /* PTFUSE: the pointwise's B base for this block */
if (0) {
const __m256i pv = _mm256_set1_epi32((int)NTP);
const u32 R2v = (u32)(((u64)1 << 32) % NTP * (((u64)1 << 32) % NTP) % NTP);
__m256i R2 = _mm256_set1_epi32((int)R2v);
const u32 *q = Bmul + b;
for (u32 i = 0; i < g_Sb; i += 8) {
__m256i x = ldu256((p + i));
__m256i y = ldu256((q + i));
stu256((p + i), vmont(x, y, pv));
}
}
for (int s = 0; s < g_nb; s++) {
u32 S = 4u << (2 * s);
if (S >= 32) {
for (u32 i = 0; i < g_Sb; i += S)
{ u32 *pp = p + i;
switch (S) {
case 16384: merge4_dit_h<4096>(pp, 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); break;
case 4096: merge4_dit_h<1024>(pp, 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); break;
case 1024: merge4_dit_h<256>(pp, 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); break;
case 256: merge4_dit_h<64>(pp, 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); break;
case 64: merge4_dit_h<16>(pp, 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); break;
default: merge4_dit_t(pp, 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 (Bmul) merge4_dit_h4pt(p, g_Sb, qpt, 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);
} // S2=4 FUSED into the S2=16 pass below (ARM R8: DIT order is h1 THEN h4)
}
if (g_odd) merge2_dit(p, g_Sb, g_TI.b2, g_TI.b2s);
}
for (int l = g_L - 1; l >= 0; l--) {
u32 SB = g_n >> (2 * l);
u32 h = SB >> 2;
for (u32 b = 0; b < g_n; b += SB)
merge4_dit(a + b, SB, h, g_TI.l1[l], g_TI.l2[l], g_TI.l3[l], g_TI.m1[l], g_TI.m2[l], g_TI.m3[l], g_ciI, g_cisI);
_mm_sfence();
}
}
// ============================================================================
// Mixed-radix radix-5 top level for N = 5 * 2^22 (same prime p = 2013265921).
// stage: for each u < M: Y_q = sum_t x[tM+u] rho^{tq} (rho = w^M, order 5)
// a[qM+u] = Y_q * w^{q u}
// then 5 size-M transforms. The composite equals the true DFT up to a
// digit-reversed permutation (harmless for pointwise multiplication), exactly
// as for the radix-3 level (verified numerically over F_p on a small N).
// DFT-5 with a=rho+rho^4, b=rho^2+rho^3, c=rho-rho^4, d=rho^2-rho^3 and
// P1=x1+x4, M1=x1-x4, P2=x2+x3, M2=x2-x3, all /2 at the end:
// Y0 = x0 + P1 + P2
// Y1 = x0 + ((aP1 + cM1) + (bP2 + dM2))/2
// Y2 = x0 + ((bP1 + dM1) + (aP2 - cM2))/2
// Y3 = x0 + ((bP1 - dM1) + (aP2 + cM2))/2
// Y4 = x0 + ((aP1 - cM1) + (bP2 - dM2))/2
// The inverse DFT-5 (root rho^{-1}) has the SAME shape with c,d -> -c,-d,
// so one routine serves both directions.
// ============================================================================
static u32 g_M5;
static u32 g_R = 1;
static u32 *g_outC = 0;
static u32 g_ninvS = 0, g_limS = 0;
static int g_fuse = 0, g_nt = 0;
static u32 *g_t5[4], *g_u5[4];
static u32 g_c5a, g_c5b, g_c5c, g_c5d, g_c5ci, g_c5di, g_inv2;
static u32 g_c5a_s, g_c5b_s, g_c5c_s, g_c5d_s, g_c5ci_s, g_c5di_s, g_inv2_s;
static void init_r5(u32 N) {
u32 M = N / 5;
g_M5 = M;
u32 w = powmod(31u, (NTP - 1) / N); // primitive N-th root
u32 wi = powmod(w, NTP - 2);
u32 rho = powmod(w, M);
u32 rho2 = mulmod(rho, rho), rho3 = mulmod(rho2, rho), rho4 = mulmod(rho3, rho);
g_c5a = (u32)(((u64)rho + rho4) % NTP);
g_c5b = (u32)(((u64)rho2 + rho3) % NTP);
g_c5c = (u32)(((u64)rho + NTP - rho4) % NTP);
g_c5d = (u32)(((u64)rho2 + NTP - rho3) % NTP);
g_c5ci = (u32)((NTP - (u64)g_c5c) % NTP);
g_c5di = (u32)((NTP - (u64)g_c5d) % NTP);
g_inv2 = (u32)(((u64)NTP + 1) / 2);
// ARM ch: fold the DFT-5 shared /2 into the CONSTANTS (no on-chain work, 4 fewer vshoup)
g_c5a = mulmod(g_c5a, g_inv2); g_c5b = mulmod(g_c5b, g_inv2);
g_c5c = mulmod(g_c5c, g_inv2); g_c5d = mulmod(g_c5d, g_inv2);
g_c5ci = mulmod(g_c5ci, g_inv2); g_c5di = mulmod(g_c5di, g_inv2);
g_c5a_s = shoup_c(g_c5a); g_c5b_s = shoup_c(g_c5b);
g_c5c_s = shoup_c(g_c5c); g_c5d_s = shoup_c(g_c5d);
g_c5ci_s = shoup_c(g_c5ci); g_c5di_s = shoup_c(g_c5di);
g_inv2_s = shoup_c(g_inv2);
u32 *pf = g_tabF + (sizeof(g_tabF) / 4 - (M + 1024));
u32 *pi = g_tabI + (sizeof(g_tabI) / 4 - (M + 1024));
for (int q = 0; q < 4; q++) { g_t5[q] = pf + (size_t)q * M; g_u5[q] = pi + (size_t)q * M; }
// Only series 0 (w^u, wi^u) is materialised -- 16.8 MB instead of 67 MB per direction,
// i.e. 100 MB less first-touch + generation + DRAM read. Series q is (series 0)^(q+1)
// and is rebuilt in-register with q vector Shoup multiplies per 8 lanes at the use site.
gen_geo_val(g_t5[0], M, 1, w);
gen_geo_val(g_u5[0], M, 1, wi);
for (int q = 1; q < 4; q++) { g_t5[q] = g_t5[0]; g_u5[q] = g_u5[0]; }
}
// vectorized DFT-5 (both directions; the caller supplies a,b and +/-c, +/-d)
TGT static inline void dft5v(__m256i x0, __m256i x1, __m256i x2, __m256i x3, __m256i x4,
__m256i &Y0, __m256i &Y1, __m256i &Y2, __m256i &Y3, __m256i &Y4,
__m256i pv, __m256i cA, __m256i cAs, __m256i cB, __m256i cBs,
__m256i cC, __m256i cCs, __m256i cD, __m256i cDs,
__m256i iv2, __m256i iv2s) {
__m256i P1 = vaddm(x1, x4, pv), M1 = vsubm(x1, x4, pv);
__m256i P2 = vaddm(x2, x3, pv), M2 = vsubm(x2, x3, pv);
__m256i aP1 = vshoup(P1, cA, cAs, pv), cM1 = vshoup(M1, cC, cCs, pv);
__m256i bP2 = vshoup(P2, cB, cBs, pv), dM2 = vshoup(M2, cD, cDs, pv);
__m256i bP1 = vshoup(P1, cB, cBs, pv), dM1 = vshoup(M1, cD, cDs, pv);
__m256i aP2 = vshoup(P2, cA, cAs, pv), cM2 = vshoup(M2, cC, cCs, pv);
Y0 = vaddm(vaddm(x0, P1, pv), P2, pv);
Y1 = vaddm(x0, vaddm(vaddm(aP1, cM1, pv), vaddm(bP2, dM2, pv), pv), pv);
Y2 = vaddm(x0, vaddm(vaddm(bP1, dM1, pv), vsubm(aP2, cM2, pv), pv), pv);
Y3 = vaddm(x0, vaddm(vsubm(bP1, dM1, pv), vaddm(aP2, cM2, pv), pv), pv);
Y4 = vaddm(x0, vaddm(vsubm(aP1, cM1, pv), vsubm(bP2, dM2, pv), pv), pv);
}
// apply out[q] = Y_q * w^{q u} for q=1..4 and store (forward twiddles)
#define ST5_FWD(a, j, Y0, Y1, Y2, Y3, Y4) \
{ stu256(((a) + (j)), (Y0)); \
__m256i w1 = ldu256((g_t5[0] + (j))); \
__m256i s1 = der_ws(w1); \
__m256i w2 = vshoup(w1, w1, s1, pv5); \
__m256i w3 = vshoup(w2, w1, s1, pv5); \
__m256i w4 = vshoup(w3, w1, s1, pv5); \
stu256(((a) + (j) + (size_t)g_M5), vshoup((Y1), w1, s1, pv5)); \
stu256(((a) + (j) + (size_t)2 * g_M5), vshoup((Y2), w2, der_ws(w2), pv5)); \
stu256(((a) + (j) + (size_t)3 * g_M5), vshoup((Y3), w3, der_ws(w3), pv5)); \
stu256(((a) + (j) + (size_t)4 * g_M5), vshoup((Y4), w4, der_ws(w4), pv5)); \
}
// scalar 5-point DFT-5 + twiddles for range tails (Shoup multiplication)
static inline void st5_scalar(u32 *a, const u32 x[5], u32 u, u32 M) {
u32 P1 = (u32)(((u64)x[1] + x[4]) % NTP), M1 = (u32)(((u64)x[1] + NTP - x[4]) % NTP);
u32 P2 = (u32)(((u64)x[2] + x[3]) % NTP), M2 = (u32)(((u64)x[2] + NTP - x[3]) % NTP);
u32 aP1 = shoup1s(P1, g_c5a, g_c5a_s), cM1 = shoup1s(M1, g_c5c, g_c5c_s);
u32 bP2 = shoup1s(P2, g_c5b, g_c5b_s), dM2 = shoup1s(M2, g_c5d, g_c5d_s);
u32 bP1 = shoup1s(P1, g_c5b, g_c5b_s), dM1 = shoup1s(M1, g_c5d, g_c5d_s);
u32 aP2 = shoup1s(P2, g_c5a, g_c5a_s), cM2 = shoup1s(M2, g_c5c, g_c5c_s);
u32 Y[5];
Y[0] = (u32)(((u64)x[0] + P1 + P2) % NTP);
Y[1] = (u32)(((u64)x[0] + (u32)(((u64)aP1 + cM1 + bP2 + dM2) % NTP)) % NTP);
Y[2] = (u32)(((u64)x[0] + (u32)(((u64)bP1 + dM1 + aP2 + NTP - cM2) % NTP)) % NTP);
Y[3] = (u32)(((u64)x[0] + (u32)(((u64)bP1 + NTP - dM1 + aP2 + cM2) % NTP)) % NTP);
Y[4] = (u32)(((u64)x[0] + (u32)(((u64)aP1 + NTP - cM1 + bP2 + NTP - dM2) % NTP)) % NTP);
a[u] = Y[0];
{ u32 w1 = g_t5[0][u], w2 = mulmod(w1, w1), w3 = mulmod(w2, w1), w4 = mulmod(w3, w1);
a[u + (size_t)M] = shoup1s(Y[1], w1, shoup_c(w1));
a[u + (size_t)2 * M] = shoup1s(Y[2], w2, shoup_c(w2));
a[u + (size_t)3 * M] = shoup1s(Y[3], w3, shoup_c(w3));
a[u + (size_t)4 * M] = shoup1s(Y[4], w4, shoup_c(w4)); }
}
// M180R2: 5-point DIF with x2 = x3 = x4 = 0 (Range 2 of stage5_dif_src). With those inputs
// P2 = x2+x3 = 0 and M2 = x2-x3 = 0 EXACTLY, so bP2, dM2, aP2, cM2 are all products with a
// structurally-zero operand and every term containing them vanishes identically.
// P1 = x1, M1 = x1 (vaddm(x1,0), vsubm(x1,0) are exact). Values are identical, so the
// result is byte-for-byte the same; only the arithmetic that was multiplying by zero is gone.
TGT static inline void dft5v_r2(__m256i x0, __m256i x1,
__m256i &Y0, __m256i &Y1, __m256i &Y2, __m256i &Y3, __m256i &Y4,
__m256i pv, __m256i cA, __m256i cAs, __m256i cB, __m256i cBs,
__m256i cC, __m256i cCs, __m256i cD, __m256i cDs) {
__m256i aP1 = vshoup(x1, cA, cAs, pv);
__m256i cM1 = vshoup(x1, cC, cCs, pv);
__m256i bP1 = vshoup(x1, cB, cBs, pv);
__m256i dM1 = vshoup(x1, cD, cDs, pv);
Y0 = vaddm(x0, x1, pv);
Y1 = vaddm(x0, vaddm(aP1, cM1, pv), pv);
Y2 = vaddm(x0, vaddm(bP1, dM1, pv), pv);
Y3 = vaddm(x0, vsubm(bP1, dM1, pv), pv);
Y4 = vaddm(x0, vsubm(aP1, cM1, pv), pv);
}
TGT static void stage5_dif_src(u32 *a, const u32 *src, u32 len, u32 M) {
const __m256i pv5 = _mm256_set1_epi32((int)NTP);
const __m256i cA = _mm256_set1_epi32((int)g_c5a), cAs = _mm256_set1_epi32((int)g_c5a_s);
const __m256i cB = _mm256_set1_epi32((int)g_c5b), cBs = _mm256_set1_epi32((int)g_c5b_s);
const __m256i cC = _mm256_set1_epi32((int)g_c5c), cCs = _mm256_set1_epi32((int)g_c5c_s);
const __m256i cD = _mm256_set1_epi32((int)g_c5d), cDs = _mm256_set1_epi32((int)g_c5d_s);
const __m256i iv2 = _mm256_set1_epi32((int)g_inv2), iv2s = _mm256_set1_epi32((int)g_inv2_s);
const __m256i zz = _mm256_setzero_si256();
// len is an INCLUSIVE index bound (coefficients 0..len): row t is valid for u + t*M <= len
u32 l4 = (len >= 4 * M) ? (len - 4 * M + 1) : 0; if (l4 > M) l4 = M;
u32 l3 = (len >= 3 * M) ? (len - 3 * M + 1) : 0; if (l3 > M) l3 = M;
u32 l2 = (len >= 2 * M) ? (len - 2 * M + 1) : 0; if (l2 > M) l2 = M;
u32 l1 = (len >= M) ? (len - M + 1) : 0; if (l1 > M) l1 = M;
u32 j = 0;
// Range 5: rows 0..4 all valid
for (; j + 8 <= l4; j += 8) {
pf5(src + j + PF_DIST); pf5(src + j + M + PF_DIST); pf5(src + j + 2 * M + PF_DIST);
pf5(src + j + 3 * M + PF_DIST); pf5(src + j + 4 * M + PF_DIST);
pf5(g_t5[0] + j + PF_DIST); pf5(g_t5[1] + j + PF_DIST); pf5(g_t5[2] + j + PF_DIST); pf5(g_t5[3] + j + PF_DIST);
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i x2 = ldu256((src + j + 2 * M));
__m256i x3 = ldu256((src + j + 3 * M));
__m256i x4 = ldu256((src + j + 4 * M));
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v(x0, x1, x2, x3, x4, Y0, Y1, Y2, Y3, Y4, pv5, cA, cAs, cB, cBs, cC, cCs, cD, cDs, iv2, iv2s);
ST5_FWD(a, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l4; j++) { u32 xx[5]; for (int t = 0; t < 5; t++) xx[t] = src[j + (size_t)t * M]; st5_scalar(a, xx, j, M); }
// Range 4: x4 = 0
for (; j + 8 <= l3; j += 8) {
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i x2 = ldu256((src + j + 2 * M));
__m256i x3 = ldu256((src + j + 3 * M));
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v(x0, x1, x2, x3, zz, Y0, Y1, Y2, Y3, Y4, pv5, cA, cAs, cB, cBs, cC, cCs, cD, cDs, iv2, iv2s);
ST5_FWD(a, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l3; j++) { u32 xx[5]; for (int t = 0; t < 5; t++) xx[t] = (t < 4 && j + (size_t)t * M <= len) ? src[j + (size_t)t * M] : 0; st5_scalar(a, xx, j, M); }
// Range 3: x3 = x4 = 0
for (; j + 8 <= l2; j += 8) {
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i x2 = ldu256((src + j + 2 * M));
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v(x0, x1, x2, zz, zz, Y0, Y1, Y2, Y3, Y4, pv5, cA, cAs, cB, cBs, cC, cCs, cD, cDs, iv2, iv2s);
ST5_FWD(a, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l2; j++) { u32 xx[5]; for (int t = 0; t < 5; t++) xx[t] = (t < 3 && j + (size_t)t * M <= len) ? src[j + (size_t)t * M] : 0; st5_scalar(a, xx, j, M); }
// Range 2: x2 = x3 = x4 = 0
for (; j + 8 <= l1; j += 8) {
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v_r2(x0, x1, Y0, Y1, Y2, Y3, Y4, pv5, cA, cAs, cB, cBs, cC, cCs, cD, cDs);
ST5_FWD(a, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l1; j++) { u32 xx[5]; for (int t = 0; t < 5; t++) xx[t] = (t < 2 && j + (size_t)t * M <= len) ? src[j + (size_t)t * M] : 0; st5_scalar(a, xx, j, M); }
// Range 1: only x0 (Y_q = x0 for all q)
for (; j + 8 <= M; j += 8) {
__m256i x0 = (j + 8 <= len + 1) ? ldu256((src + j)) : zz;
stu256((a + j), x0);
{
__m256i w1 = ldu256((g_t5[0] + j));
__m256i s1 = der_ws(w1);
__m256i w2 = vshoup(w1, w1, s1, pv5), w3 = vshoup(w2, w1, s1, pv5), w4 = vshoup(w3, w1, s1, pv5);
stu256((a + j + (size_t)M), vshoup(x0, w1, s1, pv5));
stu256((a + j + (size_t)2 * M), vshoup(x0, w2, der_ws(w2), pv5));
stu256((a + j + (size_t)3 * M), vshoup(x0, w3, der_ws(w3), pv5));
stu256((a + j + (size_t)4 * M), vshoup(x0, w4, der_ws(w4), pv5));
}
}
for (; j < M; j++) {
u32 x0 = (j <= len) ? src[j] : 0;
a[j] = x0;
u32 w1 = g_t5[0][j], w2 = mulmod(w1, w1), w3 = mulmod(w2, w1), w4 = mulmod(w3, w1);
a[j + (size_t)M] = shoup1s(x0, w1, shoup_c(w1));
a[j + (size_t)2 * M] = shoup1s(x0, w2, shoup_c(w2));
a[j + (size_t)3 * M] = shoup1s(x0, w3, shoup_c(w3));
a[j + (size_t)4 * M] = shoup1s(x0, w4, shoup_c(w4));
}
}
// inverse radix-5 DIT merge: t_q = a[qM+u]*wi^{qu}, then the inverse DFT-5
TGT static void merge5_dit(u32 *a, u32 M) {
const __m256i pv5 = _mm256_set1_epi32((int)NTP);
const __m256i cA = _mm256_set1_epi32((int)g_c5a), cAs = _mm256_set1_epi32((int)g_c5a_s);
const __m256i cB = _mm256_set1_epi32((int)g_c5b), cBs = _mm256_set1_epi32((int)g_c5b_s);
const __m256i cCi = _mm256_set1_epi32((int)g_c5ci), cCis = _mm256_set1_epi32((int)g_c5ci_s);
const __m256i cDi = _mm256_set1_epi32((int)g_c5di), cDis = _mm256_set1_epi32((int)g_c5di_s);
const __m256i iv2 = _mm256_set1_epi32((int)g_inv2), iv2s = _mm256_set1_epi32((int)g_inv2_s);
u32 s = 0;
for (; s + 8 <= M; s += 8) {
pf5(a + s + PF_DIST); pf5(a + s + M + PF_DIST); pf5(a + s + 2 * M + PF_DIST);
pf5(a + s + 3 * M + PF_DIST); pf5(a + s + 4 * M + PF_DIST);
pf5(g_u5[0] + s + PF_DIST); pf5(g_u5[1] + s + PF_DIST); pf5(g_u5[2] + s + PF_DIST); pf5(g_u5[3] + s + PF_DIST);
__m256i t0 = ldu256((a + s));
__m256i t[5]; t[0] = t0;
{
__m256i w1 = ldu256((g_u5[0] + s));
__m256i s1 = der_ws(w1);
__m256i w2 = vshoup(w1, w1, s1, pv5), w3 = vshoup(w2, w1, s1, pv5), w4 = vshoup(w3, w1, s1, pv5);
t[1] = vshoup(ldu256((a + s + (size_t)M)), w1, s1, pv5);
t[2] = vshoup(ldu256((a + s + (size_t)2 * M)), w2, der_ws(w2), pv5);
t[3] = vshoup(ldu256((a + s + (size_t)3 * M)), w3, der_ws(w3), pv5);
t[4] = vshoup(ldu256((a + s + (size_t)4 * M)), w4, der_ws(w4), pv5);
}
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v(t[0], t[1], t[2], t[3], t[4], Y0, Y1, Y2, Y3, Y4, pv5, cA, cAs, cB, cBs, cCi, cCis, cDi, cDis, iv2, iv2s);
if (g_fuse) {
// fold the final n^-1 scaling into this last stage and stream straight to the caller's c:
// removes the whole A->c pass (80 MB read) and the RFO on the A stores (84 MB).
const __m256i nv = _mm256_set1_epi32((int)g_ninvS), ns = _mm256_set1_epi32((int)shoup_c(g_ninvS));
Y0 = vshoup(Y0, nv, ns, pv5); Y1 = vshoup(Y1, nv, ns, pv5); Y2 = vshoup(Y2, nv, ns, pv5);
Y3 = vshoup(Y3, nv, ns, pv5); Y4 = vshoup(Y4, nv, ns, pv5);
u32 *C = g_outC + s;
if (g_nt) {
_mm256_stream_si256((__m256i *)(C), Y0);
_mm256_stream_si256((__m256i *)(C + M), Y1);
_mm256_stream_si256((__m256i *)(C + 2 * M), Y2);
_mm256_stream_si256((__m256i *)(C + 3 * M), Y3);
if (s < g_limS - 4 * M) _mm256_stream_si256((__m256i *)(C + 4 * M), Y4);
} else {
stu256((C), Y0);
stu256((C + M), Y1);
stu256((C + 2 * M), Y2);
stu256((C + 3 * M), Y3);
if (s < g_limS - 4 * M) stu256((C + 4 * M), Y4);
}
} else {
stu256((a + s), Y0);
stu256((a + s + M), Y1);
stu256((a + s + 2 * M), Y2);
stu256((a + s + 3 * M), Y3);
stu256((a + s + 4 * M), Y4);
}
}
for (; s < M; s++) {
u32 t[5]; t[0] = a[s];
{ u32 w1 = g_u5[0][s], w2 = mulmod(w1, w1), w3 = mulmod(w2, w1), w4 = mulmod(w3, w1);
t[1] = shoup1s(a[s + (size_t)M], w1, shoup_c(w1));
t[2] = shoup1s(a[s + (size_t)2 * M], w2, shoup_c(w2));
t[3] = shoup1s(a[s + (size_t)3 * M], w3, shoup_c(w3));
t[4] = shoup1s(a[s + (size_t)4 * M], w4, shoup_c(w4)); }
u32 P1 = (u32)(((u64)t[1] + t[4]) % NTP), M1 = (u32)(((u64)t[1] + NTP - t[4]) % NTP);
u32 P2 = (u32)(((u64)t[2] + t[3]) % NTP), M2 = (u32)(((u64)t[2] + NTP - t[3]) % NTP);
u32 Y[5];
Y[0] = (u32)(((u64)t[0] + P1 + P2) % NTP);
Y[1] = (u32)(((u64)t[0] + (u32)(((u64)mulmod(g_c5a, P1) + mulmod(g_c5ci, M1)) % NTP) + (u32)(((u64)mulmod(g_c5b, P2) + mulmod(g_c5di, M2)) % NTP)) % NTP);
Y[2] = (u32)(((u64)t[0] + (u32)(((u64)mulmod(g_c5b, P1) + mulmod(g_c5di, M1)) % NTP) + (u32)(((u64)mulmod(g_c5a, P2) + NTP - mulmod(g_c5ci, M2)) % NTP)) % NTP);
Y[3] = (u32)(((u64)t[0] + (u32)(((u64)mulmod(g_c5b, P1) + NTP - mulmod(g_c5di, M1)) % NTP) + (u32)(((u64)mulmod(g_c5a, P2) + mulmod(g_c5ci, M2)) % NTP)) % NTP);
Y[4] = (u32)(((u64)t[0] + (u32)(((u64)mulmod(g_c5a, P1) + NTP - mulmod(g_c5ci, M1)) % NTP) + (u32)(((u64)mulmod(g_c5b, P2) + NTP - mulmod(g_c5di, M2)) % NTP)) % NTP);
if (g_fuse) {
u32 *C = g_outC + s;
C[0] = (u32)((u64)Y[0] * g_ninvS % NTP);
C[M] = (u32)((u64)Y[1] * g_ninvS % NTP);
C[2 * M] = (u32)((u64)Y[2] * g_ninvS % NTP);
C[3 * M] = (u32)((u64)Y[3] * g_ninvS % NTP);
if (s < g_limS - 4 * M) C[4 * M] = (u32)((u64)Y[4] * g_ninvS % NTP);
} else {
a[s] = Y[0];
a[s + M] = Y[1]; a[s + 2 * M] = Y[2]; a[s + 3 * M] = Y[3]; a[s + 4 * M] = Y[4];
}
}
}
// ===================== mixed-radix driver =====================
static void ntt5_fwd(u32 *a, const u32 *src, u32 len) {
stage5_dif_src(a, src, len, g_M5);
for (int q = 0; q < 5; q++) ntt_fwd(a + (size_t)q * g_M5, 0, 0);
}
static void ntt5_inv(u32 *a, const u32 *Bmul) {
for (int q = 0; q < 5; q++) ntt_inv(a + (size_t)q * g_M5, Bmul ? (Bmul + (size_t)q * g_M5) : 0);
merge5_dit(a, g_M5);
}
// ===================== 1002e7: poly_multiply, n=m=1e7 (mixed radix 3*2^23) ====
static u32 g_dataA[5 * (1u << 22)] __attribute__((aligned(64)));
static u32 g_dataB[5 * (1u << 22)] __attribute__((aligned(64)));
static int g_ready = 0;
__attribute__((target("avx2"))) void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
g_R = (u32)(((u64)1 << 32) % NTP);
const u32 M = 1u << 22;
const u32 N = 5u * M;
if (!g_ready) { ntt_init(M); init_r5(N); g_ready = 1; }
u32 *A = g_dataA, *B = g_dataB;
u32 na = (u32)n, nb = (u32)m;
ntt5_fwd(A, a, na);
ntt5_fwd(B, b, nb);
{
u32 ninv = mulmod(powmod(N % NTP, NTP - 2), g_R);
g_ninvS = ninv; g_limS = na + nb + 1; g_outC = c;
g_nt = (((uintptr_t)c & 31) == 0);
g_fuse = (g_limS >= 4 * M);
}
ntt5_inv(A, B);
if (g_fuse) { _mm_sfence(); return; }
{
const __m256i pv = _mm256_set1_epi32((int)NTP);
u32 ninv = g_ninvS;
__m256i nv = _mm256_set1_epi32((int)ninv), ns = _mm256_set1_epi32((int)shoup_c(ninv));
u32 lim = na + nb + 1;
u32 i = 0;
while (i < lim && (((uintptr_t)(c + i)) & 31)) { c[i] = (u32)((u64)A[i] * ninv % NTP); i++; } // 32B-align the NT stores
for (; i + 8 <= lim; i += 8) {
__m256i x = vshoup(ldu256((A + i)), nv, ns, pv);
_mm256_stream_si256((__m256i *)(c + i), x); // c is write-only for us: no RFO traffic
}
for (; i < lim; i++) c[i] = (u32)((u64)A[i] * ninv % NTP);
_mm_sfence();
}
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 401.442 ms | 300 MB + 464 KB | Accepted | Score: 100 | 显示更多 |