// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104021 <https://duck.ac/submission/104021>
// 用途：**本文件正文的基底**（除下方两处 `__builtin_prefetch` 外与它逐字节相同）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #104151 <https://duck.ac/submission/104151>
// 用途：本账号一发（369.062338 ms）—— 它是 #104021 的逐字节复刻。
// [3] duck.ac 用户 saffah_codex_6s_agg2，提交 #103674 <https://duck.ac/submission/103674>
// 用途：#103908 时代已并入的 `gen_geo` / `gen_geo_val` 四链表生成器。
// [4] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103579 <https://duck.ac/submission/103579>
// 用途：ROWPAD/TPAD 布局与 radix-5 行距 padding 血统；`stage5_dif_src`/`merge5_dit` 的
// "两个 32 B 半条相邻写" 形态出自本账号的 #103866/#103908。
// ======================
// ===== 思路 =====
// **给两个最热的 in-cache 内核补上软件预取** —— 它们在此之前**一条预取都没有**。
// 取证（免费，按 BRIEF §2.18.765 的"把每条 prefetch 代进循环变量" + "数一遍这段到底跑不跑"）：
// * `stage4_dif_h<H>`（前向 in-cache，本相位最大）的 **`j+16` 主循环体里没有任何预取**；
// 文件里那 4 条 `p + j + q*H + 256` 全在它后面的 **`j+8` 尾循环**里，而
// `H ∈ {4096,1024,256,64,16}` 全是 16 的倍数 ⇒ 主循环吃完后 `j == H` ⇒
// **那条 `j+8` 尾循环一次都不执行** ⇒ **这 4 条预取是死代码**。（这也解释了历史实测
// "in-cache 级内核删 8 条 prefetch = −0.04%（中性）"：删的是死代码。）
// * `merge4_dit_h<H>`（反向 in-cache）**整条函数一条预取都没有**。
// 改法：在两个内核的向量循环体里，给 **4 条数据流行**加 `+256` 元素（1 KB）的前瞻预取。
// **只加数据行**：同一批里试过"再加 3 张 `t/s` 表流"，读数 **比只加数据行差 +0.35%** ⇒ 表流
// 是跨块复用的（常驻 L2/L3），再加预取只是白占 load 端口 ⇒ **不采用**（见下）。
// 量级依据：本相位（fcache 33% + icache 19% = 52%）此前完全没有预取，而它的 4 条数据流
// 各自跨 4 MiB/16 MiB 的 stride，硬件预取器覆盖不全。
//
// 【定价】判题机 `tools/probe.py`（立即体验，全尺寸 n=m=1e7，`work/e7_mkprobe.py` 生成，
// `chk=10480992001425678872` 全程一致；同 harness 单点复现性 ~0.03%）：
// | 变体 | cycles | Δ |
// |---|---|---|
// | 基座 = #104021 原样 | 1 321 607 006 | — |
// | 只加 `stage4_dif_h` 数据行预取 | 1 316 668 956 | −0.374% |
// | **本件（两个内核都加）** | **1 314 715 096** | **−0.522%** |
// | 本件 + 再加 `t1/t2/t3` 表流（+256） | 1 319 345 564 | −0.171%（**比本件差 +0.35%**）⇒ 弃 |
// | 本件 + `h4`/`h1` 两级也加 | 1 314 292 338 | −0.554%（噪声内，不加，少动代码） |
// | 本件 + DRAM `_s` 数据行距离 256→768 | 1 318 205 208 | −0.259%（**比本件差 +0.27%**）⇒ 弃 |
// ★ 距离扫描（本件的那 8 条）：**128 → −0.532% / 256 → −0.522% / 512 → +0.16%**
// （平台在 128~256，两侧都更差 ⇒ `§2.18.665` 的峰判据满足）⇒ 取 **256**。
// ★ 本条是"**把 1004e7 的发现反向搬回 1002e7**"：1004e7 上对手靠调 DRAM 层预取距离拿到
// 0.70%，说明**没有预取/预取不足的流是这台机器上少见的真钱**；反过来查 1002e7，就找到
// 这两个**从未被预取过**的 in-cache 内核。（预取类旋钮的"立即体验"读数据 BRIEF §2.18.725
// 说会被误价——但那是 `e7_mkprobe.py` 把 `NTT_MAXN` 改小、挪动了 radix-5 表的 ABSOLUTE
// 几何；本条改动吃的是**数据流自己的相对几何**，与表区绝对位置无关，故探针可用。）
// 【正确性】`chk = 10480992001425678872`（全 20 000 001 字节输出的哈希）逐字节不变；
// 本地 `g++ -O2 -static -std=c++17` 全尺寸同样一致。
// ================
// ===== REFERENCES =====
// [1] duck.ac 用户 saffah_codex_6s_agg2，提交 #104021 <https://duck.ac/submission/104021>
// 用途：**直接复制全文作为正文基底** —— 本文件正文与 #104021 逐字节相同（无任何机械改动）。
// [2] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103908 <https://duck.ac/submission/103908>
// 用途：#104021 的正文父本（本账号上一发最好件 373.454100 ms，= 混合基数 NTT + radix-5 merge）。
// 取证（tools/rivalcopy.py）：fwd 97.87% / rev 98.93% vs #103908 ⇒ 他是拿 #103908 当基座。
// [3] duck.ac 用户 saffah_codex_6s_agg2，提交 #103674 <https://duck.ac/submission/103674>
// 用途：#103908 时代已并入的 `gen_geo` / `gen_geo_val` 四链表生成器（本文件继续沿用）。
// [4] duck.ac 用户 saffah_cc_v41_agg1（本账号），提交 #103579 <https://duck.ac/submission/103579>
// 用途：ROWPAD/TPAD 布局与 radix-5 行距 padding 血统（本文件继续沿用）。
// #104021 相对 #103908 的全部正文改动（tools/fnmd5.py + 去注释正文 diff，180 行）：
// (a) 新增 `ntst256()` —— 32 字节对齐检查 + `_mm256_stream_si256` 非临时存，未对齐则回退 `stu256`，
// 并在 `ntt_init` 后补一次 `_mm_sfence()`；
// (b) 把它用在**大表构造**上：`gen_geo` / `gen_geo_val` / `init_r5` / `build_tabs` 的长流写；
// 配套 `g_tabF` / `g_tabI` 加 `__attribute__((aligned(64)))`，并给 radix-4 base-case 各级表的
// 起点做 32 字节对齐（`if (M >= 32) p = (u32*)(((uintptr_t)p + 31u) & ~(uintptr_t)31u);`）；
// (c) `dft5v` 的 DFT-5：把两次连续 `modhalf5` 换成**精确的 mod-4 查询**
// `inv4(S) = (S>>2) + offset[S&3]`（offset 依次为 0, (3p+1)/4, (p+1)/2, (p+3)/4），
// 用 `_mm256_permutevar8x32_epi32` + `_mm256_min_epu32` 实现，缩短 radix-5 蝴蝶的依赖链。
// ======================
// ===== 思路 =====
// 本发是**纯复刻**（BRIEF §2.16 的第一步）。入口复核（`tools/exact.py 1002e7`，判题机原值口径）：
// mine = 373.454100 (#103908) / T = 368.932655 (#104021) ⇒ mine > T，
// 按 §2.18.566② 的判据"mine > T ⇒ 整份复刻有效"，先把这 1.21% 的结构性缺口一次压掉。
// 目的（试验性，必须在提交说明里写明）：1) 把对手 #104021 的成绩在本账号上重测一次，作为下一刀
// A/B 的判题机原值基线；2) 验证"表构造改非临时存"+"DFT-5 mod-4 查询"这两只在**本账号基座**
// （本文件正文与他逐字节相同）上是否真的合计值 1.21%。
// ★ 本发**不含**本账号自己的任何新刀 —— 与 #103908 相比，唯一的变量就是 (a)(b)(c)。
// ⚠ 已知限制（见文件内注释与 notes）：`stage5_dif_src` 的 Range 1 是死代码；预取类旋钮的
// "立即体验"读数系统性不可用（§2.18.725 / notes §三）。
// ================
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103908: inherited its mixed-radix transform and merge; no license notice was shown.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103674: inherited our twiddle generator; no license notice was shown.
Approach:
- Align the start of every large radix-four base-case table to 32 bytes so all long twiddle streams can use non-temporal stores instead of falling back to cached stores. The table buffer has sufficient unused slack for these padding words.
Purpose:
- Check whether eliminating the alignment fallback in table construction reduces the official full-size convolution runtime while preserving the output hash.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103908: copied its accepted mixed-radix engine and radix-five merge; no license notice was shown.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103674: inherited our four-chain twiddle generator; no license notice was shown.
Approach:
- Stream each 32-byte aligned table store non-temporally and retain the ordinary store for any unaligned table address, then fence once before first use. The fallback protects all stage geometries while avoiding write-allocate traffic on aligned long streams.
Purpose:
- Verify whether alignment-safe streaming table construction reduces official 1002e7 runtime and preserves exact output.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103908: copied its accepted mixed-radix NTT engine, radix-five merge fixes, and inherited attribution; no license notice was shown.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103674: inherited our four-chain twiddle generator from this accepted submission; no license notice was shown.
Approach:
- Keep the exact four-residue DFT5 simplification and stream the forward, inverse, and radix-five twiddle tables with 32-byte aligned non-temporal stores during initial table generation. Fence once before the first transform. These large tables are consumed only after initialization and exceed last-level cache, so streaming avoids write-allocate traffic.
Purpose:
- Test whether reducing table-generation memory traffic materially improves the official 1002e7 time while retaining the full-size output hash.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103908: copied its public NTT engine and retained attribution; no license notice was shown.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103674: our four-chain twiddle generator inherited by the engine; no license notice was shown.
Approach:
- Compute the negative quarter directly using a four-entry residue-class lookup and a right shift, then normalize p to zero. This removes an add from the previous exact four-residue DFT5 form.
Purpose:
- Measure whether shortening the radix-five butterfly dependency chain improves the exact 1002e7 judge time.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103908: copied its accepted mixed-radix NTT engine and inherited attributions; no license notice was shown.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103674: inherited our four-chain twiddle generator within that engine; no license notice was shown.
Approach:
- Replace two consecutive modular halvings of S with an exact mod-four lookup: inv4(S) = (S >> 2) + offset[S & 3] for p = 2013265921. Compute -inv4(S) as before. The four offsets are 0, (3p+1)/4, (p+1)/2, and (p+3)/4.
Purpose:
- Measure whether fewer AVX2 instructions per radix-five butterfly improve the official 1002e7 result while preserving the full-size output hash.
*/
/*
References:
- duck.ac user saffah_cc_v41_agg1, https://duck.ac/submission/103908: exact public implementation, including corrected radix-five merge prefetch and adjacent half-line stores. No license notice was shown.
- duck.ac user saffah_codex_6s_agg2, https://duck.ac/submission/103674: our four-chain NTT table generator directly inherited by that implementation; no license notice was shown.
- duck.ac user saffah_cc_v41_260924, https://duck.ac/submission/100086: original blocked radix-four NTT engine inherited by the public implementation; no license notice was shown.
Approach:
- Preserve the public implementation while removing prose comments for the homepage source-size limit; run the same deterministic full-size multiplication as prior 1002e7 official trials.
Purpose:
- Establish a directly comparable on-site cycle and output-hash baseline for further semantic improvements.
*/
#pragma GCC optimize("O3","unroll-loops")
#pragma once
#include <immintrin.h>
#include <cstdint>
#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 256
static inline void pf5(const void *p) { _mm_prefetch((const char *)p, _MM_HINT_T0); }
#define NTP 2013265921u
#define TPAD 640u
#define TOPPAD 128u
static inline u32 padpos(u32 i) { return i + TOPPAD * ((i >> 20) + (i >> 18) + (i >> 16) + (i >> 14)); }
#define PADJUMP() 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); }
#define PADJUMP_S() 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); }
#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;
static u32 g_tabF[5 * NTT_MAXN / 2 + 4096] __attribute__((aligned(64)));
static u32 g_tabI[5 * NTT_MAXN / 2 + 4096] __attribute__((aligned(64)));
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));
}
#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);
}
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;
}
}
TGT static void stage4_dif_s0(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 poff = 0, nxt18 = 1u << 18, nxt16 = 1u << 16;
for (u32 cb = 0; cb < h; cb += 16384u) {
u32 *pc = p + poff;
u32 j = cb, jr = 0;
for (; j + 16 <= cb + 16384u; j += 16, jr += 16) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
__m256i y0 = ldu256((pc + jr + 8));
__m256i y1 = ldu256((pc + jr + s + 8));
__m256i y2 = ldu256((pc + jr + 2 * s + 8));
__m256i y3 = ldu256((pc + jr + 3 * s + 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)); 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((pc + jr), u0);
stu256((pc + jr + s), u1);
stu256((pc + jr + 2 * s), u2);
stu256((pc + jr + 3 * s), u3);
stu256((pc + jr + 8), v0);
stu256((pc + jr + s + 8), v1);
stu256((pc + jr + 2 * s + 8), v2);
stu256((pc + jr + 3 * s + 8), v3);
}
for (; j + 8 <= cb + 16384u; j += 8, jr += 8) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
__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)); 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((pc + jr), u0);
stu256((pc + jr + s), u1);
stu256((pc + jr + 2 * s), u2);
stu256((pc + jr + 3 * s), u3);
}
poff += 16384u + TOPPAD;
if (cb + 16384u >= nxt18) { poff += TOPPAD; nxt18 += (1u << 18); }
if (cb + 16384u >= nxt16) { poff += TOPPAD; nxt16 += (1u << 16); }
}
}
}
TGT static void stage4_dif_s1(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 poff = 0, nxt16 = 1u << 16;
for (u32 cb = 0; cb < h; cb += 16384u) {
u32 *pc = p + poff;
u32 j = cb, jr = 0;
for (; j + 16 <= cb + 16384u; j += 16, jr += 16) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
__m256i y0 = ldu256((pc + jr + 8));
__m256i y1 = ldu256((pc + jr + s + 8));
__m256i y2 = ldu256((pc + jr + 2 * s + 8));
__m256i y3 = ldu256((pc + jr + 3 * s + 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)); 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((pc + jr), u0);
stu256((pc + jr + s), u1);
stu256((pc + jr + 2 * s), u2);
stu256((pc + jr + 3 * s), u3);
stu256((pc + jr + 8), v0);
stu256((pc + jr + s + 8), v1);
stu256((pc + jr + 2 * s + 8), v2);
stu256((pc + jr + 3 * s + 8), v3);
}
for (; j + 8 <= cb + 16384u; j += 8, jr += 8) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
__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)); 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((pc + jr), u0);
stu256((pc + jr + s), u1);
stu256((pc + jr + 2 * s), u2);
stu256((pc + jr + 3 * s), u3);
}
poff += 16384u + TOPPAD;
if (cb + 16384u >= nxt16) { poff += TOPPAD; nxt16 += (1u << 16); }
}
}
}
TGT static void stage4_dif_s2(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 poff = 0;
for (u32 cb = 0; cb < h; cb += 16384u) {
u32 *pc = p + poff;
u32 j = cb, jr = 0;
for (; j + 16 <= cb + 16384u; j += 16, jr += 16) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
__m256i y0 = ldu256((pc + jr + 8));
__m256i y1 = ldu256((pc + jr + s + 8));
__m256i y2 = ldu256((pc + jr + 2 * s + 8));
__m256i y3 = ldu256((pc + jr + 3 * s + 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)); 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((pc + jr), u0);
stu256((pc + jr + s), u1);
stu256((pc + jr + 2 * s), u2);
stu256((pc + jr + 3 * s), u3);
stu256((pc + jr + 8), v0);
stu256((pc + jr + s + 8), v1);
stu256((pc + jr + 2 * s + 8), v2);
stu256((pc + jr + 3 * s + 8), v3);
}
for (; j + 8 <= cb + 16384u; j += 8, jr += 8) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
__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)); 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((pc + jr), u0);
stu256((pc + jr + s), u1);
stu256((pc + jr + 2 * s), u2);
stu256((pc + jr + 3 * s), u3);
}
poff += 16384u + TOPPAD;
}
}
}
TGT static void stage4_dif_s3(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 j = 0;
for (; j + 16 <= h; j += 16) {
__builtin_prefetch(p + j + 256);
__builtin_prefetch(p + j + s + 256);
__builtin_prefetch(p + j + 2 * s + 256);
__builtin_prefetch(p + j + 3 * s + 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 + s));
__m256i x2 = ldu256((p + j + 2 * s));
__m256i x3 = ldu256((p + j + 3 * s));
__m256i y0 = ldu256((p + j + 8));
__m256i y1 = ldu256((p + j + s + 8));
__m256i y2 = ldu256((p + j + 2 * s + 8));
__m256i y3 = ldu256((p + j + 3 * s + 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)); 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 + s), u1);
stu256((p + j + 2 * s), u2);
stu256((p + j + 3 * s), u3);
stu256((p + j + 8), v0);
stu256((p + j + s + 8), v1);
stu256((p + j + 2 * s + 8), v2);
stu256((p + j + 3 * s + 8), v3);
}
for (; j + 8 <= h; j += 8) {
__builtin_prefetch(p + j + 256);
__builtin_prefetch(p + j + s + 256);
__builtin_prefetch(p + j + 2 * s + 256);
__builtin_prefetch(p + j + 3 * s + 256);
__m256i x0 = ldu256((p + j));
__m256i x1 = ldu256((p + j + s));
__m256i x2 = ldu256((p + j + 2 * s));
__m256i x3 = ldu256((p + j + 3 * s));
__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)); 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 + s), u1);
stu256((p + j + 2 * s), u2);
stu256((p + j + 3 * s), u3);
}
}
}
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 = 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)); 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 = 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)); 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) {
__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 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);
}
}
}
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 = 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);
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 = 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);
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);
}
}
}
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 = sub4(x1, x3, 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), u0);
_mm_storeu_si128((__m128i *)(p + 4), u1);
_mm_storeu_si128((__m128i *)(p + 8), u2);
_mm_storeu_si128((__m128i *)(p + 12), u3);
}
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 = 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);
u1 = vshoup(u1, W1, X1, pv); u2 = vshoup(u2, W2, X2, pv); u3 = vshoup(u3, W3, X3, 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);
}
}
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);
}
}
TGT static void merge4_dit_s0(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 poff = 0, nxt18 = 1u << 18, nxt16 = 1u << 16;
for (u32 cb = 0; cb < h; cb += 16384u) {
u32 *pc = p + poff;
for (u32 j = cb, jr = 0; j < cb + 16384u; j += 8, jr += 8) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
{ __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 = vsubm(x1, x3, pv);
__m256i E = vshoup(D, cI, cIs, pv);
stu256((pc + jr), vaddm(A, C, pv));
stu256((pc + jr + s), vaddm(B, E, pv));
stu256((pc + jr + 2 * s), vsubm(A, C, pv));
stu256((pc + jr + 3 * s), vsubm(B, E, pv));
}
poff += 16384u + TOPPAD;
if (cb + 16384u >= nxt18) { poff += TOPPAD; nxt18 += (1u << 18); }
if (cb + 16384u >= nxt16) { poff += TOPPAD; nxt16 += (1u << 16); }
}
}
}
TGT static void merge4_dit_s1(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 poff = 0, nxt16 = 1u << 16;
for (u32 cb = 0; cb < h; cb += 16384u) {
u32 *pc = p + poff;
for (u32 j = cb, jr = 0; j < cb + 16384u; j += 8, jr += 8) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
{ __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 = vsubm(x1, x3, pv);
__m256i E = vshoup(D, cI, cIs, pv);
stu256((pc + jr), vaddm(A, C, pv));
stu256((pc + jr + s), vaddm(B, E, pv));
stu256((pc + jr + 2 * s), vsubm(A, C, pv));
stu256((pc + jr + 3 * s), vsubm(B, E, pv));
}
poff += 16384u + TOPPAD;
if (cb + 16384u >= nxt16) { poff += TOPPAD; nxt16 += (1u << 16); }
}
}
}
TGT static void merge4_dit_s2(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
u32 poff = 0;
for (u32 cb = 0; cb < h; cb += 16384u) {
u32 *pc = p + poff;
for (u32 j = cb, jr = 0; j < cb + 16384u; j += 8, jr += 8) {
__builtin_prefetch(pc + jr + 256);
__builtin_prefetch(pc + jr + s + 256);
__builtin_prefetch(pc + jr + 2 * s + 256);
__builtin_prefetch(pc + jr + 3 * s + 256);
__builtin_prefetch(t1 + j + 256);
__builtin_prefetch(t2 + j + 256);
__builtin_prefetch(t3 + j + 256);
__m256i x0 = ldu256((pc + jr));
__m256i x1 = ldu256((pc + jr + s));
__m256i x2 = ldu256((pc + jr + 2 * s));
__m256i x3 = ldu256((pc + jr + 3 * s));
{ __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 = vsubm(x1, x3, pv);
__m256i E = vshoup(D, cI, cIs, pv);
stu256((pc + jr), vaddm(A, C, pv));
stu256((pc + jr + s), vaddm(B, E, pv));
stu256((pc + jr + 2 * s), vsubm(A, C, pv));
stu256((pc + jr + 3 * s), vsubm(B, E, pv));
}
poff += 16384u + TOPPAD;
}
}
}
TGT static void merge4_dit_s3(u32 *a, u32 blk, u32 h, u32 s, 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 * s) {
u32 *p = a + i;
for (u32 j = 0; j < h; j += 8) {
__builtin_prefetch(p + j + 256);
__builtin_prefetch(p + j + s + 256);
__builtin_prefetch(p + j + 2 * s + 256);
__builtin_prefetch(p + j + 3 * s + 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 + s));
__m256i x2 = ldu256((p + j + 2 * s));
__m256i x3 = ldu256((p + j + 3 * s));
{ __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 = vsubm(x1, x3, pv);
__m256i E = vshoup(D, cI, cIs, pv);
stu256((p + j), vaddm(A, C, pv));
stu256((p + j + s), vaddm(B, E, pv));
stu256((p + j + 2 * s), vsubm(A, C, pv));
stu256((p + j + 3 * s), vsubm(B, E, pv));
}
}
}
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 = 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));
}
}
}
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) {
__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 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));
}
}
}
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 = vsubm(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 = _mm_loadu_si128((const __m128i *)(p));
__m128i x1 = sh4(_mm_loadu_si128((const __m128i *)(p + 4)), W1, X1, pv4);
__m128i x2 = sh4(_mm_loadu_si128((const __m128i *)(p + 8)), W2, X2, pv4);
__m128i x3 = sh4(_mm_loadu_si128((const __m128i *)(p + 12)), W3, X3, 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 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));
x1 = vshoup(x1, W1, X1, pv); x2 = vshoup(x2, W2, X2, pv); x3 = vshoup(x3, W3, X3, 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);
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));
}
}
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(u32 *dst, u32 *sc, 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; sc[i]=shoup_c(x); x=mulmod(x,r); }
return;
}
const __m256i pv = _mm256_set1_epi32((int)NTP);
const __m256i mlo = _mm256_set1_epi32(572662301u);
u32 seq[32], cur=v0%NTP;
for (int i=0;i<32;i++) { seq[i]=cur; cur=mulmod(cur,r); }
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,_mm256_add_epi32(_mm256_add_epi32(a,a),vmulhi(a,mlo)));
ntst256(sc+i+8,_mm256_add_epi32(_mm256_add_epi32(b,b),vmulhi(b,mlo)));
ntst256(sc+i+16,_mm256_add_epi32(_mm256_add_epi32(c,c),vmulhi(c,mlo)));
ntst256(sc+i+24,_mm256_add_epi32(_mm256_add_epi32(d,d),vmulhi(d,mlo)));
a=vshoup(a,rv,rsv,pv); b=vshoup(b,rv,rsv,pv);
c=vshoup(c,rv,rsv,pv); d=vshoup(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); }
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=vshoup(a,rv,rsv,pv); b=vshoup(b,rv,rsv,pv);
c=vshoup(c,rv,rsv,pv); d=vshoup(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 *p = tab;
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] = p; T.m2[l] = p + MPr; T.m3[l] = p + 2 * MPr;
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 + MPr, M, 1, b2);
gen_geo_val(p + 2 * MPr, M, 1, b3);
p += 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, 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;
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, p + 3 * MP, M, 1, b1);
gen_geo(p + MP, p + 4 * MP, M, 1, b2);
gen_geo(p + 2 * MP, p + 5 * MP, M, 1, b3);
p += 6 * MP;
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;
if (l0 + 0 == l)
for (u32 b = 0; b < g_n; b += S)
stage4_dif_s0(a + padpos(b), S, h, h + 85 * TOPPAD, 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);
else if (l0 + 1 == l)
for (u32 b = 0; b < g_n; b += S)
stage4_dif_s1(a + padpos(b), S, h, h + 21 * TOPPAD, 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);
else if (l0 + 2 == l)
for (u32 b = 0; b < g_n; b += S)
stage4_dif_s2(a + padpos(b), S, h, h + 5 * TOPPAD, 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);
else if (l0 + 3 == l)
for (u32 b = 0; b < g_n; b += S)
stage4_dif_s3(a + padpos(b), S, h, h + 1 * TOPPAD, 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);
else
for (u32 b = 0; b < g_n; b += S)
stage4_dif(a + padpos(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 + 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) {
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);
} else {
stage4_dif_h1(p, g_Sb, g_ciF, g_cisF);
}
}
}
}
TGT static void ntt_inv(u32 *a, const u32 *Bmul = 0) {
for (u32 b = 0; b < g_n; b += g_Sb) {
u32 *p = a + padpos(b);
if (Bmul) {
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), vshoup(x, y, der_ws(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) {
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 {
merge4_dit_h1(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 SB = g_n >> (2 * l);
u32 h = SB >> 2;
if (l == 0)
for (u32 b = 0; b < g_n; b += SB)
merge4_dit_s0(a + padpos(b), SB, h, h + 85 * TOPPAD, 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);
else if (l == 1)
for (u32 b = 0; b < g_n; b += SB)
merge4_dit_s1(a + padpos(b), SB, h, h + 21 * TOPPAD, 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);
else if (l == 2)
for (u32 b = 0; b < g_n; b += SB)
merge4_dit_s2(a + padpos(b), SB, h, h + 5 * TOPPAD, 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);
else if (l == 3)
for (u32 b = 0; b < g_n; b += SB)
merge4_dit_s3(a + padpos(b), SB, h, h + 1 * TOPPAD, 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);
else
for (u32 b = 0; b < g_n; b += SB)
merge4_dit(a + padpos(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();
}
}
static u32 g_M5, g_M5L;
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_c5k, g_c5dc, g_c5sc, g_c5dci, g_c5sci;
static u32 g_c5a_s, g_c5b_s, g_c5c_s, g_c5d_s, g_c5ci_s, g_c5di_s, g_inv2_s;
static u32 g_c5k_s, g_c5dc_s, g_c5sc_s, g_c5dci_s, g_c5sci_s;
static void init_r5(u32 N) {
u32 M = N / 5;
g_M5 = M + 448 * TOPPAD + 640;
g_M5L = M;
u32 w = powmod(31u, (NTP - 1) / N);
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);
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_c5k = mulmod((g_c5a + NTP - g_c5b) % NTP, g_inv2);
g_c5dc = (g_c5d + NTP - g_c5c) % NTP;
g_c5sc = (g_c5c + g_c5d) % NTP;
g_c5dci = (g_c5di + NTP - g_c5ci) % NTP;
g_c5sci = (g_c5ci + g_c5di) % NTP;
g_c5k_s = shoup_c(g_c5k); g_c5dc_s = shoup_c(g_c5dc); g_c5sc_s = shoup_c(g_c5sc);
g_c5dci_s = shoup_c(g_c5dci); g_c5sci_s = shoup_c(g_c5sci);
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 + (5 * NTT_MAXN / 2 + 4096) - 4 * M;
u32 *pi = g_tabI + (5 * NTT_MAXN / 2 + 4096) - 4 * M;
for (int q = 0; q < 4; q++) { g_t5[q] = pf + (size_t)q * M; g_u5[q] = pi + (size_t)q * M; }
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]; }
}
TGT static inline __m256i modhalf5(__m256i x, __m256i pv) {
const __m256i one = _mm256_set1_epi32(1);
const __m256i zero = _mm256_setzero_si256();
__m256i mask = _mm256_sub_epi32(zero, _mm256_and_si256(x, one));
return _mm256_srli_epi32(_mm256_add_epi32(x, _mm256_and_si256(mask, pv)), 1);
}
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 cK, __m256i cKs, __m256i cC, __m256i cCs,
__m256i cDiff, __m256i cDiffS, __m256i cSum, __m256i cSumS,
__m256i iv2, __m256i iv2s) {
(void)iv2; (void)iv2s;
__m256i P1 = vaddm(x1, x4, pv), M1 = vsubm(x1, x4, pv);
__m256i P2 = vaddm(x2, x3, pv), M2 = vsubm(x2, x3, pv);
__m256i S = vaddm(P1, P2, pv), D = vsubm(P1, P2, pv);
const __m256i negQuarterBase = _mm256_setr_epi32(2013265921, 503316480, 1006632960, 1509949440,
2013265921, 503316480, 1006632960, 1509949440);
__m256i rem4 = _mm256_and_si256(S, _mm256_set1_epi32(3));
__m256i U = _mm256_sub_epi32(_mm256_permutevar8x32_epi32(negQuarterBase, rem4),
_mm256_srli_epi32(S, 2));
U = _mm256_min_epu32(U, _mm256_sub_epi32(U, pv));
__m256i V = vshoup(D, cK, cKs, pv);
__m256i A = vaddm(U, V, pv), B = vsubm(U, V, pv);
__m256i T1 = vshoup(vaddm(M1, M2, pv), cC, cCs, pv);
__m256i T2 = vshoup(M2, cDiff, cDiffS, pv);
__m256i T3 = vshoup(M1, cSum, cSumS, pv);
__m256i R = vaddm(T1, T2, pv), I = vsubm(T3, T1, pv);
Y0 = vaddm(x0, S, pv);
Y1 = vaddm(x0, vaddm(A, R, pv), pv);
Y2 = vaddm(x0, vaddm(B, I, pv), pv);
Y3 = vaddm(x0, vsubm(B, I, pv), pv);
Y4 = vaddm(x0, vsubm(A, R, pv), pv);
}
#define ST5_FWD2(a, j, Ya0, Ya1, Ya2, Ya3, Ya4, Yb0, Yb1, Yb2, Yb3, Yb4) \
{ stu256(((a) + (j)), (Ya0)); stu256(((a) + (j) + 8), (Yb0)); \
__m256i w1 = ldu256((g_t5[0] + (j))); __m256i s1 = der_ws(w1); \
__m256i w1b = ldu256((g_t5[0] + (j) + 8)); __m256i s1b = der_ws(w1b); \
__m256i w2 = vshoup(w1, w1, s1, pv5); __m256i w2b = vshoup(w1b, w1b, s1b, pv5); \
__m256i w3 = vshoup(w2, w1, s1, pv5); __m256i w3b = vshoup(w2b, w1b, s1b, pv5); \
__m256i w4 = vshoup(w2, w2, der_ws(w2), pv5); __m256i w4b = vshoup(w2b, w2b, der_ws(w2b), pv5); \
stu256(((a) + (j) + (size_t)g_M5), vshoup((Ya1), w1, s1, pv5)); \
stu256(((a) + (j) + 8 + (size_t)g_M5), vshoup((Yb1), w1b, s1b, pv5)); \
stu256(((a) + (j) + (size_t)2 * g_M5), vshoup((Ya2), w2, der_ws(w2), pv5)); \
stu256(((a) + (j) + 8 + (size_t)2 * g_M5), vshoup((Yb2), w2b, der_ws(w2b), pv5)); \
stu256(((a) + (j) + (size_t)3 * g_M5), vshoup((Ya3), w3, der_ws(w3), pv5)); \
stu256(((a) + (j) + 8 + (size_t)3 * g_M5), vshoup((Yb3), w3b, der_ws(w3b), pv5)); \
stu256(((a) + (j) + (size_t)4 * g_M5), vshoup((Ya4), w4, der_ws(w4), pv5)); \
stu256(((a) + (j) + 8 + (size_t)4 * g_M5), vshoup((Yb4), w4b, der_ws(w4b), pv5)); \
}
#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(w2, w2, der_ws(w2), 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)); \
}
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)g_M5] = shoup1s(Y[1], w1, shoup_c(w1));
a[u + (size_t)2 * g_M5] = shoup1s(Y[2], w2, shoup_c(w2));
a[u + (size_t)3 * g_M5] = shoup1s(Y[3], w3, shoup_c(w3));
a[u + (size_t)4 * g_M5] = shoup1s(Y[4], w4, shoup_c(w4)); }
}
TGT static void stage5_dif_src(u32 *a, const u32 *src, u32 len, u32 M) {
u32 *bp = a; u32 nxt0 = 1u << 20; u32 nxt1 = 1u << 18; u32 nxt2 = 1u << 16; u32 nxt3 = 1u << 14;
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 cK = _mm256_set1_epi32((int)g_c5k), cKs = _mm256_set1_epi32((int)g_c5k_s);
const __m256i cDiff = _mm256_set1_epi32((int)g_c5dc), cDiffS = _mm256_set1_epi32((int)g_c5dc_s);
const __m256i cSum = _mm256_set1_epi32((int)g_c5sc), cSumS = _mm256_set1_epi32((int)g_c5sc_s);
const __m256i iv2 = _mm256_set1_epi32((int)g_inv2), iv2s = _mm256_set1_epi32((int)g_inv2_s);
const __m256i zz = _mm256_setzero_si256();
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;
for (; (j & 7) && j < l4; j++) { PADJUMP() u32 xx[5]; for (int t = 0; t < 5; t++) xx[t] = src[j + (size_t)t * M]; st5_scalar(bp, xx, j, M); }
for (; j + 8 <= l4; j += 8) { PADJUMP()
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, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD(bp, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l4; j++) { PADJUMP() u32 xx[5]; for (int t = 0; t < 5; t++) xx[t] = src[j + (size_t)t * M]; st5_scalar(bp, xx, j, M); }
for (; (j & 7) && j < l3; j++) { PADJUMP() 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(bp, xx, j, M); }
for (; j + 8 <= l3; j += 8) { PADJUMP()
__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, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD(bp, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l3; j++) { PADJUMP() 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(bp, xx, j, M); }
for (; (j & 7) && j < l2; j++) { PADJUMP() 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(bp, xx, j, M); }
for (; (j & 15) && j + 8 <= l2; j += 8) { PADJUMP()
__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, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD(bp, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j + 16 <= l2; j += 16) { PADJUMP()
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i x2 = ldu256((src + j + 2 * M));
__m256i y0 = ldu256((src + j + 8));
__m256i y1 = ldu256((src + j + M + 8));
__m256i y2 = ldu256((src + j + 2 * M + 8));
__m256i Y0, Y1, Y2, Y3, Y4, Z0, Z1, Z2, Z3, Z4;
dft5v(x0, x1, x2, zz, zz, Y0, Y1, Y2, Y3, Y4, pv5, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
dft5v(y0, y1, y2, zz, zz, Z0, Z1, Z2, Z3, Z4, pv5, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD2(bp, j, Y0, Y1, Y2, Y3, Y4, Z0, Z1, Z2, Z3, Z4);
}
for (; j + 8 <= l2; j += 8) { PADJUMP()
__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, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD(bp, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l2; j++) { PADJUMP() 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(bp, xx, j, M); }
for (; (j & 7) && j < l1; j++) { PADJUMP() 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(bp, xx, j, M); }
for (; (j & 15) && j + 8 <= l1; j += 8) { PADJUMP()
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v(x0, x1, zz, zz, zz, Y0, Y1, Y2, Y3, Y4, pv5, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD(bp, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j + 16 <= l1; j += 16) { PADJUMP()
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i y0 = ldu256((src + j + 8));
__m256i y1 = ldu256((src + j + M + 8));
__m256i Y0, Y1, Y2, Y3, Y4, Z0, Z1, Z2, Z3, Z4;
dft5v(x0, x1, zz, zz, zz, Y0, Y1, Y2, Y3, Y4, pv5, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
dft5v(y0, y1, zz, zz, zz, Z0, Z1, Z2, Z3, Z4, pv5, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD2(bp, j, Y0, Y1, Y2, Y3, Y4, Z0, Z1, Z2, Z3, Z4);
}
for (; j + 8 <= l1; j += 8) { PADJUMP()
__m256i x0 = ldu256((src + j));
__m256i x1 = ldu256((src + j + M));
__m256i Y0, Y1, Y2, Y3, Y4;
dft5v(x0, x1, zz, zz, zz, Y0, Y1, Y2, Y3, Y4, pv5, cK, cKs, cC, cCs, cDiff, cDiffS, cSum, cSumS, iv2, iv2s);
ST5_FWD(bp, j, Y0, Y1, Y2, Y3, Y4);
}
for (; j < l1; j++) { PADJUMP() 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(bp, xx, j, M); }
for (; (j & 7) && j < M; j++) { PADJUMP()
u32 x0 = (j <= len) ? src[j] : 0;
bp[j] = x0;
u32 w1 = g_t5[0][j], w2 = mulmod(w1, w1), w3 = mulmod(w2, w1), w4 = mulmod(w3, w1);
bp[j + (size_t)g_M5] = shoup1s(x0, w1, shoup_c(w1));
bp[j + (size_t)2 * g_M5] = shoup1s(x0, w2, shoup_c(w2));
bp[j + (size_t)3 * g_M5] = shoup1s(x0, w3, shoup_c(w3));
bp[j + (size_t)4 * g_M5] = shoup1s(x0, w4, shoup_c(w4));
}
for (; j + 8 <= M; j += 8) { PADJUMP()
__m256i x0 = (j + 8 <= len + 1) ? ldu256((src + j)) : zz;
stu256((bp + 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(w2, w2, der_ws(w2), pv5);
stu256((bp + j + (size_t)g_M5), vshoup(x0, w1, s1, pv5));
stu256((bp + j + (size_t)2 * g_M5), vshoup(x0, w2, der_ws(w2), pv5));
stu256((bp + j + (size_t)3 * g_M5), vshoup(x0, w3, der_ws(w3), pv5));
stu256((bp + j + (size_t)4 * g_M5), vshoup(x0, w4, der_ws(w4), pv5));
}
}
for (; j < M; j++) { PADJUMP()
u32 x0 = (j <= len) ? src[j] : 0;
bp[j] = x0;
u32 w1 = g_t5[0][j], w2 = mulmod(w1, w1), w3 = mulmod(w2, w1), w4 = mulmod(w3, w1);
bp[j + (size_t)g_M5] = shoup1s(x0, w1, shoup_c(w1));
bp[j + (size_t)2 * g_M5] = shoup1s(x0, w2, shoup_c(w2));
bp[j + (size_t)3 * g_M5] = shoup1s(x0, w3, shoup_c(w3));
bp[j + (size_t)4 * g_M5] = shoup1s(x0, w4, shoup_c(w4));
}
}
TGT static void merge5_dit(u32 *a, u32 M) {
u32 *bp = a; u32 nxt0 = 1u << 20; u32 nxt1 = 1u << 18; u32 nxt2 = 1u << 16; u32 nxt3 = 1u << 14;
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 cK = _mm256_set1_epi32((int)g_c5k), cKs = _mm256_set1_epi32((int)g_c5k_s);
const __m256i cDiffI = _mm256_set1_epi32((int)g_c5dci), cDiffIs = _mm256_set1_epi32((int)g_c5dci_s);
const __m256i cSumI = _mm256_set1_epi32((int)g_c5sci), cSumIs = _mm256_set1_epi32((int)g_c5sci_s);
const __m256i iv2 = _mm256_set1_epi32((int)g_inv2), iv2s = _mm256_set1_epi32((int)g_inv2_s);
u32 s = 0;
for (; s + 16 <= M; s += 16) { PADJUMP_S()
pf5(bp + s + PF_DIST); pf5(bp + s + (size_t)g_M5 + PF_DIST); pf5(bp + s + (size_t)2 * g_M5 + PF_DIST);
pf5(bp + s + (size_t)3 * g_M5 + PF_DIST); pf5(bp + s + (size_t)4 * g_M5 + 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 ta[5], tb[5];
ta[0] = ldu256((bp + s)); tb[0] = ldu256((bp + s + 8));
{
__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(w2, w2, der_ws(w2), pv5);
__m256i w1b = ldu256((g_u5[0] + s + 8));
__m256i s1b = der_ws(w1b);
__m256i w2b = vshoup(w1b, w1b, s1b, pv5), w3b = vshoup(w2b, w1b, s1b, pv5), w4b = vshoup(w2b, w2b, der_ws(w2b), pv5);
ta[1] = vshoup(ldu256((bp + s + (size_t)g_M5)), w1, s1, pv5);
tb[1] = vshoup(ldu256((bp + s + 8 + (size_t)g_M5)), w1b, s1b, pv5);
ta[2] = vshoup(ldu256((bp + s + (size_t)2 * g_M5)), w2, der_ws(w2), pv5);
tb[2] = vshoup(ldu256((bp + s + 8 + (size_t)2 * g_M5)), w2b, der_ws(w2b), pv5);
ta[3] = vshoup(ldu256((bp + s + (size_t)3 * g_M5)), w3, der_ws(w3), pv5);
tb[3] = vshoup(ldu256((bp + s + 8 + (size_t)3 * g_M5)), w3b, der_ws(w3b), pv5);
ta[4] = vshoup(ldu256((bp + s + (size_t)4 * g_M5)), w4, der_ws(w4), pv5);
tb[4] = vshoup(ldu256((bp + s + 8 + (size_t)4 * g_M5)), w4b, der_ws(w4b), pv5);
}
__m256i Y0, Y1, Y2, Y3, Y4, Z0, Z1, Z2, Z3, Z4;
dft5v(ta[0], ta[1], ta[2], ta[3], ta[4], Y0, Y1, Y2, Y3, Y4, pv5, cK, cKs, cCi, cCis, cDiffI, cDiffIs, cSumI, cSumIs, iv2, iv2s);
dft5v(tb[0], tb[1], tb[2], tb[3], tb[4], Z0, Z1, Z2, Z3, Z4, pv5, cK, cKs, cCi, cCis, cDiffI, cDiffIs, cSumI, cSumIs, iv2, iv2s);
if (g_fuse) {
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);
Z0 = vshoup(Z0, nv, ns, pv5); Z1 = vshoup(Z1, nv, ns, pv5); Z2 = vshoup(Z2, nv, ns, pv5);
Z3 = vshoup(Z3, nv, ns, pv5); Z4 = vshoup(Z4, nv, ns, pv5);
u32 *C = g_outC + s;
if (g_nt) {
_mm256_stream_si256((__m256i *)(C), Y0); _mm256_stream_si256((__m256i *)(C + 8), Z0);
_mm256_stream_si256((__m256i *)(C + M), Y1); _mm256_stream_si256((__m256i *)(C + M + 8), Z1);
_mm256_stream_si256((__m256i *)(C + 2 * M), Y2); _mm256_stream_si256((__m256i *)(C + 2 * M + 8), Z2);
_mm256_stream_si256((__m256i *)(C + 3 * M), Y3); _mm256_stream_si256((__m256i *)(C + 3 * M + 8), Z3);
if (s < g_limS - 4 * M) { _mm256_stream_si256((__m256i *)(C + 4 * M), Y4); _mm256_stream_si256((__m256i *)(C + 4 * M + 8), Z4); }
} else {
stu256((C), Y0); stu256((C + 8), Z0);
stu256((C + M), Y1); stu256((C + M + 8), Z1);
stu256((C + 2 * M), Y2); stu256((C + 2 * M + 8), Z2);
stu256((C + 3 * M), Y3); stu256((C + 3 * M + 8), Z3);
if (s < g_limS - 4 * M) { stu256((C + 4 * M), Y4); stu256((C + 4 * M + 8), Z4); }
}
} else {
stu256((bp + s), Y0); stu256((bp + s + 8), Z0);
stu256((bp + s + (size_t)g_M5), Y1); stu256((bp + s + 8 + (size_t)g_M5), Z1);
stu256((bp + s + (size_t)2 * g_M5), Y2); stu256((bp + s + 8 + (size_t)2 * g_M5), Z2);
stu256((bp + s + (size_t)3 * g_M5), Y3); stu256((bp + s + 8 + (size_t)3 * g_M5), Z3);
stu256((bp + s + (size_t)4 * g_M5), Y4); stu256((bp + s + 8 + (size_t)4 * g_M5), Z4);
}
}
for (; s + 8 <= M; s += 8) { PADJUMP_S()
pf5(bp + s + PF_DIST); pf5(bp + s + (size_t)g_M5 + PF_DIST); pf5(bp + s + (size_t)2 * g_M5 + PF_DIST);
pf5(bp + s + (size_t)3 * g_M5 + PF_DIST); pf5(bp + s + (size_t)4 * g_M5 + 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((bp + 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(w2, w2, der_ws(w2), pv5);
t[1] = vshoup(ldu256((bp + s + (size_t)g_M5)), w1, s1, pv5);
t[2] = vshoup(ldu256((bp + s + (size_t)2 * g_M5)), w2, der_ws(w2), pv5);
t[3] = vshoup(ldu256((bp + s + (size_t)3 * g_M5)), w3, der_ws(w3), pv5);
t[4] = vshoup(ldu256((bp + s + (size_t)4 * g_M5)), 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, cK, cKs, cCi, cCis, cDiffI, cDiffIs, cSumI, cSumIs, iv2, iv2s);
if (g_fuse) {
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((bp + s), Y0);
stu256((bp + s + g_M5), Y1);
stu256((bp + s + 2 * g_M5), Y2);
stu256((bp + s + 3 * g_M5), Y3);
stu256((bp + s + 4 * g_M5), Y4);
}
}
for (; s < M; s++) { PADJUMP_S()
u32 t[5]; t[0] = bp[s];
{ u32 w1 = g_u5[0][s], w2 = mulmod(w1, w1), w3 = mulmod(w2, w1), w4 = mulmod(w3, w1);
t[1] = shoup1s(bp[s + (size_t)g_M5], w1, shoup_c(w1));
t[2] = shoup1s(bp[s + (size_t)2 * g_M5], w2, shoup_c(w2));
t[3] = shoup1s(bp[s + (size_t)3 * g_M5], w3, shoup_c(w3));
t[4] = shoup1s(bp[s + (size_t)4 * g_M5], 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 {
bp[s] = Y[0];
bp[s + g_M5] = Y[1]; bp[s + 2 * g_M5] = Y[2]; bp[s + 3 * g_M5] = Y[3]; bp[s + 4 * g_M5] = Y[4];
}
}
}
static void ntt5_fwd(u32 *a, const u32 *src, u32 len) {
stage5_dif_src(a, src, len, g_M5L);
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_M5L);
}
static u32 g_dataA[5 * ((1u << 22) + 448 * TOPPAD + 896)] __attribute__((aligned(64)));
static u32 g_dataB[5 * ((1u << 22) + 448 * TOPPAD + 896)] __attribute__((aligned(64)));
static int g_ready = 0;
__attribute__((target("avx2"))) void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c) {
const u32 M = 1u << 22;
const u32 N = 5u * M;
if (!g_ready) { ntt_init(M); init_r5(N); _mm_sfence(); 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 = powmod(N % NTP, NTP - 2);
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[padpos(i)] * ninv % NTP); i++; }
for (; i + 8 <= lim; i += 8) {
__m256i x = vshoup(ldu256((A + padpos(i))), nv, ns, pv);
_mm256_stream_si256((__m256i *)(c + i), x);
}
for (; i < lim; i++) c[i] = (u32)((u64)A[padpos(i)] * ninv % NTP);
_mm_sfence();
}
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 367.784 ms | 302 MB + 392 KB | Accepted | Score: 100 | 显示更多 |