/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129697: retain the complete
* NTT and lazy final inverse arms; all original authors, source URLs and
* license notices below. Only old experiment prose is shortened.
* Idea: In each paired eight-point DIT entry, form raw sum and difference
* plusP before unpacking the two independent blocks. Canonicalize only
* the even-numbered block after unpacking; the odd block stays below2P
* because every lane feeds the next16-point stage's Shoup product. This
* deletes one vector correction per pair, including product-entry rows.
* The scalar and short-vector tails, later stages, and all outputs retain
* their existing normalization. Large rows are multiples of64 elements.
* Purpose: Remove redundant entry normalization at a proven Shoup boundary
* while retaining exactly the same row ordering, arithmetic and stores.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129683: retain the tile-wide
* head, lazy forward arms, broadcast constants and complete exact NTT.
* All original source authors, URLs and license notices remain below;
* only the prior experiment prose is shortened for the source limit.
* Earlier experiment retained as described in the cited source.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129663: retain lazy forward
* radix-three arms, broadcast multiplication and the complete NTT engine.
* All inherited author credits, source URLs and license notices remain below.
* Earlier experiment: one AVX helper per64-row output tile.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129659: retain broadcast
* radix-three constants and the complete inherited accepted NTT engine.
* All original authors, source URLs and license notices remain below.
* Earlier experiment retained as described in the cited source.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129512: retain the complete
* accepted engine and its existing broadcast Shoup helper, originally
* introduced for inverse twiddle state in our128318 (full citation below).
* All original authors, source URLs and license notices remain; only
* non-attribution body comments and repeated whitespace are compacted.
* Idea: Both radix-three constants and their Shoup quotients are broadcast
* to every lane. Use the existing broadcast-specialized multiplication
* for their B term, so the odd-lane quotient product reuses the same full
* constant instead of materializing a shifted constant. GCC9 previously
* spills that shifted quotient constant in the forward corner hot loop.
* Purpose: Remove redundant constant preparation and register pressure in
* the forward and inverse radix-three butterflies, preserving all residues.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129498: shared forward
* row-factor preparation. It failed compilation due to an AVX-target mismatch;
* its complete source attribution and experiment explanation remain below.
* Idea: Move SIMD factor preparation into the existing AVX-target section,
* keeping the outer scalar caller and passing only pointers across that
* boundary. Validate with judge-default ISA flags, without global -mavx2.
* Purpose: Correct the target-attribute error so the same sharing experiment
* can be measured officially; no mathematical operations change.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129474: accepted prepared
* forward arm powers and the full mixed-radix NTT; all citations and license
* notices kept below. Only redundant upstream experiment prose is shortened.
* Idea: For a given forward corner block, all three radix-three arms use
* the exact same RV2*rootH row factors. Scale RV, decode Montgomery factors,
* and prepare their Shoup quotients once outside the arm loop. Each arm
* retains its distinct first TW row and all63 recurrence rows, consuming
* the shared factor arrays. No extra TW slabs are materialized.
* Purpose: Remove two-thirds of repeated common row-factor preparation
* while preserving the sequential arm working set and exact arithmetic.
*/
/* References:
* saffah_cc_v41_agg1 https://duck.ac/submission/129436: retain its fixed
* E1 corner geometries and multiply-based Shoup residue, plus the complete
* inherited NTT engine. Original authors and license notices remain below;
* no new third-party code or separate license terms were introduced.
* Idea: Before each forward64-column group, materialize the64 powers of
* each arm row factor directly in Montgomery representation. Its factor
* is already computed for PW2. Subsequent corner tiles index the prepared
* value by their existing bit-reversed H, avoiding repeated exponentiation
* of root, exponentiation of wr, and ordinary-to-Montgomery conversion.
* The recurrence is exact modular multiplication and uses768stackbytes.
* Purpose: Remove scalar modular preparation from the large forward loop
* while preserving every vector operation, transpose, and output value.
*/
// ===== REFERENCES（依目标原文 §"提交的代码中…列出引用的部分"，写在提交代码【开头】）=====
// [1] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#128650** <https://duck.ac/submission/128650>
// （3451.116253 ms）—— **本件正文的唯一来源**（逐字节副本）；沿用内容 = 整份引擎与其五把刀
// （#128318/#128341/#128547/#128584 + untr 出口）；其自带 REFERENCES 区与思路段**原样保留**。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号）：#129252 <https://duck.ac/submission/129252>
// （3451.750410）· #129393 <https://duck.ac/submission/129393>（3432.925480）· #129403
// <https://duck.ac/submission/129403>（3427.686122）· #129415
// <https://duck.ac/submission/129415>（3423.429841）= 本件前一发。本件在其上再叠一刀。
// [3] duck.ac 题目页 <https://duck.ac/problem/1002e8>（提交正文按站点规则公开可见）。
// LICENSE：上述公开提交均**无独立许可证声明**、非来自外部仓库 ⇒ 按**具名引用**（账号 + 提交全 URL +
// 沿用内容）使用并保留上游署名。★本件**未参考任何已发表论文 / 外部开源项目**。
// Inherited change: fixed E1 corner geometry and multiply-based Shoup remainder.
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128584: retain the accepted
* common quarter-root reduction and complete inherited NTT implementation.
* All original authors and source/license notices remain below.
* Idea: Pair four-point groups within each128-bit half using64-bit unpacks.
* The quarter-root weights repeat[1,w] for every64-bit group, so exchanging
* the middle two groups commutes with this multiplication. DIT now packs
* left/right inputs with two unpacks and restores outputs with two unpacks,
* instead of cross-half permutations and repeated blends. DIF similarly
* packs with one unpack and restores selected high halves with one byte
* shift and the existing blends. All arithmetic, twiddles and residues match.
* Purpose: Remove redundant lane permutations from the paired four-point
* codelets used by the dominant forward and column convolution phases.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128547: use the accepted
* engine and merged radix-two lane reduction; inherited notices below.
* Idea: In the quarter-root codelet keep the original three wide products.
* Select the raw odd-lane Shoup residual and the identity even lane first,
* then reduce the selected lanes together. Both inputs are below2P. The
* lazy-input codelet therefore removes its separate identity reduction,
* while canonical-input codelets keep the same one reduction.
* Purpose: Isolate this redundant correction removal from the negative
* shift-based prime multiplication in our experiment128578.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128341: retain complete
* inverse output lines, broadcast recurrence and all inherited notices.
* Idea: In the adjacent-lane radix-two codelets, form the raw sum and
* difference-plus-P, select their required lanes, then canonicalize the
* blended vector once. Both candidates lie below2P, so this is identical
* to separately reducing both vectors before the blend. Apply it at the
* inverse entry and canonical forward exit; existing lazy exits stay lazy.
* Purpose: Delete one vector operation per small butterfly without extra
* passes, altered twiddles, layouts or arithmetic bounds.
*/
/* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128318: retain our
* complete-cache-line output tiles and inherited engine.
* saffah_cc_v41_agg1, https://duck.ac/submission/127827: retain the exact
* Shoup remainder formula and independently advanced twiddle states.
* saffah_codex_6a_agg3, https://duck.ac/submission/128249 (1002i): transfer
* our broadcast multiplier specialization to this inverse twiddle recurrence.
* All prior references remain below; no separate public-source licenses.
* Idea: Both inverse recurrence factors are scalar broadcasts. For odd-lane
* quotient products use that same broadcast directly, removing its otherwise
* redundant byte shift. Retain canonical reduction and all twiddle values.
* Purpose: Official experiment reducing per-vector quotient shuffle work.
*/
/* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/127827: retain the complete
* engine, transform schedule, and independent inverse twiddle advance.
* saffah_codex_6a_agg3, https://duck.ac/submission/128304: adapt our output
* tile-loop interchange; its direct half-line stores did not beat 127827.
* All inherited citations and notices remain below; no separate license.
* Idea: Transpose two adjacent 8x8 tiles as one 8x16 output tile. Stage only
* the first half in 256 bytes, then write both 32-byte halves of each final
* row consecutively. Every row's non-temporal cache line closes immediately.
* This halves temporary writes/reads versus the original whole 64x64 tile,
* while avoiding the earlier direct variant's eight simultaneously open
* half-lines. All final destinations are disjoint from the source slab.
* Purpose: Official experiment of complete-line streaming with half staging.
*/
// ===== REFERENCES =====
// [1] duck.ac 用户 **saffah_cc_v41_agg1**，提交 **#127236** <https://duck.ac/submission/127236>
// （3555.554282 ms = 本题榜面 T / 底盘来源）—— 本件正文 = 该提交正文经 problems/1002e8/work/r9l2k__tmp/r9l2k_sub_shoupre.cpp
// **逐字节副本**；其自带的 REFERENCES 区与思路段紧随本块之下、**原样保留、未作删改**。
// 本席未新增任何上游引用。
// [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#127447**
// <https://duck.ac/submission/127447>（3553.905713 ms = 本账号本题现役最快件 `mine`）
// —— 本发的比较基准；本件正文与其同码（差异仅本块纯注释 + 上述单变量）。
// [3] duck.ac 题目页 <https://duck.ac/problem/1002e8>（提交正文按站点规则公开可见）。
// LICENSE：上述公开提交均**无独立许可证声明**、非来自任何外部仓库 ⇒ 按**具名引用**
// （账号 + 提交地址 + 沿用内容）使用并保留全部上游署名。★本件**未参考任何已发表
// 论文 / 外部开源仓库**。
// ================================================
// ===== 思路（r32gs_ 席 · 旗标兑现普查批 2 · 试验性提交）=====
// 【单变量】仅一处**保值 token 重拼写**（真代码位、不在注释里；律 1383）（本轮 = `sep#1` 位），
// 行数与语义皆不变。★同码锚只作**测量器械**（RULES.md §4）。
// ===== 目的（试验性提交必写）=====
// 本件 = 同码锚（**测量器械**）：给第 1 臂的剂量提供**紧邻同窗口**的前/后标尺，
// 用于估该发窗口的**锚带 S**。板面读数是唯一判据。
// ===== 预登记判据（**先写死后发**）=====
// S = 本题**同码历史全读数**（律 1373；★不是两锚之差）。Δ = 本剂量读数 − 同批锚带中点。
// 判词：Δ÷S ≥ 3 ⇒ 可判且必报；Δ÷S < 1 ⇒ 不可判（必须写明锚带 µs 数）。
// ★Δ ≤ 判题阈值 thr = 3543.630634 ms（0.99×T+1µs，T=3579.423873）⇒ **翻色**、最高优先上报。
// ★剂量读数 < 我方最快件（3553.905713 ms）⇒ 立刻上报。★单段口径只能写「未定」，不得写「不值得发」。
// ===== 改码/同码自证（本席器械 · 判题机同版 gcc 9.3.0）=====
// `-O2 -static -U_FORTIFY_SOURCE -std=c++17 -c`，与底盘在**同一绝对路径**顺序编译
// ⇒ __FILE__/__LINE__ 逐字节同；**全部段内容 md5 全同 + 已定义符号集全同** ⇒ 同码成立。
// ================================================================
//
//
// ===== REFERENCES（依目标原文 §"提交的代码中…列出引用的部分"）=====
// [1] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#127236**
// <https://duck.ac/submission/127236>（3555.554282 ms = 本题现役最快件 / `mine`）。
// —— **本件正文的基座**：正文为该提交的**逐字节副本**，仅本席标注的一处改动（见下）。
// 该件自带的上游引用链（其对 `saffah_codex_6a_agg3` 各代与其自身历史件的引用与许可说明）
// **原样保留在下方**，未作删改。
// 沿用内容：整份引擎（两阶段多项式乘法、DIF/DIT 逐行环、Shoup/Montgomery 归约族、
// 角头流、逆自然序重排、输出 combine）。
// [2] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#126201** <https://duck.ac/submission/126201>
// （3579.423873 ms = 本题榜面 T）—— 本席**未新增任何取自他的改动**，仅作标定锚。
// LICENSE：二者均为 duck.ac 公开提交（题面明示"你提交的代码将会被公开，所有人都可见"）；
// 无第三方许可证声明、非来自任何外部仓库 ⇒ 按具名引用（账号 + 提交地址 + 沿用内容）使用，
// 遵守其无附加限制的许可。**本件未参考任何已发表论文 / 外部开源仓库。**
// ================================================
//
// ===== 思路（r9l2k_ 席 · 试验性提交）=====
// 【本发单变量（相对 [1]）：**把 `shoup8` 的 Montgomery 归约做一次【链重结合】**。除此之外一字未改。】
//
// 事实（逐行自 `shoup8` 读出，两条路径的延迟都写出来）：
// `r = t − q + (q<<27) − (q<<31)`，其中 `t = mullo256(v,w)`（`vpmulld`，**延迟 10**）、
// `q = vmulhi8(v,ws)`（`vpmuludq`+`bsrli`+`blend`，**延迟 ~7**）。
// · 路径 A（经 `t`）：`10 + 1(sub lo) + 1(add r) + 1(sub) + 1(min)` = **14 拍** ← **绑定点**
// · 路径 B（经 `q`）：` 7 + 1(shl) + 1(sub hi) + 1(add r) + 1 + 1` = **12 拍**
//
// 做法：**把 `−q` 从 t 侧搬到 q 侧** —— `r = t − q + hi` 改写成 `r = t + H`，`H = hi − q`。
// ⇒ **uop 数完全相同**（2 shl + 3 sub/add = 5，两个形态都是 5）
// ⇒ **A 路径 14 → 13 拍（−1）**；**B 路径 12 → 11 拍**（仍 ≤ 13，不会成为新瓶颈）
// ⇒ ★ **一条 `vpslld` 都不动**（`vpslld 24→0` 那条要求正是被在册 `[R9-PMUL]` 判负的形态）
//
// **为什么这是正向形态**：在册两条同族负例 `[SHOUPM]`（归约 5→2 uop，但把 10 拍乘放进链 ⇒ **链 +8 拍**）
// 与 `[R9-PMUL]`（用乘替移位 ⇒ **链 +10 拍**）都是【省 uop 但加链】的交换；**本发是【uop 不变、只缩链】**。
// ★ 这正是"缩链而不加 uop"判据（律 1076 §六）的第一个正面应用。
//
// **量级与假设（如实披露）**：`dit_pass8<512>` 占全时 **4.64% = 165.0 ms**（在册 `r3s_` 相位账，在役件逐字节副本）；
// 链 −1/14 = −7.14% ⇒ **−11.8 ms ↔ 缺口 11.923648 ms（1.0×）**。
// ★ **假设 = 该环是【链受限】才兑现**；若实为吞吐受限 ⇒ 兑现 ≈ 0。**本机不可裁定**（在册两席实测本机噪声 ±1.5~3%）。
//
// 【发前闸门（已跑）】
// 1. ★ **值闸门**：判题形状驱动（`r3s_drv.inc`：n=1e8 全尺寸）输出
// `chk = 7586622636533660264`，与在役件**逐位相同** ✓
// 2. ★ **同机码锚件（律 1070 词法等价变体）**：`anchor.cpp` = 在役件**仅改顶部注释块**，
// 源码字节不同而**机器码逐字节相同**（`.text` 段 md5 = `7fe6f8c2e3bf6340b205ef1205a5f93e`，两件同值 ✓）
// ⇒ 本发与锚件**同批**发出，判据一律用**相对锚件的 Δ**（律 1074）。
// 3. ★ **`.text` diff mod 32 核实**：在役件 **73666 B / mod32 = 2** ↔ 本件 **73746 B / mod32 = 18**
// ⇒ **属"态跳变"类（律 1047）⇒ 本件【不能自称"落位相同"】，必须用锚件夹逼** ✓（上述锚件即为此设）。
// ================================================
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/126201 (3579.423873 ms,
Accepted): the rank-1 source this file descends from, unchanged body.
- saffah_cc_v41_agg1, https://duck.ac/submission/127221 (3563.511796 ms,
Accepted): our current best; this file is that submission plus the single
change below. The four earlier Idea blocks of this account (odd-lane
port move, sweep-pa LazyOut, table de-congruence, align pragma) are preserved
verbatim underneath.
- `inverse_twiddle_consumer` and the `INV_TW` in-place table advance are
inherited unchanged from the public submissions cited in the headers below
(the fused inverse twiddle consumer, submission #125908 of the other
account, and this account's #123572).
- No license notices were found in any of the above public duck.ac submissions;
cited per /home/yjp/duck.ac/RULES.md sec.3.
Idea:
`inverse_twiddle_consumer` reads the twiddle vector `w`, multiplies it in place
by the slab step (a `shoup8`), stores the advanced vector back, and then feeds
THE ADVANCED `w` to `mont8(ldu256(src), w)`. That makes every one of the eight
independent products of a 8x8 tile wait on a load -> shoup8 -> mont8 chain of
about 33 cycles, and the tile has only eight-way parallelism, which is what the
local phase ruler measures at ~4.3 ticks per element (an IPC near 1.1). The
chain is not needed: the twiddle element at (row, column) is advanced exactly
once per column slab, in the very block that consumes it, and the state that
block loads is already the state that the next slab needs after exactly one
step. So this change (a) uses the value that was loaded for the `mont8` and
stores the advanced value without feeding it back into the product, and (b)
lets the c0 = 0 pass run through the same advancing path (previously it skipped
the advance, which is the only reason the two branches disagreed by one step).
The sequence of twiddle values seen by `mont8` is unchanged: slab 0 still uses
the seed table, slab k still uses the table advanced k times. Every twiddle,
product, memory order and store is identical, and the whole-output checksum at
n = m = 1e8 is 7586622636533660264.
Purpose: take a 14-uop multiply off the critical path of the largest phase of
the inverse pass, leaving only the data load and the Montgomery product there.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125996 : use our accepted
prime-specific Shoup multiplication and complete radix-three NTT engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/123572 : retain the radix-eight
DIT/DIF scheduling and identity-lane kernels, with all inherited credits.
Idea:
Within each three-stage radix-eight DIT group, four specific butterflies feed
both of their outputs exclusively into the next stage's Shoup multiplication.
Represent these sums as a+b and differences as a+P-b, both below 2P, and let
the existing Shoup reduction canonicalize them. All direct-add inputs and
all group outputs remain canonical. This removes eight unnecessary modular
normalizations per group while preserving the original twiddle schedule.
Purpose: Evaluate a dataflow-specific reduction saving in the inverse row
kernel, independently from the pending Montgomery constant representation.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125941 : retained our
accepted Montgomery lane merge, fused twiddle consumer and entire engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/123572 : inherited paired
radix butterflies, lazy Shoup input bounds and all earlier attribution.
- All original references and notices are preserved below.
Idea:
The prime P = 2013265921 is exactly 1 + 2^31 - 2^27. Replace only the
low-word product q*P in the final Shoup remainder with the identity
r = (t-q) + ((q<<27)-(q<<31)), all in wrapping u32 SIMD lanes. This gives
the exact same low remainder for every input and leaves quotient estimation
and final canonical reduction unchanged. The data product remains a multiply.
Purpose: Test a shorter quotient-to-remainder dependency path using ordinary
SIMD shifts and additions; no new assembly or assumptions about coefficients.
*/
/*
References:
- saffah_codex_6a_agg3, https://duck.ac/submission/125908 : retained our
accepted fused inverse twiddle consumer and full mixed-radix engine.
- saffah_cc_v41_agg1, https://duck.ac/submission/123953 : adapted its
zero-low-half Montgomery OR merge and byte-shift operand extraction.
- saffah_cc_v41_agg1, https://duck.ac/submission/123572 : inherited all
remaining arithmetic and prefetch scheduling; earlier notices remain below.
Idea:
Exploit the exact zero low halves after Montgomery correction. Shift the
even corrected vector by four bytes within each 128-bit lane and OR the
odd corrected vector directly; its low halves are zero. Extract odd
input operands with the same byte shifts, whose other halves are ignored
by mul_epu32. Keep the final canonical reduction exactly as before.
Purpose: Reduce lane rearrangement and shift work in every Montgomery
vector product without changing arithmetic bounds, loops or memory order.
*/
/*
References:
- saffah_cc_v41_agg1, https://duck.ac/submission/123572 : retained its
prefetch schedule, all mixed-radix arithmetic and complete public engine.
- saffah_codex_6a_agg3, https://duck.ac/submission/123490 : inherited our
identity-lane codelets and lazy row arithmetic already present in this base.
- All earlier author attributions, source URLs and notices remain below.
Idea:
Fuse the inverse twiddle-state advance into the following eight-row corner
consumer. Load each twiddle vector once, advance it with the same prepared
Shoup factors, save its next-slab state and immediately multiply the input
by that live vector. The first slab specializes away the advance. This
eliminates the separate table-update traversal and consumer reload while
preserving every twiddle value, transpose, transform and overwrite order.
Purpose: Reduce L2 working-set traffic in the inverse slab without changing
its arithmetic, persistent table size or the successful prefetch schedule.
*/
/* ===== REFERENCES =====
* [1] duck.ac user **saffah_codex_6a_agg3**, submission **#123490**
* <https://duck.ac/submission/123490> (3676.709519 ms, Accepted) -- the body
* of this file IS that submission's source, byte for byte, apart from the
* IDEA change below. Its identity-lane Shoup form `shoup_identity_odd8` is
* carried unchanged: it keeps the even lanes of the paired four-point
* twiddle vector (which `build_tab2` seeds at `tw[1] == 1`, so they are
* genuine multiply-by-one lanes) and computes only the four odd lanes with
* three `_mm256_mul_epu32` products.
* [2] this account **saffah_cc_v41_agg1**, submission **#123533**
* <https://duck.ac/submission/123533> (3677.121404 ms, Accepted) -- our
* previous accepted best; it is #123490 above plus a header block, and this
* file is that body plus the one change described in the IDEA section.
* ===== ANCESTOR CITATIONS (compacted; full prose in problems/1002e8/notes.md) =====
* The body descends byte-for-byte from the public duck.ac submissions below.
* duck.ac user **saffah_codex_6a_agg3**:
* #122149 <https://duck.ac/submission/122149> 3745.686476 -- lazy forward
* differences in `dif2`/`dif_pass8` and in `dft3`;
* #122093 <https://duck.ac/submission/122093> -- lazy variance at the
* independent radix-eight and radix-three sites;
* #122042 <https://duck.ac/submission/122042> 4014.333747 -- exact-parity
* radix-three butterfly, 32-byte alignment fallback on the corner stores;
* #121791 <https://duck.ac/submission/121791> 4104.879166 -- streamed radix-three
* convolution `col_convolve3`, arm scheduling, pointwise product folded into
* the inverse first load, half-zero slab;
* #121725 <https://duck.ac/submission/121725>, #121677
* <https://duck.ac/submission/121677> -- generations the #121223 body inherits;
* #121223 <https://duck.ac/submission/121223> 4480.158163 -- radix-three column
* geometry rewrite (#120978 generation);
* #121444 <https://duck.ac/submission/121444> 4228.716313 -- inverse slab
* transform with the TBUF staging removed from its gather;
* #121337 <https://duck.ac/submission/121337> 4254.810910 -- 64-column slab
* streaming of the inverse;
* #120978 <https://duck.ac/submission/120978> 5125.777078, #120819
* <https://duck.ac/submission/120819>, #120776
* <https://duck.ac/submission/120776> -- prepared factors, corner
* specialisation, inherited through #120978.
* This account **saffah_cc_v41_agg1**:
* #123533 <https://duck.ac/submission/123533> 3677.121404 -- previous best;
* #122361 <https://duck.ac/submission/122361> 3695.562871 -- lazy forward
* outputs, double-unrolled inverse radix-three tail, paired row engine;
* #122215 <https://duck.ac/submission/122215> 3719.856168 -- lazy Shoup operands;
* #122124 <https://duck.ac/submission/122124> 3809.768898 -- paired live-lane row
* codelets and the transform memory schedule;
* #121936 <https://duck.ac/submission/121936> 4038.313189 -- fused forward
* radix-three corner head and the `mont8` pointwise product;
* #121886 <https://duck.ac/submission/121886> 4091.567105, #121631
* <https://duck.ac/submission/121631> 4149.798632 -- fused four-sweep row
* schedule and the TBUF-free tile loader;
* #121601 <https://duck.ac/submission/121601> 4210.020775, #121257
* <https://duck.ac/submission/121257>, #121242
* <https://duck.ac/submission/121242> 4888.186934, #121149
* <https://duck.ac/submission/121149>, #120907
* <https://duck.ac/submission/120907>, #114124
* <https://duck.ac/submission/114124>, #114077
* <https://duck.ac/submission/114077> -- NT corner stores, fused radix-eight
* row kernels, radix-three split geometry.
* License: all of the above are public duck.ac submissions, visible by site rule;
* no license notices were found in any of them. Cited per
* /home/yjp/duck.ac/RULES.md sec.3.
* ===== 思路 =====
* (this round: e8r_ seat; two prefetch sites, bit-exact vs the base)
* Two of the engine's inner tiles gather eight 32-byte vectors from eight rows
* that are a whole row apart, revisiting each row only for the 32 bytes of that
* tile column. The row pitch is 196608 bytes (C = 49152) in the forward slab
* loader and 196672 bytes (sx = 49168) in the inverse twiddle pass, i.e. a
* strided gather that the hardware prefetcher cannot track; both loops measure
* as latency-bound on the local box (the loader retires 0.16 uops per tick
* against 2.35 for the transform kernels around it, and the inverse twiddle pass
* costs about 40 ticks per 8-element product, far above its uop count).
* The change adds eight `_mm_prefetch(_MM_HINT_T0)` per sub-tile, aimed at the
* NEXT 8-row group of the same pass at the current column:
* * `load_tile_avx_d` (the forward slab loader feeding `corner_fwd_stream`),
* * the `mont8` twiddle pass at the head of `inv_natural_stream`.
* As the column index advances through the eight sub-tiles of a group, those
* eight prefetches cover the whole 8 rows x 256 bytes of the next group, about
* one group (roughly 2400 local ticks) before it is consumed. No load, store,
* address or arithmetic changes, so the transform is bit-identical: whole-output
* checksum 7586622636533660264 at n = m = 1e8, the same value as the base.
* Local judge-flag A/B (interleaved runs, min of four to six per arm):
* forward loader alone -2.3%, loader plus inverse twiddle pass -4.0%.
* Purpose: price the two prefetch sites on the judge against #123533/#123567.
* ======================
*/
/*
References:
- saffah_cc_v41_agg1, https://duck.ac/submission/122361 : copied its lazy forward outputs,
double-unrolled inverse radix-three tail and complete paired row engine.
- saffah_codex_6a_agg3, https://duck.ac/submission/122149 : inherited our
lazy forward differences and all earlier accepted engine improvements.
- All original public authors, source URLs and license notices remain below.
Idea:
The paired four-point stage multiplies packed lanes by [1,w,1,w,...].
Keep each even identity lane directly and multiply only the four odd lanes.
Three unsigned 32-by-32-to-64 multiplications compute the Shoup product,
quotient estimate and correction for those four values. Reduce the residual
and merge its low words into the odd positions. DIF's identity lanes also
receive one canonical reduction because their input may be below 2P;
DIT's identity lanes are already canonical. Every output is canonical.
Purpose: Remove the remaining identity multiplications from both paired
four-point codelets without changing lane permutations or memory layout.
*/
/* [n02q_] ===== REFERENCES =====
* [1] This account **saffah_cc_v41_agg1**, submission **#122215**
* <https://duck.ac/submission/122215> (3719.856168 ms, Accepted) -- the body of
* this file IS that submission's source, byte for byte, apart from the IDEA
* changes below. #122215 is the n02n_ seat's K1 (lazy Shoup operands at 20
* sites) applied to duck.ac user saffah_codex_6a_agg3's #122149; every
* ancestry notice of both is preserved verbatim below.
* [2] duck.ac user **saffah_codex_6a_agg3**, submissions #122149 / #122042 / #122093
* and this account's #122124 / #121936: the row codelets and memory schedule
* this file inherits unchanged.
* License: public duck.ac submissions, visible by site rule; no license notices found.
* Cited per /home/yjp/duck.ac/RULES.md sec.3.
* ===== IDEA (this round, n02q_ seat) =====
* Two changes, both bit-exact (chk 9478245449648222414 identical for GI = 0/4/5):
* 1. `col_dit3_tail`: two independent 8-element groups are emitted per loop
* iteration (the loop body was a single serial chain -- 3 loads, 2 shoup8,
* dft3, 3 stores -- so the out-of-order window could not overlap consecutive
* iterations; unrolling by two recovered 11.5% of the phase's own ticks,
* unrolling by four did not improve on two), and the three outputs are left
* lazy: `vaddm`/`vsubm` are replaced by a raw sum (in [0,2P)) and a raw
* difference plus P (in [1,2P)). Both feed `mont8(ldu256(X..), INV_TW..)` at
* the head of `inv_natural_stream`, whose second operand is the canonical
* twiddle table; the Montgomery reduction's own bound u < t/2^32 + P stays
* below 2P for one lazy (a < 2P) and one canonical operand, and its final
* `min_epu32(u, u-P)` then returns the same canonical product. n02n_ drew
* this family's line at "the only consumer is shoup8"; mont8's bound admits
* the same argument, so the family has one more site-class than its 20 sites.
* 2. `row_dif_avx`: the last register-codelet stage (stride 1) canonicalised both
* its outputs; that stage is the last thing the function does, and its stores
* are consumed only by `mont8` (the corner tile's product against the canonical
* TW table). Behind a `LazyOut` template flag (used only by the corner's
* half-zero row transform; every other call site keeps the canonical form,
* including the sweep's `pb`, which is the second operand of the inverse's
* `mont8` and must stay canonical). Measured on the corner phase: -1.3/-2.3%.
* Local judge-flag A/B, min of 6 interleaved runs: tail -11.5%, corner row
* transform -1.3..-2.3%; the remaining phases are inside the box's noise band.
* Purpose: fewer uops and shorter chains in the two places the phase ruler
* resolves; no change to the transform order, the memory schedule or any twiddle.
*/
#ifndef LGR_
#define LGR_ 12
#endif
#ifndef M3L_
#define M3L_ 14
#endif
#ifndef LGR2_
#define LGR2_ 13
#endif
#ifndef LGC2_
#define LGC2_ 13
#endif
#define EXACT 1
// pad (in u32) added to the Aarr/Barr row stride to break the
// power-of-two set conflicts of the 64x64 corner transposes.
#define PADC_ 16
#define PADC2_ 16
// C = 3 * 2^M3L_ (radix-3 column dimension), R = 2^LGR_ (pure power of two,
// so the corner's bit-reversal algebra over the R index stays exact).
// N1 = R*C = 3*2^26 = 201 326 592 >= n+m+1 = 200 000 001 => M < 0: the wrap
// region is EMPTY and the wrap-correction engine is never used.
#define N3_ 3
#define M3_ (1u << M3L_)
#define NR_ (1u << LGR_)
#define NC_ (N3_ * M3_)
#define N1_ ((u32)(NR_ * NC_))
#define CSX_ (NC_ + PADC_)
#define AOFF(k) (((long)(k) / (long)NC_) * CSX_ + ((long)(k) % (long)NC_))
#define NR2_ (1u << LGR2_)
#define NC2_ (1u << LGC2_)
#define N2_ (1u << (LGR2_ + LGC2_))
#define CSX2_ (NC2_ + PADC2_)
// ============================================================================
// 1002e8 core: NTT mod p = 2013265921 = 15*2^27+1 (primitive root 31)
//
// Transform of length N = R*C ("four step"):
// input (natural flat): x[j1*C + j0], j1 in [0,R), j0 in [0,C)
// scratch (transposed) : S[j0*R + j1]
// spectral output : X[k0*C + k1] holds A_hat[k1*R + k0]
// inverse output : S[j0*R + j1] holds x[j1*C + j0]
// ============================================================================
#include <cstdint>
#define NO_AVX2_UNUSED 0
#include <cstring>
#ifndef NO_AVX2
#include <immintrin.h>
#endif
typedef uint16_t u16;
typedef uint32_t u32;
typedef uint64_t u64;
typedef int64_t i64;
typedef int32_t i32;
typedef unsigned __int128 u128;
#ifndef BSTAIL_SPLIT
#define BSTAIL_SPLIT 188805120L
#endif
static u32 *g_B = 0;
static u32 *g_BT = 0;
static inline u32 *BROW(long el) {
return (el < BSTAIL_SPLIT) ? g_B + el : g_BT + (el - BSTAIL_SPLIT);
}
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __attribute__((always_inline)) __m256i ldu256(const void *p) {
__m256i r; __asm__("vmovdqu %1, %0" : "=x"(r) : "m"(*(const __m256i *)p)); return r;
}
static inline __attribute__((always_inline)) void stu256(void *p, __m256i v) {
__asm__ volatile("vmovdqu %1, %0" : "=m"(*(__m256i *)p) : "x"(v));
}
static inline __attribute__((always_inline)) __m256i vaddm256(__m256i a, __m256i b, __m256i Pv) {
__m256i s = _mm256_add_epi32(a, b);
return _mm256_min_epu32(s, _mm256_sub_epi32(s, Pv));
}
static inline __attribute__((always_inline)) __m256i vsubm256(__m256i a, __m256i b, __m256i Pv) {
__m256i d = _mm256_sub_epi32(a, b);
return _mm256_min_epu32(d, _mm256_add_epi32(d, Pv));
}
#undef CMBSTEP
#pragma GCC pop_options
#pragma GCC optimize("O3","rename-registers","live-range-shrinkage","ira-loop-pressure","no-stack-protector")
#define TS 64
#define PADB_ 16
#define NRMAX_ ((NR_ > NR2_) ? NR_ : NR2_)
static const u32 P = 2'013'265'921u;
static const u32 G = 31u;
static u64 IM;
static u32 TBUF[TS * TS];
static int g_fuse4 = 0;
static u32 TW[TS * TS];
static inline u32 mm(u32 a, u32 b) {
u64 z = (u64)a * b;
u64 q = (u64)(((u128)z * IM) >> 64);
u64 r = z - q * (u64)P;
return (u32)(r >= P ? r - P : r);
}
static inline u32 shoup_c(u32 w) { return (u32)(((u64)w << 32) / (u64)P); }
static u32 pw(u32 a, u32 e) {
u32 r = 1;
while (e) { if (e & 1) r = (u32)((u64)r * a % P); a = (u32)((u64)a * a % P); e >>= 1; }
return r;
}
static void build_tab(u32 *tw, int L, u32 w) {
int off = 0;
for (int len = 2; len <= L; len <<= 1) {
u32 wl = pw(w, (u32)(L / len));
u32 c = 1;
int h = len >> 1;
for (int j = 0; j < h; j++) { tw[off + j] = c; c = mm(c, wl); }
off += h;
}
}
static void build_br(u16 *br, int L) {
int bits = 0; while ((1 << bits) < L) bits++;
for (int i = 0; i < L; i++) {
int r = 0;
for (int b = 0; b < bits; b++) if (i & (1 << b)) r |= 1 << (bits - 1 - b);
br[i] = (u16)r;
}
}
#ifdef NO_AVX2
#define HAVE_AVX2 0
#else
#define HAVE_AVX2 1
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __m256i vmulhi8(__m256i a, __m256i b) {
__m256i e = _mm256_mul_epu32(a, b);
__m256i o = _mm256_mul_epu32(_mm256_bsrli_epi128(a, 4), _mm256_bsrli_epi128(b, 4));
e = _mm256_bsrli_epi128(e, 4);
return _mm256_blend_epi32(e, o, 0xAA);
}
static inline __attribute__((always_inline)) __m256i mullo256(__m256i a, __m256i b) {
__m256i r; __asm__("vpmulld %2, %1, %0" : "=x"(r) : "x"(a), "x"(b)); return r;
}
static inline __m256i shoup8(__m256i v, __m256i w, __m256i ws) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i t = mullo256(v, w);
__m256i q = vmulhi8(v, ws);
__m256i qp = mullo256(q, Pv);
__m256i r = _mm256_sub_epi32(t, qp);
return _mm256_min_epu32(r, _mm256_sub_epi32(r, Pv));
}
static inline __m256i shoup8_broadcast(__m256i v, __m256i w, __m256i ws) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i t = mullo256(v, w);
const __m256i qe = _mm256_mul_epu32(v,ws);
const __m256i qo = _mm256_mul_epu32(_mm256_bsrli_epi128(v,4),ws);
__m256i q = _mm256_blend_epi32(_mm256_bsrli_epi128(qe,4),qo,0xaa);
__m256i qp = mullo256(q, Pv);
__m256i r = _mm256_sub_epi32(t, qp);
return _mm256_min_epu32(r, _mm256_sub_epi32(r, Pv));
}
template<bool Lazy> static inline __m256i shoup_identity_odd8(__m256i v,__m256i w,__m256i ws){
const __m256i pv=_mm256_set1_epi32((int)P);
const __m256i p64=_mm256_set1_epi64x((long long)P);
const __m256i a=_mm256_bsrli_epi128(v,4);
const __m256i prod=_mm256_mul_epu32(a,w);
const __m256i quot=_mm256_bsrli_epi128(_mm256_mul_epu32(a,ws),4);
__m256i r=_mm256_sub_epi64(prod,_mm256_mul_epu32(quot,p64));
__m256i mixed=_mm256_blend_epi32(v,_mm256_bslli_epi128(r,4),0xaa);
return _mm256_min_epu32(mixed,_mm256_sub_epi32(mixed,pv));
}
#define SP(w, ws) _mm256_and_si256(w,M32), _mm256_srli_epi64(w,32), _mm256_and_si256(ws,M32), _mm256_srli_epi64(ws,32)
template<int H> static void dit_pass8(u32 *x,int L,const u32 *tw,const u32 *tws){
const __m256i Pv=_mm256_set1_epi32((int)P);
const u32 *t1=tw+H-1,*s1=tws+H-1;
const u32 *t2=tw+2*H-1,*s2=tws+2*H-1;
const u32 *t3=tw+4*H-1,*s3=tws+4*H-1;
for(int i=0;i<L;i+=8*H){
u32 *p=x+i;
for(int j=0;j<H;j+=8){
__m256i v0=ldu256(p+j),v1=ldu256(p+j+H),
v2=ldu256(p+j+2*H),v3=ldu256(p+j+3*H),
v4=ldu256(p+j+4*H),v5=ldu256(p+j+5*H),
v6=ldu256(p+j+6*H),v7=ldu256(p+j+7*H);
auto dit2=[&](__m256i &a,__m256i &b,__m256i w,__m256i ws,bool lazy=false){
__m256i bm=shoup8(b,w,ws),old=a;
if(lazy){a=_mm256_add_epi32(old,bm);b=_mm256_add_epi32(_mm256_sub_epi32(old,bm),Pv);}
else {a=vaddm256(old,bm,Pv);b=vsubm256(old,bm,Pv);}
};
__m256i w1=ldu256(t1+j),ws1=ldu256(s1+j);
dit2(v0,v1,w1,ws1);dit2(v2,v3,w1,ws1,true);
dit2(v4,v5,w1,ws1);dit2(v6,v7,w1,ws1,true);
__m256i w20=ldu256(t2+j),ws20=ldu256(s2+j);
__m256i w21=ldu256(t2+j+H),ws21=ldu256(s2+j+H);
dit2(v0,v2,w20,ws20);dit2(v4,v6,w20,ws20,true);
dit2(v1,v3,w21,ws21);dit2(v5,v7,w21,ws21,true);
dit2(v0,v4,ldu256(t3+j),ldu256(s3+j));
dit2(v1,v5,ldu256(t3+j+H),ldu256(s3+j+H));
dit2(v2,v6,ldu256(t3+j+2*H),ldu256(s3+j+2*H));
dit2(v3,v7,ldu256(t3+j+3*H),ldu256(s3+j+3*H));
stu256(p+j,v0);stu256(p+j+H,v1);
stu256(p+j+2*H,v2);stu256(p+j+3*H,v3);
stu256(p+j+4*H,v4);stu256(p+j+5*H,v5);
stu256(p+j+6*H,v6);stu256(p+j+7*H,v7);
}
}
}
static inline __m256i shoup_of(__m256i w);
static inline __m256i mont8(__m256i a, __m256i b);
template<bool Product=false,bool LazyTail=false>
static void row_dit_avx(u32 *x, int L, const u32 *tw, const u32 *tws, const u32 *other=nullptr) {
auto initial_load=[&](int offset) {
__m256i v=ldu256(x+offset);
if constexpr(Product){__m256i y=ldu256(other+offset);v=mont8(v,y);}
return v;
};
if constexpr(Product){if(L<8)for(int i=0;i<L;++i)x[i]=mm(x[i],other[i]);}
__m256i Pv = _mm256_set1_epi32((int)P), M1 = _mm256_set1_epi32(-1), Z = _mm256_setzero_si256();
int off = 0;
if (L >= 8) {
const __m256i V4 = _mm256_setr_epi32((int)tw[1], (int)tw[2], 1, 1, (int)tw[1], (int)tw[2], 1, 1);
const __m256i V4s = _mm256_setr_epi32((int)tws[1], (int)tws[2], 2, 2, (int)tws[1], (int)tws[2], 2, 2);
const __m256i V8 = _mm256_setr_epi32((int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6], 1, 1, 1, 1);
const __m256i V8s = _mm256_setr_epi32((int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6], 2, 2, 2, 2);
const __m256i V8Q = _mm256_setr_epi32((int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6],
(int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6]);
const __m256i V8Qs = _mm256_setr_epi32((int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6],
(int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6]);
const __m256i V4P=_mm256_set1_epi64x((long long)tw[2]);
const __m256i V4Ps=_mm256_set1_epi64x((long long)tws[2]);
const __m256i IDX8 = _mm256_setr_epi32(4, 5, 6, 7, 0, 1, 2, 3);
const __m256i IDXQ4 = _mm256_setr_epi32(0, 1, 0, 1, 4, 5, 4, 5);
const __m256i IDXQ8 = _mm256_setr_epi32(0, 1, 2, 3, 0, 1, 2, 3);
int i4 = 0;
for (; i4 + 64 <= L; i4 += 64) {
__m256i vr0 = initial_load(i4 + 0);
__m256i vr1 = initial_load(i4 + 8);
__m256i vr2 = initial_load(i4 + 16);
__m256i vr3 = initial_load(i4 + 24);
__m256i vr4 = initial_load(i4 + 32);
__m256i vr5 = initial_load(i4 + 40);
__m256i vr6 = initial_load(i4 + 48);
__m256i vr7 = initial_load(i4 + 56);
{
{
__m256i y = _mm256_shuffle_epi32(vr0, 0xB1);
__m256i sa = _mm256_add_epi32(vr0,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr0),Pv);
vr0 = _mm256_blend_epi32(sa,d,0xAA);
vr0 = _mm256_min_epu32(vr0,_mm256_sub_epi32(vr0,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr1, 0xB1);
__m256i sa = _mm256_add_epi32(vr1,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr1),Pv);
vr1 = _mm256_blend_epi32(sa,d,0xAA);
vr1 = _mm256_min_epu32(vr1,_mm256_sub_epi32(vr1,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr2, 0xB1);
__m256i sa = _mm256_add_epi32(vr2,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr2),Pv);
vr2 = _mm256_blend_epi32(sa,d,0xAA);
vr2 = _mm256_min_epu32(vr2,_mm256_sub_epi32(vr2,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr3, 0xB1);
__m256i sa = _mm256_add_epi32(vr3,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr3),Pv);
vr3 = _mm256_blend_epi32(sa,d,0xAA);
vr3 = _mm256_min_epu32(vr3,_mm256_sub_epi32(vr3,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr4, 0xB1);
__m256i sa = _mm256_add_epi32(vr4,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr4),Pv);
vr4 = _mm256_blend_epi32(sa,d,0xAA);
vr4 = _mm256_min_epu32(vr4,_mm256_sub_epi32(vr4,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr5, 0xB1);
__m256i sa = _mm256_add_epi32(vr5,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr5),Pv);
vr5 = _mm256_blend_epi32(sa,d,0xAA);
vr5 = _mm256_min_epu32(vr5,_mm256_sub_epi32(vr5,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr6, 0xB1);
__m256i sa = _mm256_add_epi32(vr6,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr6),Pv);
vr6 = _mm256_blend_epi32(sa,d,0xAA);
vr6 = _mm256_min_epu32(vr6,_mm256_sub_epi32(vr6,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr7, 0xB1);
__m256i sa = _mm256_add_epi32(vr7,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,vr7),Pv);
vr7 = _mm256_blend_epi32(sa,d,0xAA);
vr7 = _mm256_min_epu32(vr7,_mm256_sub_epi32(vr7,Pv));
}
}
{
{
__m256i PP = _mm256_unpackhi_epi64(vr0,vr1);
__m256i MM = shoup_identity_odd8<false>(PP, V4P, V4Ps);
__m256i LL = _mm256_unpacklo_epi64(vr0,vr1);
__m256i SS = _mm256_min_epu32(_mm256_add_epi32(LL, MM), _mm256_sub_epi32(_mm256_add_epi32(LL, MM), Pv));
__m256i DD = _mm256_min_epu32(_mm256_sub_epi32(LL, MM), _mm256_add_epi32(_mm256_sub_epi32(LL, MM), Pv));
vr0 = _mm256_unpacklo_epi64(SS,DD);
vr1 = _mm256_unpackhi_epi64(SS,DD);
}
}
{
{
__m256i PP = _mm256_unpackhi_epi64(vr2,vr3);
__m256i MM = shoup_identity_odd8<false>(PP, V4P, V4Ps);
__m256i LL = _mm256_unpacklo_epi64(vr2,vr3);
__m256i SS = _mm256_min_epu32(_mm256_add_epi32(LL, MM), _mm256_sub_epi32(_mm256_add_epi32(LL, MM), Pv));
__m256i DD = _mm256_min_epu32(_mm256_sub_epi32(LL, MM), _mm256_add_epi32(_mm256_sub_epi32(LL, MM), Pv));
vr2 = _mm256_unpacklo_epi64(SS,DD);
vr3 = _mm256_unpackhi_epi64(SS,DD);
}
}
{
{
__m256i PP = _mm256_unpackhi_epi64(vr4,vr5);
__m256i MM = shoup_identity_odd8<false>(PP, V4P, V4Ps);
__m256i LL = _mm256_unpacklo_epi64(vr4,vr5);
__m256i SS = _mm256_min_epu32(_mm256_add_epi32(LL, MM), _mm256_sub_epi32(_mm256_add_epi32(LL, MM), Pv));
__m256i DD = _mm256_min_epu32(_mm256_sub_epi32(LL, MM), _mm256_add_epi32(_mm256_sub_epi32(LL, MM), Pv));
vr4 = _mm256_unpacklo_epi64(SS,DD);
vr5 = _mm256_unpackhi_epi64(SS,DD);
}
}
{
{
__m256i PP = _mm256_unpackhi_epi64(vr6,vr7);
__m256i MM = shoup_identity_odd8<false>(PP, V4P, V4Ps);
__m256i LL = _mm256_unpacklo_epi64(vr6,vr7);
__m256i SS = _mm256_min_epu32(_mm256_add_epi32(LL, MM), _mm256_sub_epi32(_mm256_add_epi32(LL, MM), Pv));
__m256i DD = _mm256_min_epu32(_mm256_sub_epi32(LL, MM), _mm256_add_epi32(_mm256_sub_epi32(LL, MM), Pv));
vr6 = _mm256_unpacklo_epi64(SS,DD);
vr7 = _mm256_unpackhi_epi64(SS,DD);
}
}
{
{
__m256i YX = _mm256_permute2x128_si256(vr0, vr1, 0x31);
__m256i MM = shoup8(YX, V8Q, V8Qs);
__m256i XY = _mm256_permute2x128_si256(vr0, vr1, 0x20);
__m256i SS = _mm256_add_epi32(XY,MM);
__m256i DD = _mm256_add_epi32(_mm256_sub_epi32(XY,MM),Pv);
vr0 = _mm256_permute2x128_si256(SS, DD, 0x20);
vr0 = _mm256_min_epu32(vr0,_mm256_sub_epi32(vr0,Pv));
vr1 = _mm256_permute2x128_si256(SS, DD, 0x31);
}
}
{
{
__m256i YX = _mm256_permute2x128_si256(vr2, vr3, 0x31);
__m256i MM = shoup8(YX, V8Q, V8Qs);
__m256i XY = _mm256_permute2x128_si256(vr2, vr3, 0x20);
__m256i SS = _mm256_add_epi32(XY,MM);
__m256i DD = _mm256_add_epi32(_mm256_sub_epi32(XY,MM),Pv);
vr2 = _mm256_permute2x128_si256(SS, DD, 0x20);
vr2 = _mm256_min_epu32(vr2,_mm256_sub_epi32(vr2,Pv));
vr3 = _mm256_permute2x128_si256(SS, DD, 0x31);
}
}
{
{
__m256i YX = _mm256_permute2x128_si256(vr4, vr5, 0x31);
__m256i MM = shoup8(YX, V8Q, V8Qs);
__m256i XY = _mm256_permute2x128_si256(vr4, vr5, 0x20);
__m256i SS = _mm256_add_epi32(XY,MM);
__m256i DD = _mm256_add_epi32(_mm256_sub_epi32(XY,MM),Pv);
vr4 = _mm256_permute2x128_si256(SS, DD, 0x20);
vr4 = _mm256_min_epu32(vr4,_mm256_sub_epi32(vr4,Pv));
vr5 = _mm256_permute2x128_si256(SS, DD, 0x31);
}
}
{
{
__m256i YX = _mm256_permute2x128_si256(vr6, vr7, 0x31);
__m256i MM = shoup8(YX, V8Q, V8Qs);
__m256i XY = _mm256_permute2x128_si256(vr6, vr7, 0x20);
__m256i SS = _mm256_add_epi32(XY,MM);
__m256i DD = _mm256_add_epi32(_mm256_sub_epi32(XY,MM),Pv);
vr6 = _mm256_permute2x128_si256(SS, DD, 0x20);
vr6 = _mm256_min_epu32(vr6,_mm256_sub_epi32(vr6,Pv));
vr7 = _mm256_permute2x128_si256(SS, DD, 0x31);
}
}
stu256((void *)((x + i4 + 0)), vr0);
stu256((void *)((x + i4 + 8)), vr1);
stu256((void *)((x + i4 + 16)), vr2);
stu256((void *)((x + i4 + 24)), vr3);
stu256((void *)((x + i4 + 32)), vr4);
stu256((void *)((x + i4 + 40)), vr5);
stu256((void *)((x + i4 + 48)), vr6);
stu256((void *)((x + i4 + 56)), vr7);
}
for (; i4 < L; i4 += 8) {
__m256i v = initial_load(i4);
{
__m256i y = _mm256_shuffle_epi32(v, 0xB1);
__m256i sa = _mm256_add_epi32(v,y);
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(y,v),Pv);
v = _mm256_blend_epi32(sa,d,0xAA);
v = _mm256_min_epu32(v,_mm256_sub_epi32(v,Pv));
}
{
__m256i y = _mm256_shuffle_epi32(v, 0x4E);
y = shoup8(y, V4, V4s);
__m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
__m256i d = _mm256_min_epu32(_mm256_sub_epi32(v, y), _mm256_add_epi32(_mm256_sub_epi32(v, y), Pv));
d = _mm256_permutevar8x32_epi32(d, IDXQ4);
v = _mm256_blend_epi32(sa, d, 0xCC);
}
{
__m256i y = _mm256_permutevar8x32_epi32(v, IDX8);
y = shoup8(y, V8, V8s);
__m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
__m256i d = _mm256_min_epu32(_mm256_sub_epi32(v, y), _mm256_add_epi32(_mm256_sub_epi32(v, y), Pv));
d = _mm256_permutevar8x32_epi32(d, IDXQ8);
v = _mm256_blend_epi32(sa, d, 0xF0);
}
stu256((void *)((x + i4)), v);
}
off = 7;
int len = 16;
if(L==4096 || L==8192 || L==16384){
dit_pass8<8>(x,L,tw,tws);
dit_pass8<64>(x,L,tw,tws);
dit_pass8<512>(x,L,tw,tws);
len=8192;off=4095;
}
for (; (len << 1) <= L; len <<= 2) {
int h = len >> 1;
const u32 *t1 = tw + off, *ts1 = tws + off;
const u32 *t2 = t1 + h, *ts2 = ts1 + h;
const u32 *t3 = t2 + h, *ts3 = ts2 + h;
off += 3 * h;
for (int i = 0; i < L; i += (len << 1)) {
u32 *p = x + i;
for (int j = 0; j < h; j += 8) {
__m256i WV1 = ldu256((const void *)((t1 + j))), WS1 = ldu256((const void *)((ts1 + j)));
__m256i WV2 = ldu256((const void *)((t2 + j))), WS2 = ldu256((const void *)((ts2 + j)));
__m256i WV3 = ldu256((const void *)((t3 + j))), WS3 = ldu256((const void *)((ts3 + j)));
__m256i A = ldu256((const void *)(p + j));
__m256i B = ldu256((const void *)(p + j + h));
__m256i C = ldu256((const void *)(p + j + 2 * h));
__m256i D = ldu256((const void *)(p + j + 3 * h));
__m256i Bm = shoup8(B, WV1, WS1), Dm = shoup8(D, WV1, WS1);
__m256i Ap = vaddm256(A, Bm, Pv), Bp = vsubm256(A, Bm, Pv);
__m256i Cp = _mm256_add_epi32(C, Dm), Dp = _mm256_add_epi32(_mm256_sub_epi32(C, Dm), Pv);
__m256i Cm = shoup8(Cp, WV2, WS2), Dn = shoup8(Dp, WV3, WS3);
if(LazyTail && (len<<1)==L){
stu256(p+j,_mm256_add_epi32(Ap,Cm));
stu256(p+j+2*h,_mm256_add_epi32(_mm256_sub_epi32(Ap,Cm),Pv));
stu256(p+j+h,_mm256_add_epi32(Bp,Dn));
stu256(p+j+3*h,_mm256_add_epi32(_mm256_sub_epi32(Bp,Dn),Pv));
}else{
stu256((void *)(p + j), vaddm256(Ap, Cm, Pv));
stu256((void *)(p + j + 2 * h), vsubm256(Ap, Cm, Pv));
stu256((void *)(p + j + h), vaddm256(Bp, Dn, Pv));
stu256((void *)(p + j + 3 * h), vsubm256(Bp, Dn, Pv));
}
}
}
}
if (len <= L) {
int h = len >> 1;
const u32 *t = tw + off, *ts = tws + off; off += h;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j += 8) {
__m256i u = ldu256((const void *)(p + j));
__m256i v = shoup8(ldu256((const void *)(q + j)),
ldu256((const void *)(t + j)), ldu256((const void *)(ts + j)));
stu256((void *)(p + j), vaddm256(u, v, Pv));
stu256((void *)(q + j), vsubm256(u, v, Pv));
}
}
}
} else {
for (int len = 2; len <= L; len <<= 1) {
int h = len >> 1;
const u32 *t = tw + off; off += h;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = mm(q[j], t[j]);
u32 sa = u + v; if (sa >= P) sa -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = sa; q[j] = d;
}
}
}
}
}
template<int H3, bool HalfZero=false> static void dif_pass8(u32 *x, int L, const u32 *tw, const u32 *tws) {
const __m256i Pv = _mm256_set1_epi32((int)P);
const int len = H3 << 3;
const u32 *t1 = tw + (len >> 1) - 1, *s1 = tws + (len >> 1) - 1;
const u32 *t2 = tw + (len >> 2) - 1, *s2 = tws + (len >> 2) - 1;
const u32 *t3 = tw + H3 - 1, *s3 = tws + H3 - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i;
for (int j = 0; j < H3; j += 8) {
__m256i v0=ldu256(p+j),v1=ldu256(p+j+H3),
v2=ldu256(p+j+2*H3),v3=ldu256(p+j+3*H3);
__m256i v4,v5,v6,v7;
auto dif2=[&](__m256i &a,__m256i &b,__m256i w,__m256i ws) {
__m256i sm=vaddm256(a,b,Pv),df=_mm256_add_epi32(_mm256_sub_epi32(a,b),Pv);
a=sm;b=shoup8(df,w,ws);
};
if constexpr(HalfZero) {
v4=shoup8(v0,ldu256(t1+j),ldu256(s1+j));
v5=shoup8(v1,ldu256(t1+j+H3),ldu256(s1+j+H3));
v6=shoup8(v2,ldu256(t1+j+2*H3),ldu256(s1+j+2*H3));
v7=shoup8(v3,ldu256(t1+j+3*H3),ldu256(s1+j+3*H3));
} else {
v4=ldu256(p+j+4*H3);v5=ldu256(p+j+5*H3);
v6=ldu256(p+j+6*H3);v7=ldu256(p+j+7*H3);
dif2(v0,v4,ldu256(t1+j),ldu256(s1+j));
dif2(v1,v5,ldu256(t1+j+H3),ldu256(s1+j+H3));
dif2(v2,v6,ldu256(t1+j+2*H3),ldu256(s1+j+2*H3));
dif2(v3,v7,ldu256(t1+j+3*H3),ldu256(s1+j+3*H3));
}
__m256i w20=ldu256(t2+j),s20=ldu256(s2+j);
__m256i w21=ldu256(t2+j+H3),s21=ldu256(s2+j+H3);
dif2(v0,v2,w20,s20);dif2(v4,v6,w20,s20);
dif2(v1,v3,w21,s21);dif2(v5,v7,w21,s21);
__m256i w3=ldu256(t3+j),ws3=ldu256(s3+j);
dif2(v0,v1,w3,ws3);dif2(v2,v3,w3,ws3);
dif2(v4,v5,w3,ws3);dif2(v6,v7,w3,ws3);
stu256(p+j,v0);stu256(p+j+H3,v1);
stu256(p+j+2*H3,v2);stu256(p+j+3*H3,v3);
stu256(p+j+4*H3,v4);stu256(p+j+5*H3,v5);
stu256(p+j+6*H3,v6);stu256(p+j+7*H3,v7);
}
}
}
template<bool HalfZero=false, bool LazyOut=false>
static void row_dif_avx(u32 *x, int L, const u32 *tw, const u32 *tws) {
__m256i Pv = _mm256_set1_epi32((int)P), M1 = _mm256_set1_epi32(-1), Z = _mm256_setzero_si256();
if (L >= 8) {
int len = L;
if (L == 4096) {
dif_pass8<512,HalfZero>(x, L, tw, tws); dif_pass8<64>(x, L, tw, tws); dif_pass8<8>(x, L, tw, tws);
len = 8;
} else if (L == 8192) {
dif_pass8<1024>(x, L, tw, tws); dif_pass8<128>(x, L, tw, tws); dif_pass8<16>(x, L, tw, tws);
len = 16;
}
if (L == 16384) {
dif_pass8<2048>(x,L,tw,tws); dif_pass8<256>(x,L,tw,tws); dif_pass8<32>(x,L,tw,tws);
len=32;
}
if (L != 4096 && L != 8192)
for (; len >= 32; len >>= 2) {
int h = len >> 1, h2 = len >> 2;
const u32 *t1 = tw + h - 1, *ts1 = tws + h - 1;
const u32 *t2 = t1 + h2, *ts2 = ts1 + h2;
const u32 *t3 = tw + h2 - 1, *ts3 = tws + h2 - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i;
for (int j = 0; j < h2; j += 8) {
__m256i WV1 = ldu256((const void *)((t1 + j))), WS1 = ldu256((const void *)((ts1 + j)));
__m256i WV2 = ldu256((const void *)((t2 + j))), WS2 = ldu256((const void *)((ts2 + j)));
__m256i WV3 = ldu256((const void *)((t3 + j))), WS3 = ldu256((const void *)((ts3 + j)));
__m256i v0 = ldu256((const void *)((p + j)));
__m256i v1 = ldu256((const void *)((p + j + h2)));
__m256i v2 = ldu256((const void *)((p + j + h)));
__m256i v3 = ldu256((const void *)((p + j + h + h2)));
__m256i s02 = _mm256_min_epu32(_mm256_add_epi32(v0, v2), _mm256_sub_epi32(_mm256_add_epi32(v0, v2), Pv));
__m256i d02 = _mm256_sub_epi32(v0, v2);
__m256i s13 = _mm256_min_epu32(_mm256_add_epi32(v1, v3), _mm256_sub_epi32(_mm256_add_epi32(v1, v3), Pv));
__m256i d13 = _mm256_sub_epi32(v1, v3);
d02 = _mm256_add_epi32(d02, Pv);
d13 = _mm256_add_epi32(d13, Pv);
d02 = shoup8(d02, WV1, WS1);
d13 = shoup8(d13, WV2, WS2);
__m256i ra = _mm256_min_epu32(_mm256_add_epi32(s02, s13), _mm256_sub_epi32(_mm256_add_epi32(s02, s13), Pv));
__m256i rb = _mm256_add_epi32(_mm256_sub_epi32(s02, s13), Pv);
rb = shoup8(rb, WV3, WS3);
__m256i rc = _mm256_min_epu32(_mm256_add_epi32(d02, d13), _mm256_sub_epi32(_mm256_add_epi32(d02, d13), Pv));
__m256i rd = _mm256_add_epi32(_mm256_sub_epi32(d02, d13), Pv);
rd = shoup8(rd, WV3, WS3);
stu256((void *)((p + j)), ra);
stu256((void *)((p + j + h2)), rb);
stu256((void *)((p + j + h)), rc);
stu256((void *)((p + j + h + h2)), rd);
}
}
}
if (len >= 16) {
int h = len >> 1;
const u32 *t = tw + h - 1, *ts = tws + h - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j += 8) {
__m256i u = ldu256((const void *)((p + j)));
__m256i v = ldu256((const void *)((q + j)));
__m256i w = ldu256((const void *)((t + j)));
__m256i ws = ldu256((const void *)((ts + j)));
__m256i sa = _mm256_min_epu32(_mm256_add_epi32(u, v), _mm256_sub_epi32(_mm256_add_epi32(u, v), Pv));
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(u, v), Pv);
d = shoup8(d, w, ws);
stu256((void *)((p + j)), sa);
stu256((void *)((q + j)), d);
}
}
}
const __m256i V8 = _mm256_setr_epi32(1, 1, 1, 1, (int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6]);
const __m256i V8s = _mm256_setr_epi32(2, 2, 2, 2, (int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6]);
const __m256i V8N = _mm256_setr_epi32((int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6],
(int)(P - tw[3]), (int)(P - tw[4]), (int)(P - tw[5]), (int)(P - tw[6]));
const __m256i V8Ns = _mm256_setr_epi32((int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6],
(int)shoup_c(P - tw[3]), (int)shoup_c(P - tw[4]), (int)shoup_c(P - tw[5]), (int)shoup_c(P - tw[6]));
const __m256i V4 = _mm256_setr_epi32(1, 1, (int)tw[1], (int)tw[2], 1, 1, (int)tw[1], (int)tw[2]);
const __m256i V4s = _mm256_setr_epi32(2, 2, (int)tws[1], (int)tws[2], 2, 2, (int)tws[1], (int)tws[2]);
const __m256i V4P=_mm256_set1_epi64x((long long)tw[2]);
const __m256i V4Ps=_mm256_set1_epi64x((long long)tws[2]);
const __m256i IDX8 = _mm256_setr_epi32(4, 5, 6, 7, 0, 1, 2, 3);
int i4 = 0;
for (; i4 + 64 <= L; i4 += 64) {
__m256i vr0 = ldu256((const void *)((x + i4 + 0)));
__m256i vr1 = ldu256((const void *)((x + i4 + 8)));
__m256i vr2 = ldu256((const void *)((x + i4 + 16)));
__m256i vr3 = ldu256((const void *)((x + i4 + 24)));
__m256i vr4 = ldu256((const void *)((x + i4 + 32)));
__m256i vr5 = ldu256((const void *)((x + i4 + 40)));
__m256i vr6 = ldu256((const void *)((x + i4 + 48)));
__m256i vr7 = ldu256((const void *)((x + i4 + 56)));
{
{
__m256i X = _mm256_permute2x128_si256(vr0, vr1, 0x21);
__m256i Y = _mm256_permute2x128_si256(vr0, vr1, 0x30);
__m256i S = _mm256_min_epu32(_mm256_add_epi32(X, Y), _mm256_sub_epi32(_mm256_add_epi32(X, Y), Pv));
__m256i D = _mm256_add_epi32(_mm256_sub_epi32(Y, X), Pv);
D = shoup8(D, V8N, V8Ns);
vr0 = _mm256_permute2x128_si256(S, D, 0x20);
vr1 = _mm256_permute2x128_si256(S, D, 0x31);
}
}
{
{
__m256i X = _mm256_permute2x128_si256(vr2, vr3, 0x21);
__m256i Y = _mm256_permute2x128_si256(vr2, vr3, 0x30);
__m256i S = _mm256_min_epu32(_mm256_add_epi32(X, Y), _mm256_sub_epi32(_mm256_add_epi32(X, Y), Pv));
__m256i D = _mm256_add_epi32(_mm256_sub_epi32(Y, X), Pv);
D = shoup8(D, V8N, V8Ns);
vr2 = _mm256_permute2x128_si256(S, D, 0x20);
vr3 = _mm256_permute2x128_si256(S, D, 0x31);
}
}
{
{
__m256i X = _mm256_permute2x128_si256(vr4, vr5, 0x21);
__m256i Y = _mm256_permute2x128_si256(vr4, vr5, 0x30);
__m256i S = _mm256_min_epu32(_mm256_add_epi32(X, Y), _mm256_sub_epi32(_mm256_add_epi32(X, Y), Pv));
__m256i D = _mm256_add_epi32(_mm256_sub_epi32(Y, X), Pv);
D = shoup8(D, V8N, V8Ns);
vr4 = _mm256_permute2x128_si256(S, D, 0x20);
vr5 = _mm256_permute2x128_si256(S, D, 0x31);
}
}
{
{
__m256i X = _mm256_permute2x128_si256(vr6, vr7, 0x21);
__m256i Y = _mm256_permute2x128_si256(vr6, vr7, 0x30);
__m256i S = _mm256_min_epu32(_mm256_add_epi32(X, Y), _mm256_sub_epi32(_mm256_add_epi32(X, Y), Pv));
__m256i D = _mm256_add_epi32(_mm256_sub_epi32(Y, X), Pv);
D = shoup8(D, V8N, V8Ns);
vr6 = _mm256_permute2x128_si256(S, D, 0x20);
vr7 = _mm256_permute2x128_si256(S, D, 0x31);
}
}
{
{
__m256i yu = _mm256_shuffle_epi32(vr0, 0x4E);
__m256i saU = _mm256_min_epu32(_mm256_add_epi32(vr0, yu), _mm256_sub_epi32(_mm256_add_epi32(vr0, yu), Pv));
__m256i du = _mm256_sub_epi32(yu, vr0);
du = _mm256_add_epi32(du, Pv);
__m256i yv = _mm256_shuffle_epi32(vr1, 0x4E);
__m256i saV = _mm256_min_epu32(_mm256_add_epi32(vr1, yv), _mm256_sub_epi32(_mm256_add_epi32(vr1, yv), Pv));
__m256i dv = _mm256_sub_epi32(yv, vr1);
dv = _mm256_add_epi32(dv, Pv);
__m256i PP = _mm256_unpackhi_epi64(du,dv);
__m256i MM = shoup_identity_odd8<true>(PP, V4P, V4Ps);
vr0 = _mm256_blend_epi32(saU, _mm256_slli_si256(MM,8), 0xCC);
vr1 = _mm256_blend_epi32(saV, MM, 0xCC);
}
}
{
{
__m256i yu = _mm256_shuffle_epi32(vr2, 0x4E);
__m256i saU = _mm256_min_epu32(_mm256_add_epi32(vr2, yu), _mm256_sub_epi32(_mm256_add_epi32(vr2, yu), Pv));
__m256i du = _mm256_sub_epi32(yu, vr2);
du = _mm256_add_epi32(du, Pv);
__m256i yv = _mm256_shuffle_epi32(vr3, 0x4E);
__m256i saV = _mm256_min_epu32(_mm256_add_epi32(vr3, yv), _mm256_sub_epi32(_mm256_add_epi32(vr3, yv), Pv));
__m256i dv = _mm256_sub_epi32(yv, vr3);
dv = _mm256_add_epi32(dv, Pv);
__m256i PP = _mm256_unpackhi_epi64(du,dv);
__m256i MM = shoup_identity_odd8<true>(PP, V4P, V4Ps);
vr2 = _mm256_blend_epi32(saU, _mm256_slli_si256(MM,8), 0xCC);
vr3 = _mm256_blend_epi32(saV, MM, 0xCC);
}
}
{
{
__m256i yu = _mm256_shuffle_epi32(vr4, 0x4E);
__m256i saU = _mm256_min_epu32(_mm256_add_epi32(vr4, yu), _mm256_sub_epi32(_mm256_add_epi32(vr4, yu), Pv));
__m256i du = _mm256_sub_epi32(yu, vr4);
du = _mm256_add_epi32(du, Pv);
__m256i yv = _mm256_shuffle_epi32(vr5, 0x4E);
__m256i saV = _mm256_min_epu32(_mm256_add_epi32(vr5, yv), _mm256_sub_epi32(_mm256_add_epi32(vr5, yv), Pv));
__m256i dv = _mm256_sub_epi32(yv, vr5);
dv = _mm256_add_epi32(dv, Pv);
__m256i PP = _mm256_unpackhi_epi64(du,dv);
__m256i MM = shoup_identity_odd8<true>(PP, V4P, V4Ps);
vr4 = _mm256_blend_epi32(saU, _mm256_slli_si256(MM,8), 0xCC);
vr5 = _mm256_blend_epi32(saV, MM, 0xCC);
}
}
{
{
__m256i yu = _mm256_shuffle_epi32(vr6, 0x4E);
__m256i saU = _mm256_min_epu32(_mm256_add_epi32(vr6, yu), _mm256_sub_epi32(_mm256_add_epi32(vr6, yu), Pv));
__m256i du = _mm256_sub_epi32(yu, vr6);
du = _mm256_add_epi32(du, Pv);
__m256i yv = _mm256_shuffle_epi32(vr7, 0x4E);
__m256i saV = _mm256_min_epu32(_mm256_add_epi32(vr7, yv), _mm256_sub_epi32(_mm256_add_epi32(vr7, yv), Pv));
__m256i dv = _mm256_sub_epi32(yv, vr7);
dv = _mm256_add_epi32(dv, Pv);
__m256i PP = _mm256_unpackhi_epi64(du,dv);
__m256i MM = shoup_identity_odd8<true>(PP, V4P, V4Ps);
vr6 = _mm256_blend_epi32(saU, _mm256_slli_si256(MM,8), 0xCC);
vr7 = _mm256_blend_epi32(saV, MM, 0xCC);
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr0, 0xB1);
__m256i sa=_mm256_add_epi32(vr0,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr0),Pv);
vr0=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr0=_mm256_min_epu32(vr0,_mm256_sub_epi32(vr0,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr1, 0xB1);
__m256i sa=_mm256_add_epi32(vr1,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr1),Pv);
vr1=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr1=_mm256_min_epu32(vr1,_mm256_sub_epi32(vr1,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr2, 0xB1);
__m256i sa=_mm256_add_epi32(vr2,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr2),Pv);
vr2=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr2=_mm256_min_epu32(vr2,_mm256_sub_epi32(vr2,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr3, 0xB1);
__m256i sa=_mm256_add_epi32(vr3,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr3),Pv);
vr3=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr3=_mm256_min_epu32(vr3,_mm256_sub_epi32(vr3,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr4, 0xB1);
__m256i sa=_mm256_add_epi32(vr4,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr4),Pv);
vr4=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr4=_mm256_min_epu32(vr4,_mm256_sub_epi32(vr4,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr5, 0xB1);
__m256i sa=_mm256_add_epi32(vr5,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr5),Pv);
vr5=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr5=_mm256_min_epu32(vr5,_mm256_sub_epi32(vr5,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr6, 0xB1);
__m256i sa=_mm256_add_epi32(vr6,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr6),Pv);
vr6=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr6=_mm256_min_epu32(vr6,_mm256_sub_epi32(vr6,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr7, 0xB1);
__m256i sa=_mm256_add_epi32(vr7,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,vr7),Pv);
vr7=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)vr7=_mm256_min_epu32(vr7,_mm256_sub_epi32(vr7,Pv));
}
}
stu256((void *)((x + i4 + 0)), vr0);
stu256((void *)((x + i4 + 8)), vr1);
stu256((void *)((x + i4 + 16)), vr2);
stu256((void *)((x + i4 + 24)), vr3);
stu256((void *)((x + i4 + 32)), vr4);
stu256((void *)((x + i4 + 40)), vr5);
stu256((void *)((x + i4 + 48)), vr6);
stu256((void *)((x + i4 + 56)), vr7);
}
for (; i4 < L; i4 += 8) {
__m256i v = ldu256((const void *)((x + i4)));
{
__m256i y = _mm256_permutevar8x32_epi32(v, IDX8);
__m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
__m256i d0 = _mm256_sub_epi32(y, v);
d0 = _mm256_add_epi32(d0, Pv);
__m256i d = shoup8(d0, V8, V8s);
v = _mm256_blend_epi32(sa, d, 0xF0);
}
{
__m256i y = _mm256_shuffle_epi32(v, 0x4E);
__m256i sa = _mm256_min_epu32(_mm256_add_epi32(v, y), _mm256_sub_epi32(_mm256_add_epi32(v, y), Pv));
__m256i d = _mm256_sub_epi32(y, v);
d = _mm256_add_epi32(d, Pv);
d = shoup8(d, V4, V4s);
v = _mm256_blend_epi32(sa, d, 0xCC);
}
{
__m256i y = _mm256_shuffle_epi32(v, 0xB1);
__m256i sa=_mm256_add_epi32(v,y);
__m256i d=_mm256_add_epi32(_mm256_sub_epi32(y,v),Pv);
v=_mm256_blend_epi32(sa,d,0xAA);
if constexpr(!LazyOut)v=_mm256_min_epu32(v,_mm256_sub_epi32(v,Pv));
}
stu256((void *)((x + i4)), v);
}
} else {
for (int len = L; len >= 2; len >>= 1) {
int h = len >> 1;
const u32 *t = tw + h - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = q[j];
u32 sa = u + v; if (sa >= P) sa -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = sa; q[j] = mm(d, t[j]);
}
}
}
}
}
static inline __m256i shoup_of(__m256i w) {
const __m256i ML = _mm256_set1_epi32(572662301);
__m256i lo = _mm256_mul_epu32(w, ML);
__m256i hi = _mm256_mul_epu32(_mm256_srli_epi64(w, 32), ML);
__m256i m = _mm256_blend_epi32(_mm256_srli_epi64(lo, 32), _mm256_slli_epi64(_mm256_srli_epi64(hi, 32), 32), 0xAA);
return _mm256_add_epi32(_mm256_add_epi32(w, w), m);
}
static inline void t8x8(__m256i *r) {
__m256i t0 = _mm256_unpacklo_epi32(r[0], r[1]), t1 = _mm256_unpackhi_epi32(r[0], r[1]);
__m256i t2 = _mm256_unpacklo_epi32(r[2], r[3]), t3 = _mm256_unpackhi_epi32(r[2], r[3]);
__m256i t4 = _mm256_unpacklo_epi32(r[4], r[5]), t5 = _mm256_unpackhi_epi32(r[4], r[5]);
__m256i t6 = _mm256_unpacklo_epi32(r[6], r[7]), t7 = _mm256_unpackhi_epi32(r[6], r[7]);
__m256i u0 = _mm256_unpacklo_epi64(t0, t2), u1 = _mm256_unpackhi_epi64(t0, t2);
__m256i u2 = _mm256_unpacklo_epi64(t1, t3), u3 = _mm256_unpackhi_epi64(t1, t3);
__m256i u4 = _mm256_unpacklo_epi64(t4, t6), u5 = _mm256_unpackhi_epi64(t4, t6);
__m256i u6 = _mm256_unpacklo_epi64(t5, t7), u7 = _mm256_unpackhi_epi64(t5, t7);
r[0] = _mm256_permute2x128_si256(u0, u4, 0x20); r[1] = _mm256_permute2x128_si256(u1, u5, 0x20);
r[2] = _mm256_permute2x128_si256(u2, u6, 0x20); r[3] = _mm256_permute2x128_si256(u3, u7, 0x20);
r[4] = _mm256_permute2x128_si256(u0, u4, 0x31); r[5] = _mm256_permute2x128_si256(u1, u5, 0x31);
r[6] = _mm256_permute2x128_si256(u2, u6, 0x31); r[7] = _mm256_permute2x128_si256(u3, u7, 0x31);
}
static void pointwise_mul_avx(u32 *A, const u32 *B, u32 n, u32 sc, u32 sx) {
(void)sc;
long k = 0;
for (int r = 0; k < (long)n; r++) {
long lim = (long)sx - PADC_;
if (lim > (long)n - k) lim = (long)n - k;
u32 *a = A + (long)r * sx;
const u32 *b = BROW((long)r * sx);
long i = 0;
for (; i + 8 <= lim; i += 8) {
__m256i x = ldu256((const void *)(a + i));
__m256i y = ldu256((const void *)(b + i));
stu256((void *)(a + i), shoup8(x, y, shoup_of(y)));
}
for (; i < lim; i++) a[i] = mm(a[i], b[i]);
k += lim;
}
}
static u32 NPI = 0;
static const u32 R32 = (u32)((1ull << 32) % P);
static inline u32 mont1(u32 a, u32 b) {
u64 z = (u64)a * (u64)b;
u32 m = (u32)((u64)(u32)z * (u64)NPI);
u64 t = (z + (u64)m * (u64)P) >> 32;
return (u32)(t >= P ? t - P : t);
}
static inline __m256i mont8(__m256i a, __m256i b) {
__m256i P64 = _mm256_set1_epi64x((long long)P), NP = _mm256_set1_epi64x((long long)NPI);
__m256i Pv = _mm256_set1_epi32((int)P);
__m256i ahi = _mm256_bsrli_epi128(a, 4), bhi = _mm256_bsrli_epi128(b, 4);
__m256i tlo = _mm256_mul_epu32(a, b), thi = _mm256_mul_epu32(ahi, bhi);
__m256i e2 = _mm256_mul_epu32(tlo, NP), o2 = _mm256_mul_epu32(thi, NP);
__m256i e4 = _mm256_mul_epu32(e2, P64), o4 = _mm256_mul_epu32(o2, P64);
__m256i e5 = _mm256_add_epi64(tlo, e4), o5 = _mm256_add_epi64(thi, o4);
__m256i u = _mm256_or_si256(_mm256_bsrli_epi128(e5, 4), o5);
return _mm256_min_epu32(u, _mm256_sub_epi32(u, Pv));
}
static void scale_row_prepared(u32 *dst, const u32 *src, u32 k, u32 ks) {
const __m256i kv = _mm256_set1_epi32((int)k);
const __m256i ksv = _mm256_set1_epi32((int)ks);
for (int v = 0; v < TS; v += 8)
stu256(dst + v, shoup8(ldu256(src + v), kv, ksv));
}
static void scale_row_mont(u32 *dst, const u32 *src, u32 kM) {
const u32 k = mont1(kM, 1u);
const __m256i kv = _mm256_set1_epi32((int)k);
const __m256i ks = _mm256_set1_epi32((int)shoup_c(k));
for (int v = 0; v < TS; v += 8)
stu256(dst + v, shoup8(ldu256(src + v), kv, ks));
}
static void twiddle_rows_mont(u32 *TW, const u32 *RVM) {
__m256i factor[TS / 8], factor_shoup[TS / 8];
const __m256i one = _mm256_set1_epi32(1);
for (int v = 0; v < TS; v += 8) {
factor[v / 8] = mont8(ldu256(RVM + v), one);
factor_shoup[v / 8] = shoup_of(factor[v / 8]);
}
for (int u = 1; u < TS; ++u) {
const u32 *prev = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v += 8)
stu256(cur + v, shoup8(ldu256(prev + v),
factor[v / 8], factor_shoup[v / 8]));
}
}
static void prepare_shared_row_factors(const u32 *RV,__m256i *shared_factor,__m256i *shared_shoup){
for(int v=0;v<TS;v+=8){
shared_factor[v/8]=mont8(ldu256(RV+v),_mm256_set1_epi32(1));
shared_shoup[v/8]=shoup_of(shared_factor[v/8]);
}
}
static void twiddle_rows_prepared(u32 *TW,const __m256i *factor,const __m256i *factor_shoup) {
for (int u = 1; u < TS; ++u) {
const u32 *prev = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v += 8)
stu256(cur + v, shoup8(ldu256(prev + v),
factor[v / 8], factor_shoup[v / 8]));
}
}
static void scale_row_avx(u32 *dst, const u32 *src, u32 k) {
const __m256i kv = _mm256_set1_epi32((int)k), ks = _mm256_set1_epi32((int)shoup_c(k));
for (int v = 0; v < TS; v += 8)
stu256((void *)((dst + v)), shoup8(ldu256((const void *)((src + v))), kv, ks));
}
static void twiddle_rows_avx(u32 *TW, const u32 *RV) {
__m256i RVS[TS / 8];
for (int v = 0; v < TS; v += 8) RVS[v >> 3] = shoup_of(ldu256((const void *)((RV + v))));
for (int u = 1; u < TS; u++) {
const u32 *prv = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v += 8) {
__m256i r = ldu256((const void *)((RV + v)));
stu256((void *)((cur + v)), shoup8(ldu256((const void *)((prv + v))), r, RVS[v >> 3]));
}
}
}
static void corner_tile_avx(const u32 *in, u32 *out, int A0, int B0, int strIn, int strOut, int split) {
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 8) {
__m256i r[8];
for (int k = 0; k < 8; k++) {
const u32 *s = in + (u64)(A0 + u0 + k) * strIn + B0 + v0;
const u32 *tw = TW + (u0 + k) * TS + v0;
__m256i x = ldu256((const void *)(s));
__m256i w = ldu256((const void *)(tw));
r[k] = mont8(x, w);
}
t8x8(r);
for (int k = 0; k < 8; k++)
stu256((void *)(TBUF + (v0 + k) * TS + u0), r[k]);
}
}
for (int v = 0; v < TS; v++) {
u32 *d = split ? BROW((long)(B0 + v) * strOut + A0) : out + (u64)(B0 + v) * strOut + A0;
const u32 *s = TBUF + v * TS;
if ((((uintptr_t)d) & 31u) == 0u) {
for (int u = 0; u < TS; u += 8)
_mm256_stream_si256((__m256i *)(d + u), ldu256((const void *)(s + u)));
} else {
for (int u = 0; u < TS; u += 8) stu256((void *)(d + u), ldu256((const void *)(s + u)));
}
}
}
static void corner_tile_avx_blk(const u32 *in, u32 *out, int A0c, int B0, int strIn, int strOut, int split) {
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 8) {
__m256i r[8];
for (int k = 0; k < 8; k++) {
const u32 *s = in + (u64)(u0 + k) * strIn + B0 + v0;
const u32 *tw = TW + (u0 + k) * TS + v0;
r[k] = mont8(ldu256((const void *)(s)), ldu256((const void *)(tw)));
}
t8x8(r);
for (int k = 0; k < 8; k++)
stu256((void *)(TBUF + (v0 + k) * TS + u0), r[k]);
}
}
for (int v = 0; v < TS; v++) {
u32 *d = split ? BROW((long)(B0 + v) * strOut + A0c) : out + (u64)(B0 + v) * strOut + A0c;
const u32 *s = TBUF + v * TS;
if ((((uintptr_t)d) & 31u) == 0u) {
for (int u = 0; u < TS; u += 8)
_mm256_stream_si256((__m256i *)(d + u), ldu256((const void *)(s + u)));
} else {
for (int u = 0; u < TS; u += 8) stu256((void *)(d + u), ldu256((const void *)(s + u)));
}
}
}
static void corner_tile_avx_to(const u32 *in, u32 *tile, int B0, int strIn) {
(void)strIn;
const long SRF_ = (long)(NR_ + PADB_);
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 8) {
__m256i r[8];
{
const u32 *sb = in + B0 + v0 + (u64)u0 * SRF_;
const u32 *tb = TW + (u64)u0 * TS + v0;
for (int k = 0; k < 8; k++)
r[k] = mont8(ldu256((const void *)(sb + (u64)k * SRF_)),
ldu256((const void *)(tb + (long)k * TS)));
}
t8x8(r);
for (int k = 0; k < 8; k++)
stu256((void *)(tile + (v0 + k) * TS + u0), r[k]);
}
}
}
static void load_tile_avx(const u32 *src, u32 *S, int r0, int c0, int R, int C, int ctr, int allzero) {
for (int rb0 = 0; rb0 < TS; rb0 += 8) {
for (int cb0 = 0; cb0 < TS; cb0 += 16) {
__m256i r[8];
if (allzero) {
for (int k = 0; k < 8; k++) r[k] = _mm256_setzero_si256();
} else {
for (int k = 0; k < 8; k++) {
const u32 *s = src + (u64)(r0 + rb0 + k) * C + c0 + cb0;
__m256i x = ldu256((const void *)(s));
if (ctr) {
__m256i y = _mm256_sub_epi32(x, _mm256_set1_epi32(4));
__m256i z = _mm256_add_epi32(x, _mm256_set1_epi32((int)(P - 4)));
x = _mm256_min_epu32(y, z);
}
r[k] = x;
}
}
t8x8(r);
for (int j = 0; j < 8; j++)
stu256((void *)(TBUF + (cb0 + j) * TS + rb0), r[j]);
{
int cb0b = cb0 + 8;
__m256i r2[8];
if (allzero) {
for (int k = 0; k < 8; k++) r2[k] = _mm256_setzero_si256();
} else {
for (int k = 0; k < 8; k++) {
const u32 *s = src + (u64)(r0 + rb0 + k) * C + c0 + cb0b;
__m256i x = ldu256((const void *)(s));
if (ctr) {
__m256i y = _mm256_sub_epi32(x, _mm256_set1_epi32(4));
__m256i z = _mm256_add_epi32(x, _mm256_set1_epi32((int)(P - 4)));
x = _mm256_min_epu32(y, z);
}
r2[k] = x;
}
}
t8x8(r2);
for (int j = 0; j < 8; j++)
stu256((void *)(TBUF + (cb0b + j) * TS + rb0), r2[j]);
}
}
}
for (int cb = 0; cb < TS; cb++) {
u32 *d = S + (u64)(c0 + cb) * R + r0;
const u32 *s = TBUF + cb * TS;
for (int rb = 0; rb < TS; rb += 8) stu256((void *)(d + rb), ldu256((const void *)(s + rb)));
}
}
static void load_tile_avx_d(const u32 *src, u32 *S, int r0, int c0, int C, int dstr, int ctr, int allzero) {
(void)C; (void)dstr;
const long CF_ = (long)NC_, DF_ = (long)(NR_ + PADB_);
for (int rb0 = 0; rb0 < TS; rb0 += 8) {
for (int cb0 = 0; cb0 < TS; cb0 += 8) {
{ const u32 *pf = src + c0 + cb0 + (u64)(r0 + rb0 + 8) * CF_;
for (int k = 0; k < 8; k++) _mm_prefetch((const char *)(pf + (u64)k * CF_), _MM_HINT_T0); }
__m256i r[8];
if (allzero) {
for (int k = 0; k < 8; k++) r[k] = _mm256_setzero_si256();
} else {
const u32 *sb = src + (u64)(r0 + rb0) * CF_ + c0 + cb0;
for (int k = 0; k < 8; k++) {
__m256i x = ldu256((const void *)(sb + (u64)k * CF_));
if (ctr) {
__m256i y = _mm256_sub_epi32(x, _mm256_set1_epi32(4));
__m256i z = _mm256_add_epi32(x, _mm256_set1_epi32((int)(P - 4)));
x = _mm256_min_epu32(y, z);
}
r[k] = x;
}
}
t8x8(r);
{
const u32 *db = S + (u64)r0 + rb0;
for (int j = 0; j < 8; j++)
stu256((void *)(db + (long)(cb0 + j) * DF_), r[j]);
}
}
}
}
static void load_block(const u32 *src, u32 nsrc, u32 *S, int c0, int R, int C, int dstr, int ctr) {
for (int r0 = 0; r0 < R; r0 += TS) {
u64 firstidx = (u64)r0 * C + c0;
u64 lastidx = (u64)(r0 + TS - 1) * C + c0 + TS - 1;
if (lastidx < nsrc) { load_tile_avx_d(src, S, r0, c0, C, dstr, ctr, 0); continue; }
if (firstidx >= nsrc) { load_tile_avx_d(src, S, r0, c0, C, dstr, ctr, 1); continue; }
for (int rb = 0; rb < TS; rb++) {
u64 base = (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u64 id = base + cb;
u32 v = 0;
if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
TBUF[cb * TS + rb] = v;
}
}
for (int cb = 0; cb < TS; cb++) {
u32 *d = S + (u64)cb * dstr + r0;
const u32 *s = TBUF + cb * TS;
for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
}
}
}
static void copy_u32_avx(u32 *dst, const u32 *src, long n) {
long i = 0;
for (; i + 8 <= n; i += 8) stu256((void *)(dst + i), ldu256((const void *)(src + i)));
for (; i < n; i++) dst[i] = src[i];
}
static void copy_u32_nt(u32 *dst, const u32 *src, long n) {
long i = 0;
long ka = ((-(long)(((uintptr_t)dst) >> 2)) & 7);
for (; i < ka && i < n; i++) dst[i] = src[i];
for (; i + 8 <= n; i += 8) {
__m256i v = ldu256((const void *)(src + i));
_mm256_stream_si256((__m256i *)(dst + i), v);
}
for (; i < n; i++) dst[i] = src[i];
_mm_sfence();
}
static void untr_tile_avx(const u32 *in, u32 *out, int A0, int B0, int strIn, int strOut) {
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 16) {
__m256i r[8];
for (int k = 0; k < 8; k++)
r[k] = ldu256((const void *)(BROW((long)(A0 + u0 + k) * strIn + B0 + v0)));
t8x8(r);
for (int k = 0; k < 8; k++)
stu256((void *)(TBUF + (v0 + k) * TS + u0), r[k]);
{
int v0b = v0 + 8;
__m256i r2[8];
for (int k = 0; k < 8; k++)
r2[k] = ldu256((const void *)(BROW((long)(A0 + u0 + k) * strIn + B0 + v0b)));
t8x8(r2);
for (int k = 0; k < 8; k++)
stu256((void *)(TBUF + (v0b + k) * TS + u0), r2[k]);
}
}
}
for (int v = 0; v < TS; v++) {
u32 *d = out + (u64)(B0 + v) * strOut + A0;
const u32 *s = TBUF + v * TS;
for (int u = 0; u < TS; u += 8)
_mm256_stream_si256((__m256i *)(d + u), ldu256((const void *)(s + u)));
}
}
static void build_tab2(u32 *tw, u32 *tws, int L, u32 w) {
int off = 0;
for (int len = 2; len <= L; len <<= 1) {
u32 wl = pw(w, (u32)(L / len));
u32 c = 1;
int h = len >> 1;
for (int j = 0; j < h; j++) { tw[off + j] = c; tws[off + j] = shoup_c(c); c = mm(c, wl); }
off += h;
}
}
#pragma GCC pop_options
#endif
static u32 R3F[2], R3FS[2];
static u32 R3I[2], R3IS[2];
static u32 T3F_[2][M3_ + 16], T3FS_[2][M3_ + 16];
#define T3F T3F_[0]
#define T3FS T3FS_[0]
static u32 T3I[2 * M3_], T3IS[2 * M3_];
static u32 tsubFw[M3_], tsubFsw[M3_], tsubIw[M3_], tsubIsw[M3_];
static void build_radix3(u32 root) {
const u32 H2 = 1006632961u;
u32 w3 = pw(root, (u32)(NC_ / 3));
u32 w32 = mm(w3, w3);
u32 cA = (u32)(P - H2);
u32 cP = mm((u32)(((u64)w3 + P - w32) % P), H2);
R3F[0] = cA; R3F[1] = cP;
R3I[0] = cA; R3I[1] = (cP ? P - cP : 0);
for (int i = 0; i < 2; i++) { R3FS[i] = shoup_c(R3F[i]); R3IS[i] = shoup_c(R3I[i]); }
for (int s = 1; s < 3; s++) {
u32 g = pw(root, (u32)s), gi = pw(g, P - 2u), c = 1, ci = 1;
u32 *pf = T3F_[s-1], *pf2 = T3FS_[s-1];
u32 *pi = T3I + (s - 1) * M3_, *pi2 = T3IS + (s - 1) * M3_;
for (int j = 0; j < (int)M3_; j++) { pf[j] = c; pf2[j] = shoup_c(c); c = mm(c, g);
pi[j] = ci; pi2[j] = shoup_c(ci); ci = mm(ci, gi); }
}
build_tab2(tsubFw, tsubFsw, (int)M3_, pw(root, 3));
build_tab2(tsubIw, tsubIsw, (int)M3_, pw(pw(root, 3), P - 2u));
}
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __attribute__((always_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));
}
static inline __attribute__((always_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));
}
static inline __attribute__((always_inline)) void dft3(__m256i *v, const __m256i *c,
const __m256i *cs, __m256i Pv) {
__m256i A = vaddm(v[1], v[2], Pv);
__m256i B = _mm256_add_epi32(_mm256_sub_epi32(v[1],v[2]),Pv);
const __m256i parity=_mm256_slli_epi32(A,31);
const __m256i base=_mm256_castps_si256(_mm256_blendv_ps(_mm256_castsi256_ps(Pv),
_mm256_castsi256_ps(_mm256_set1_epi32((int)((P-1u)/2u))),_mm256_castsi256_ps(parity)));
__m256i T=_mm256_sub_epi32(base,_mm256_srli_epi32(A,1));
__m256i U = shoup8_broadcast(B, c[1], cs[1]);
__m256i X0 = vaddm(v[0], A, Pv);
__m256i Z = vaddm(v[0], T, Pv);
v[0] = X0; v[1] = _mm256_add_epi32(Z,U); v[2] = _mm256_add_epi32(_mm256_sub_epi32(Z,U),Pv);
}
static void col_dif3_head(u32 *x) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i c[2], cs[2];
for (int i = 0; i < 2; i++) { c[i] = _mm256_set1_epi32((int)R3F[i]); cs[i] = _mm256_set1_epi32((int)R3FS[i]); }
for (int j = 0; j < (int)M3_; j += 8) {
__m256i v[3];
for (int i = 0; i < 3; i++) v[i] = ldu256((const void *)(x + (u64)i * M3_ + j));
dft3(v, c, cs, Pv);
stu256((void *)(x + j), v[0]);
for (int s = 1; s < 3; s++) {
__m256i w = ldu256((const void *)(T3F_[s-1] + j));
__m256i ws = ldu256((const void *)(T3FS_[s-1] + j));
stu256((void *)(x + (u64)s * M3_ + j), shoup8(v[s], w, ws));
}
}
}
static void col_dif3(u32 *x) {
col_dif3_head(x);
for (int s = 0; s < 3; s++) row_dif_avx(x + (u64)s * M3_, (int)M3_, tsubFw, tsubFsw);
}
static inline __attribute__((always_inline)) void dft3z(__m256i *v, const __m256i *c,
const __m256i *cs, __m256i Pv) {
__m256i A = vaddm(v[1], v[2], Pv);
__m256i B = _mm256_add_epi32(_mm256_sub_epi32(v[1],v[2]),Pv);
const __m256i parity=_mm256_slli_epi32(A,31);
const __m256i base=_mm256_castps_si256(_mm256_blendv_ps(_mm256_castsi256_ps(Pv),
_mm256_castsi256_ps(_mm256_set1_epi32((int)((P-1u)/2u))),_mm256_castsi256_ps(parity)));
__m256i T=_mm256_sub_epi32(base,_mm256_srli_epi32(A,1));
__m256i U = shoup8_broadcast(B, c[1], cs[1]);
const __m256i v0 = v[0];
v[0] = _mm256_add_epi32(v0, A);
__m256i Z = vaddm(v0, T, Pv);
v[1] = _mm256_add_epi32(Z, U);
v[2] = _mm256_add_epi32(_mm256_sub_epi32(Z,U),Pv);
}
#define N02Q_TAILROW(J) do { \
__m256i v[3],w,ws; \
v[0] = ldu256((const void *)(x + (J))); \
for (int s = 1; s < 3; s++) { \
w = ldu256((const void *)(T3I + (s - 1) * M3_ + (J))); \
ws = ldu256((const void *)(T3IS + (s - 1) * M3_ + (J))); \
v[s] = shoup8(ldu256((const void *)(x + (u64)s * M3_ + (J))), w, ws); \
} \
dft3z(v, c, cs, Pv); \
for (int i = 0; i < 3; i++) stu256((void *)(x + (u64)i * M3_ + (J)), v[i]); \
} while (0)
static void col_dit3_tail(u32 *x) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i c[2], cs[2];
for (int i = 0; i < 2; i++) { c[i] = _mm256_set1_epi32((int)R3I[i]); cs[i] = _mm256_set1_epi32((int)R3IS[i]); }
for (int j = 0; j < (int)M3_; j += 16) {
N02Q_TAILROW(j);
N02Q_TAILROW(j + 8);
}
}
static void col_dit3(u32 *x) {
for (int s = 0; s < 3; s++) row_dit_avx(x + (u64)s * M3_, (int)M3_, tsubIw, tsubIsw);
col_dit3_tail(x);
}
static inline __attribute__((always_inline)) void corner_head_store_avx(const u32 *t0, const u32 *t1, const u32 *t2,
u32 *d0, u32 *d1, u32 *d2, int Ab, const __m256i*hc,const __m256i*hcs) {
const __m256i Pvh = _mm256_set1_epi32((int)P);
for (int u = 0; u < TS; u += 8) {
__m256i vv[3];
vv[0] = ldu256((const void *)(t0 + u));
vv[1] = ldu256((const void *)(t1 + u));
vv[2] = ldu256((const void *)(t2 + u));
dft3(vv, hc, hcs, Pvh);
__m256i w1 = ldu256((const void *)(T3F + Ab + u));
__m256i s1 = ldu256((const void *)(T3FS + Ab + u));
__m256i w2 = ldu256((const void *)(T3F_[1] + Ab + u));
__m256i s2 = ldu256((const void *)(T3FS_[1] + Ab + u));
const __m256i out1=shoup8(vv[1],w1,s1),out2=shoup8(vv[2],w2,s2);
if((((uintptr_t)d0|(uintptr_t)d1|(uintptr_t)d2)&31u)==0u){
_mm256_stream_si256((__m256i *)(d0+u),vv[0]);
_mm256_stream_si256((__m256i *)(d1+u),out1);
_mm256_stream_si256((__m256i *)(d2+u),out2);
}else{
stu256(d0+u,vv[0]);stu256(d1+u,out1);stu256(d2+u,out2);
}
}
}
static void corner_head_tile_avx(const u32*tile0,const u32*tile1,const u32*tile2,u32*X,int B0,int Ab,int split){
const int AM=NC_/3,strOut=CSX_;
__m256i hc[2],hcs[2];
for(int i=0;i<2;++i){hc[i]=_mm256_set1_epi32((int)R3F[i]);hcs[i]=_mm256_set1_epi32((int)R3FS[i]);}
for (int v = 0; v < TS; v++) {
const u32 *t0 = tile0 + v * TS, *t1 = tile1 + v * TS, *t2 = tile2 + v * TS;
const u64 rowbase = (u64)(B0 + v) * strOut + Ab;
u32 *d0 = split ? BROW((long)rowbase) : X + rowbase;
u32 *d1 = split ? BROW((long)(rowbase + AM)) : X + rowbase + AM;
u32 *d2 = split ? BROW((long)(rowbase + 2 * (long)AM)) : X + rowbase + 2 * AM;
corner_head_store_avx(t0, t1, t2, d0, d1, d2, Ab,hc,hcs);
}
}
static void col_convolve3(u32 *a,u32 *b) {
for(int arm=0;arm<3;++arm) {
u32 *pa=a+(u64)arm*M3_,*pb=b+(u64)arm*M3_;
row_dif_avx<false,true>(pa,(int)M3_,tsubFw,tsubFsw);
row_dif_avx(pb,(int)M3_,tsubFw,tsubFsw);
if(arm==0)row_dit_avx<true>(pa,(int)M3_,tsubIw,tsubIsw,pb);
else row_dit_avx<true,true>(pa,(int)M3_,tsubIw,tsubIsw,pb);
}
col_dit3_tail(a);
}
#pragma GCC pop_options
static void row_dit(u32 *x, int L, const u32 *tw, const u32 *tws) {
#if HAVE_AVX2
row_dit_avx(x, L, tw, tws);
return;
#else
(void)tws;
#endif
int off = 0;
for (int len = 2; len <= L; len <<= 1) {
int h = len >> 1;
const u32 *t = tw + off; off += h;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = mm(q[j], t[j]);
u32 s = u + v; if (s >= P) s -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = s; q[j] = d;
}
}
}
}
static void pointwise(u32 *A, const u32 *B, u32 n, u32 sc, u32 sx) {
#if HAVE_AVX2
pointwise_mul_avx(A, B, n, sc, sx);
return;
#else
for (long i = 0; i < (long)n; i++) { long o = (i / (long)(sx - PADC_)) * sx + (i % (long)(sx - PADC_)); A[o] = mm(mm(A[o], B[o]), sc); }
#endif
}
static void row_dif(u32 *x, int L, const u32 *tw, const u32 *tws) {
#if HAVE_AVX2
row_dif_avx(x, L, tw, tws);
return;
#else
(void)tws;
#endif
for (int len = L; len >= 2; len >>= 1) {
int h = len >> 1;
const u32 *t = tw + h - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = q[j];
u32 s = u + v; if (s >= P) s -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = s; q[j] = mm(d, t[j]);
}
}
}
}
static inline void bitrev_pass(u32 *x, int L, const u16 *br) {
for (int i = 0; i < L; i++) {
int k = br[i];
if (i < k) { u32 t = x[i]; x[i] = x[k]; x[k] = t; }
}
}
struct Eng {
int R, C, N, sx;
int blkr;
int r5;
u32 w;
u32 *tfR; u32 *tiR; u32 *tfC; u32 *tiC;
u32 *wsR; u32 *wsC;
u16 *brR; u16 *brC;
};
static void col_dif3_arms(u32 *x) {
for (int s = 0; s < 3; s++) row_dif_avx(x + (u64)s * M3_, (int)M3_, tsubFw, tsubFsw);
}
static void colfwd(Eng &e, u32 *row) {
#if HAVE_AVX2
if (e.r5) { col_dif3_arms(row); return; }
#endif
row_dif(row, e.C, e.tfC, e.wsC);
}
static void colinv(Eng &e, u32 *row) {
#if HAVE_AVX2
if (e.r5) { col_dit3(row); return; }
#endif
row_dit(row, e.C, e.tiC, e.wsC + e.C);
}
static void load_T(Eng &e, const u32 *src, u32 nsrc, u32 *S, int ctr) {
const int R = e.R, C = e.C;
#if HAVE_AVX2
if (nsrc <= (u32)R * (u32)C) {
for (int r0 = 0; r0 < R; r0 += TS) {
u64 tilerow = (u64)r0 * C;
for (int c0 = 0; c0 < C; c0 += TS) {
u64 firstidx = tilerow + c0;
u64 lastidx = (u64)(r0 + TS - 1) * C + c0 + TS - 1;
if (lastidx < nsrc) { load_tile_avx(src, S, r0, c0, R, C, ctr, 0); continue; }
if (firstidx >= nsrc) { load_tile_avx(src, S, r0, c0, R, C, ctr, 1); continue; }
for (int rb = 0; rb < TS; rb++) {
u64 base = (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u64 id = base + cb;
u32 v = 0;
if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
TBUF[cb * TS + rb] = v;
}
}
for (int cb = 0; cb < TS; cb++) {
u32 *d = S + (u64)(c0 + cb) * R + r0;
const u32 *s = TBUF + cb * TS;
for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
}
}
}
return;
}
#endif
for (int r0 = 0; r0 < R; r0 += TS) {
for (int c0 = 0; c0 < C; c0 += TS) {
u64 lastidx = (u64)(r0 + TS - 1) * C + (c0 + TS - 1);
if (lastidx < nsrc) {
for (int rb = 0; rb < TS; rb++) {
const u32 *s = src + (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u32 v = s[cb];
if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u));
TBUF[cb * TS + rb] = v;
}
}
} else {
for (int rb = 0; rb < TS; rb++) {
u64 base = (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u64 id = base + cb;
u32 v = 0;
if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
TBUF[cb * TS + rb] = v;
}
}
}
for (int cb = 0; cb < TS; cb++) {
u32 *d = S + (u64)(c0 + cb) * R + r0;
const u32 *s = TBUF + cb * TS;
for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
}
}
}
}
static u32 REV6[TS];
static u32 g_inv_scale = 1u;
static u32 POWT[1024], POWT64[1024];
static void build_rev6(void) {
for (int i = 0; i < TS; i++) {
int r = 0;
for (int b = 0; b < 6; b++) if (i & (1 << b)) r |= 1 << (5 - b);
REV6[i] = (u32)r;
}
}
static inline u32 revbits(u32 x, int lg) {
u32 r = 0;
for (int b = 0; b < lg; b++) if (x & (1u << b)) r |= 1u << (lg - 1 - b);
return r;
}
template<int perm>
static void corner_generic(const u32 *in, u32 *out, int nA, int nB, int strIn, int strOut, u32 root, int split) {
u32 wv[TS];
{ u32 c = 1; for (int v = 0; v < TS; v++) { wv[v] = c; c = mm(c, root); } }
const u32 w64 = pw(root, TS);
int lgB = 0, lgA = 0;
while ((1 << lgB) < nB) lgB++;
while ((1 << lgA) < nA) lgA++;
const int S1 = nB / TS, S2 = nA / TS;
const int HB = lgB - 6, HA = lgA - 6;
u32 RV2[TS], RV2b[TS], PW2[TS], cu[TS], cu64[TS], base[TS];
u32 cu64_plain[TS],cu64_shoup[TS];
if (perm == 1) {
u32 g = pw(root, (u32)S1);
{ u32 gM = mm(g, R32), c = R32; for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; c = mont1(c, gM); } }
{ u32 rootM = mm(root, R32), p = R32; for (int t = 0; t < S1; t++) { POWT[t] = p; p = mont1(p, rootM); } }
} else if (perm == 2) {
u32 g = pw(root, (u32)S2), g64 = pw(g, 64);
const u32 rootM2 = mm(root, R32), gM = mm(g, R32), g64M = mm(g64, R32);
{ u32 c = R32, c64 = R32;
for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; RV2b[REV6[t]] = c64; c = mont1(c, gM); c64 = mont1(c64, g64M); } }
{ u32 p = R32; for (int t = 0; t < S2; t++) { POWT[t] = p; p = mont1(p, rootM2); } }
{ u32 q = pw(root, 64), qM = mm(q, R32), p = R32; for (int t = 0; t < S2; t++) { POWT64[t] = p; p = mont1(p, qM); } }
}
u32 wr = 1;
const u32 INVSCALE = g_inv_scale;
for (int A0 = 0; A0 < nA; A0 += TS) {
u32 wr64 = pw(wr, TS);
if (perm == 1) {
u32 g = pw(wr, (u32)S1), gM2 = mm(g, R32), c = R32;
for (int t = 0; t < TS; t++) { PW2[REV6[t]] = c; c = mont1(c, gM2); }
} else if (perm == 2) {
int hi = A0 / TS;
u32 rH = POWT[revbits((u32)hi, HA)], r64H = POWT64[revbits((u32)hi, HA)];
const u32 INVSCALE_M = mm(INVSCALE, R32);
for (int u = 0; u < TS; u++) {
cu[u] = mont1(RV2[u], rH); cu64[u] = mont1(RV2b[u], r64H); base[u] = INVSCALE_M;
cu64_plain[u]=mont1(cu64[u],1u); cu64_shoup[u]=shoup_c(cu64_plain[u]);
}
} else {
u32 c = 1; for (int v = 0; v < TS; v++) { PW2[v] = c; c = mm(c, wr); }
}
u32 Vloc = 1;
u32 Vr = 1;
for (int B0 = 0; B0 < nB; B0 += TS) {
if (perm == 1) {
u32 H = revbits((u32)(B0 / TS), HB);
u32 wrH = pw(wr, H), rootH = POWT[H];
u32 RV[TS];
#if HAVE_AVX2
scale_row_mont(RV, RV2, rootH);
scale_row_mont(TW, PW2, mm(wrH, R32));
twiddle_rows_mont(TW, RV);
#else
for (int v = 0; v < TS; v++) RV[v] = mont1(RV2[v], rootH);
for (int v = 0; v < TS; v++) TW[v] = mont1(PW2[v], mm(wrH, R32));
for (int u = 1; u < TS; u++) {
const u32 *prv = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v++) cur[v] = mont1(prv[v], RV[v]);
}
#endif
} else if (perm == 2) {
if (B0 == 0) {
for (int u = 0; u < TS; u++) {
u32 *cur = TW + u * TS;
cur[0] = base[u];
for (int v = 1; v < TS; v++) cur[v] = mont1(cur[v - 1], cu[u]);
}
} else {
for (int u = 0; u < TS; u++) {
u32 *cur = TW + u * TS;
#if HAVE_AVX2
scale_row_prepared(cur, cur, cu64_plain[u], cu64_shoup[u]);
#else
u32 d = cu64[u];
for (int v = 0; v < TS; v++) cur[v] = mont1(cur[v], d);
#endif
}
}
} else {
u32 RV[TS];
for (int v = 0; v < TS; v++) RV[v] = mm(Vloc, wv[v]);
for (int v = 0; v < TS; v++) TW[v] = mm(Vr, PW2[v]);
for (int u = 1; u < TS; u++) {
const u32 *prv = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v++) cur[v] = mm(prv[v], RV[v]);
}
}
#if HAVE_AVX2
corner_tile_avx(in, out, A0, B0, strIn, strOut, split);
#else
for (int u = 0; u < TS; u++) {
const u32 *s = in + (u64)(A0 + u) * strIn + B0;
const u32 *tw = TW + u * TS;
for (int v = 0; v < TS; v++) TBUF[u * TS + v] = mont1(s[v], tw[v]);
}
for (int v = 0; v < TS; v++) {
u32 *d = out + (u64)(B0 + v) * strOut + A0;
for (int u = 0; u < TS; u++) d[u] = TBUF[u * TS + v];
}
#endif
Vr = mm(Vr, wr64);
Vloc = mm(Vloc, w64);
}
wr = mm(wr, w64);
}
}
static inline void corner_fwd(Eng &e, const u32 *S, u32 *X, u32 root) {
corner_generic<1>(S, X, e.C, e.R, e.R, e.sx, root, 0);
}
static inline void corner_inv(Eng &e, const u32 *X, u32 *S, u32 root) {
corner_generic<2>(X, S, e.R, e.C, e.sx, e.R, root, 1);
}
static void untranspose(Eng &e, const u32 *S, u32 *dst) {
const int R = e.R, C = e.C;
#if HAVE_AVX2
for (int c0 = 0; c0 < C; c0 += TS)
for (int r0 = 0; r0 < R; r0 += TS)
untr_tile_avx(S, dst, c0, r0, R, e.sx);
return;
#endif
for (int c0 = 0; c0 < C; c0 += TS) {
for (int r0 = 0; r0 < R; r0 += TS) {
for (int u = 0; u < TS; u++) {
const u32 *s = S + (u64)(c0 + u) * R + r0;
u32 *d = TBUF + u * TS;
for (int v = 0; v < TS; v++) d[v] = s[v];
}
for (int v = 0; v < TS; v++) {
u32 *d = dst + (u64)(r0 + v) * e.sx + c0;
for (int u = 0; u < TS; u++) d[u] = TBUF[u * TS + v];
}
}
}
}
#define N02I_SB (TS * (NRMAX_ + PADB_))
static u32 SBLK[3 * N02I_SB] __attribute__((aligned(64)));
static u32 TILE3[3][TS * TS + 32] __attribute__((aligned(64)));
static void corner_fwd_stream(Eng &e, const u32 *src, u32 nsrc, u32 *X, int split, int ctr) {
const int nA = (int)NC_, nB = (int)NR_, strIn = (int)(NR_ + PADB_), strOut = (int)CSX_;
const u32 root = e.w;
int lgB = 0; while ((1 << lgB) < nB) lgB++;
const int S1 = nB / TS, HB = lgB - 6;
u32 RV2[TS], PW2[3][TS], WRHM[3][NR_/TS];
{ u32 g = pw(root, (u32)S1), gM = mm(g, R32), c = R32;
for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; c = mont1(c, gM); } }
{ u32 rootM = mm(root, R32), p = R32;
for (int t = 0; t < S1; t++) { POWT[t] = p; p = mont1(p, rootM); } }
const bool halfzero=(nB==4096 && (u64)nsrc<=(u64)(nB/2)*nA);
const int AM = nA / 3;
for (int Ab = 0; Ab < AM; Ab += TS) {
for (int p = 0; p < 3; p++) {
u32 *S = SBLK + (u64)p * N02I_SB;
const int A0 = Ab + p * AM;
const u32 wr = pw(root, (u32)A0);
{ const u32 wrM=mm(wr,R32); u32 val=R32;
for(int h=0;h<S1;++h){WRHM[p][h]=val; val=mont1(val,wrM);}
}
{ u32 g = pw(wr, (u32)S1), gM2 = mm(g, R32), c = R32;
for (int t = 0; t < TS; t++) { PW2[p][REV6[t]] = c; c = mont1(c, gM2); } }
load_block(src, nsrc, S, A0, halfzero?nB/2:nB, nA, strIn, ctr);
for (int u=0;u<TS;++u){
if(halfzero)row_dif_avx<true,true>(S+(u64)u*strIn,nB,e.tfR,e.wsR);
else row_dif(S+(u64)u*strIn,nB,e.tfR,e.wsR);
}
}
for (int B0 = 0; B0 < nB; B0 += TS) {
u32 H = revbits((u32)(B0 / TS), HB);
u32 rootH = POWT[H];
u32 RV[TS];
scale_row_mont(RV,RV2,rootH);
__m256i shared_factor[TS/8],shared_shoup[TS/8];
prepare_shared_row_factors(RV,shared_factor,shared_shoup);
for (int p = 0; p < 3; p++) {
scale_row_mont(TW, PW2[p], WRHM[p][H]);
twiddle_rows_prepared(TW, shared_factor, shared_shoup);
corner_tile_avx_to(SBLK + (u64)p * N02I_SB, TILE3[p], B0, strIn);
}
corner_head_tile_avx(TILE3[0],TILE3[1],TILE3[2],X,B0,Ab,split);
}
}
}
static void fwd(Eng &e, const u32 *src, u32 nsrc, u32 *S, u32 *X, int split, int ctr) {
#ifdef OLD_FWD
load_T(e, src, nsrc, S, ctr);
for (int j0 = 0; j0 < e.C; j0++) {
u32 *row = S + (u64)j0 * e.R;
row_dif(row, e.R, e.tfR, e.wsR);
}
corner_fwd(e, S, X, e.w);
#else
(void)S;
{ corner_fwd_stream(e, src, nsrc, X, split, ctr); }
#endif
if (!g_fuse4) { for (int k0 = 0; k0 < e.R; k0++)
colfwd(e, split ? BROW((long)k0 * e.sx) : X + (u64)k0 * e.sx); }
}
#if 0
static void fwd_old(Eng &e, const u32 *src, u32 nsrc, u32 *S, u32 *X, int ctr) {
load_T(e, src, nsrc, S, ctr);
for (int j0 = 0; j0 < e.C; j0++) {
u32 *row = S + (u64)j0 * e.R;
row_dif(row, e.R, e.tfR, e.wsR);
}
corner_fwd(e, S, X, e.w);
for (int k0 = 0; k0 < e.R; k0++) colfwd(e, X + (u64)k0 * e.sx);
}
#endif
static void inv(Eng &e, u32 *X, u32 *S) { for (int k0 = 0; k0 < e.R; k0++) colinv(e, X + (u64)k0 * e.sx);
corner_inv(e, X, S, pw(e.w, P - 2u));
for (int j0 = 0; j0 < e.C; j0++) {
u32 *row = BROW((long)j0 * e.R);
row_dit(row, e.R, e.tiR, e.wsR + e.R);
}
}
static u32 INV_TW[NRMAX_ * TS] __attribute__((aligned(64)));
static u32 INV_STEP[NRMAX_], INV_STEP_SHOUP[NRMAX_];
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
template<bool Advance>
static inline __attribute__((always_inline)) void inverse_twiddle_consumer(Eng &e, u32 *X, int r0, int c0) {
(void)e;
const long SXF_ = (long)CSX_, BLF_ = (long)(NR_ + PADB_);
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 8) {
{ const u32 *pf = X + (u64)(r0 + u0 + 8) * SXF_ + c0 + v0;
for (int k = 0; k < 8; k++) _mm_prefetch((const char *)(pf + (u64)k * SXF_), _MM_HINT_T0); }
__m256i rows[8];
{
const u32 *sbase = X + (u64)(r0 + u0) * SXF_ + c0 + v0;
u32 *tbase = INV_TW + (r0 + u0) * (long)TS + v0;
u32 *sblkbase = SBLK + (u64)v0 * BLF_ + r0 + u0;
for (int k = 0; k < 8; ++k) {
const u32 *src = sbase + (u64)k * SXF_;
u32 *tw = tbase + (long)k * TS;
__m256i w = ldu256(tw);
if constexpr(Advance) {
const int j = r0 + u0 + k;
stu256(tw, shoup8_broadcast(w, _mm256_set1_epi32((int)INV_STEP[j]),
_mm256_set1_epi32((int)INV_STEP_SHOUP[j])));
}
rows[k] = mont8(ldu256(src), w);
}
t8x8(rows);
for (int k = 0; k < 8; ++k)
stu256(sblkbase + (u64)k * BLF_, rows[k]);
}
}
}
}
static void inv_natural_stream(Eng &e, u32 *X) {
if (!g_fuse4)
for (int k = 0; k < e.R; ++k) colinv(e, X + (u64)k * e.sx);
const u32 root = pw(e.w, P - 2u);
const u32 seed = mm(mm(g_inv_scale, R32), R32);
int lg = 0; while ((1 << lg) < e.R) ++lg;
for (int k = 0; k < e.R; ++k) {
const u32 factor = pw(root, revbits(k, lg));
const u32 factorM = mm(factor, R32);
INV_STEP[k] = pw(factor, TS);
INV_STEP_SHOUP[k] = shoup_c(INV_STEP[k]);
u32 *tw = INV_TW + k * TS;
tw[0] = seed;
for (int v = 1; v < TS; ++v) tw[v] = mont1(tw[v - 1], factorM);
}
for (int c0 = 0; c0 < e.C; c0 += TS) {
for (int r0 = 0; r0 < e.R; r0 += TS) {
inverse_twiddle_consumer<true>(e,X,r0,c0);
}
for (int v = 0; v < TS; ++v)
row_dit(SBLK + (u64)v * e.blkr, e.R, e.tiR, e.wsR + e.R);
{
const long SXF_ = (long)CSX_, BLF_ = (long)(NR_ + PADB_);
for (int r0 = 0; r0 < (int)NR_; r0 += TS) {
for (int v0 = 0; v0 < TS; v0 += 8) {
{
const u32 *sbase = SBLK + (u64)(r0 + v0);
u32 *dbase = X + (u64)(r0 + v0) * SXF_ + c0;
for (int u0 = 0; u0 < TS; u0 += 16) {
const u32 *sp = sbase + (u64)u0 * BLF_;
__m256i rows[8];
for (int k = 0; k < 8; ++k)
rows[k] = ldu256(sp + (u64)k * BLF_);
t8x8(rows);
for (int k = 0; k < 8; ++k) stu256(TBUF+k*8,rows[k]);
for (int k = 0; k < 8; ++k)
rows[k] = ldu256(sp + (u64)(8 + k) * BLF_);
t8x8(rows);
for (int k = 0; k < 8; ++k) {
u32 *dst = dbase + (u64)k * SXF_ + u0;
_mm256_stream_si256((__m256i *)dst, ldu256(TBUF+k*8));
_mm256_stream_si256((__m256i *)(dst+8), rows[k]);
}
}
}
}
}
}
}
_mm_sfence();
}
static void fuse4_row_avx(u32 *pa, const u32 *pb, int n) {
int i = 0;
for (; i + 8 <= n; i += 8) {
__m256i x = ldu256((const void *)(pa + i));
__m256i y = ldu256((const void *)(pb + i));
stu256((void *)(pa + i), shoup8(x, y, shoup_of(y)));
}
for (; i < n; i++) pa[i] = mm(pa[i], pb[i]);
}
#pragma GCC pop_options
static void eng_init(Eng &e, int lgR, int nC, u32 w, int sx,
u32 *tfR, u32 *tiR, u32 *tfC, u32 *tiC, u16 *brR, u16 *brC,
u32 *wsR, u32 *wsC) {
e.R = 1 << lgR; e.C = nC; e.N = e.R * e.C; e.w = w; e.sx = sx;
e.blkr = e.R + PADB_;
e.r5 = (nC == (int)NC_);
e.tfR = tfR; e.tiR = tiR; e.tfC = tfC; e.tiC = tiC; e.brR = brR; e.brC = brC;
e.wsR = wsR; e.wsC = wsC;
u32 wR = pw(w, (u32)e.C), wC = pw(w, (u32)e.R);
#if HAVE_AVX2
build_tab2(tfR, e.wsR, e.R, wR); build_tab2(tiR, e.wsR + e.R, e.R, pw(wR, P - 2u));
if (!e.r5) { build_tab2(tfC, e.wsC, e.C, wC); build_tab2(tiC, e.wsC + e.C, e.C, pw(wC, P - 2u)); }
#else
build_tab(tfR, e.R, wR); build_tab(tiR, e.R, pw(wR, P - 2u));
build_tab(tfC, e.C, wC); build_tab(tiC, e.C, pw(wC, P - 2u));
#endif
build_br(brR, e.R);
if (!e.r5) build_br(brC, e.C);
build_rev6();
}
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __attribute__((always_inline)) __m256i scan8(__m256i v, __m256i &carry) {
v = _mm256_add_epi32(v, _mm256_slli_si256(v, 4));
v = _mm256_add_epi32(v, _mm256_slli_si256(v, 8));
__m256i t = _mm256_permute2x128_si256(_mm256_setzero_si256(), v, 0x20);
t = _mm256_shuffle_epi32(t, 0xFF);
v = _mm256_add_epi32(v, t);
v = _mm256_add_epi32(v, carry);
carry = _mm256_permutevar8x32_epi32(v, _mm256_set1_epi32(7));
return v;
}
#define CMBSTEP() \
{ u32 d; \
if (PH == 0) d = 4u * (a[k] + b[k]) - 16u; \
else if (PH == 1) { if (nm_lt) d = 4u * (b[k] - b[k - 1 - n]); \
else d = 4u * (a[k] - a[k - 1 - m]); } \
else d = 16u - 4u * (a[k - 1 - m] + b[k - 1 - n]); \
SUMr += d; \
u32 tab; \
if (TB == 1) { u32 x = A[AOFF(k)], y = c[N1_ + k]; tab = (x >= y) ? (x - y) : (x + (P - y)); } \
else if (TB == 2) tab = A[AOFF(k)]; \
else tab = c[k]; \
i32 ts = (i32)tab; \
if (ts > (i32)(P >> 1)) ts -= (i32)P; \
c[k] = (u32)((u32)ts + SUMr); \
}
template<int PH, int TB, int NT>
static void cmbv(long k0, long k1, const u32 *a, const u32 *b, u32 *A, u32 *c,
long n, long m, long nmax, u32 &SUMr, int nm_lt)
{
if (k0 > k1) return;
long k = k0;
{
long ka = k + ((-(long)k) & 7);
long kend = (ka < k1 + 1) ? ka : k1 + 1;
if (ka > k) { const __m256i Pv0 = _mm256_set1_epi32((int)P); (void)Pv0; }
for (; k < kend; k++) { CMBSTEP(); }
}
const __m256i Pv = _mm256_set1_epi32((int)P);
const __m256i Phv = _mm256_set1_epi32((int)(P >> 1));
#if HAVE_AVX2
const __m256i c16 = _mm256_set1_epi32(16);
const __m256i idx7 = _mm256_set1_epi32(7);
__m256i carry = _mm256_set1_epi32((int)SUMr);
if (NT) {
carry = _mm256_set1_epi32((int)SUMr);
}
for (; k + 7 <= k1; k += 8) {
__m256i d;
if (PH == 0) {
__m256i va = ldu256(a + k), vb = ldu256(b + k);
d = _mm256_sub_epi32(_mm256_slli_epi32(_mm256_add_epi32(va, vb), 2), c16);
} else if (PH == 1) {
if (nm_lt) { __m256i v1 = ldu256(b + k), v0 = ldu256(b + k - 1 - n);
d = _mm256_slli_epi32(_mm256_sub_epi32(v1, v0), 2); }
else { __m256i v1 = ldu256(a + k), v0 = ldu256(a + k - 1 - m);
d = _mm256_slli_epi32(_mm256_sub_epi32(v1, v0), 2); }
} else {
__m256i va = ldu256(a + k - 1 - m), vb = ldu256(b + k - 1 - n);
d = _mm256_sub_epi32(c16, _mm256_slli_epi32(_mm256_add_epi32(va, vb), 2));
}
d = scan8(d, carry);
__m256i tab;
if (TB == 1) {
__m256i x = ldu256(A + AOFF(k)), y = ldu256(c + N1_ + k);
__m256i dd = _mm256_sub_epi32(x, y);
__m256i mk = _mm256_cmpgt_epi32(y, x);
tab = _mm256_add_epi32(dd, _mm256_and_si256(mk, Pv));
} else if (TB == 2) {
tab = ldu256(A + AOFF(k));
} else {
tab = ldu256(c + k);
}
__m256i ts = _mm256_sub_epi32(tab, _mm256_and_si256(_mm256_cmpgt_epi32(tab, Phv), Pv));
__m256i out = _mm256_add_epi32(ts, d);
if (NT && (((uintptr_t)(c + k)) & 31) == 0) _mm256_stream_si256((__m256i *)(c + k), out);
else stu256(c + k, out);
}
SUMr = (u32)_mm256_extract_epi32(carry, 0);
#endif
for (; k <= k1; k++) { CMBSTEP(); }
}
#pragma GCC pop_options
static u32 Aarr[NR_ * CSX_] __attribute__((aligned(64)));
#if (NR_ * CSX_) > BSTAIL_SPLIT
#define BTAILN_ ((int)((NR_ * CSX_) - BSTAIL_SPLIT))
#else
#define BTAILN_ 1
#endif
static u32 BTAIL[BTAILN_] __attribute__((aligned(64)));
static Eng E1, E2;
static u32 t1fR[NR_], t1iR[NR_], t1fC[1], t1iC[1];
static u16 b1R[NR_], b1C[1];
static u32 w1R[2 * NR_], w1C[2];
static u32 t2fR[NR2_], t2iR[NR2_], t2fC[NC2_], t2iC[NC2_];
static u16 b2R[NR2_], b2C[NC2_];
static u32 w2R[2 * NR2_], w2C[2 * NC2_];
static int inited = 0;
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
IM = ~0ull / (u64)P;
if (!inited) {
inited = 1;
{ u64 inv = 1; for (int i = 0; i < 5; i++) inv = inv * (2 - (u64)P * inv); NPI = (u32)(-(u64)inv); }
u32 w1 = pw(G, (u32)((P - 1) / N1_));
u32 w2 = pw(G, (u32)((P - 1) / N2_));
eng_init(E1, LGR_, (int)NC_, w1, (int)CSX_, t1fR, t1iR, t1fC, t1iC, b1R, b1C, w1R, w1C);
eng_init(E2, LGR2_, (1 << LGC2_), w2, (int)CSX2_, t2fR, t2iR, t2fC, t2iC, b2R, b2C, w2R, w2C);
E1.r5 = 1; E2.r5 = 0;
#if HAVE_AVX2
build_radix3(pw(w1, (u32)E1.R));
#endif
}
g_B = (u32 *)c; g_BT = BTAIL;
u32 *S = g_B;
u32 *A = Aarr, *B = g_B;
const int ctr = EXACT;
const u32 sc = pw((u32)(N1_ % P), P - 2u);
const u32 sc2 = pw((u32)(N2_ % P), P - 2u);
const long Sa = (long)N1_ - (long)m;
const long Sb = (long)N1_ - (long)n;
const long M = (long)n + (long)m - (long)N1_;
if (M >= 0) {
g_inv_scale = sc2;
fwd(E2, (const u32 *)a + Sa, (u32)(M + 1), S, A, 0, ctr);
fwd(E2, (const u32 *)b + Sb, (u32)(M + 1), S, B, 1, ctr);
pointwise(A, B, N2_, sc2, CSX2_);
inv(E2, A, S);
untranspose(E2, S, A);
{
long k = M, e = M + M + 1;
u32 *d = c + N1_;
while (k < e) {
long col = k % (long)NC2_;
long run = (long)NC2_ - col;
if (run > e - k) run = e - k;
const u32 *s = A + (k / (long)NC2_) * CSX2_ + col;
if ((d - c) % 8 == 0 && col % 8 == 0 && run % 8 == 0 && run > 0) copy_u32_nt(d, s, run);
else for (long i = 0; i < run; i++) d[i] = s[i];
d += run; k += run;
}
}
}
g_inv_scale = sc;
g_fuse4 = 1;
fwd(E1, (const u32 *)a, (u32)n + 1, S, A, 0, ctr);
fwd(E1, (const u32 *)b, (u32)m + 1, S, B, 1, ctr);
for (int r = 0; r < (int)NR_; r++) {
u32 *pa = A + (u64)r * CSX_;
u32 *pb = BROW((long)r * CSX_);
if(E1.r5) col_convolve3(pa,pb);
else {
colfwd(E1, pa);
colfwd(E1, pb);
fuse4_row_avx(pa, pb, (int)NC_);
colinv(E1, pa);
}
}
inv_natural_stream(E1, A);
g_fuse4 = 0;
long nmax = (long)n + (long)m;
#if EXACT
{
const long kA = ((long)n < (long)m) ? (long)n : (long)m;
const long kB = ((long)n < (long)m) ? (long)m : (long)n;
const int nm_lt = ((long)n < (long)m);
const long M1 = (M + 1 > 0) ? (M + 1) : 0;
u32 SUM = 0;
cmbv<0,1,1>(0, (M < kA) ? M : kA,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
cmbv<0,2,1>(M1, (N1_ - 1 < kA) ? (long)(N1_ - 1) : kA,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
cmbv<0,3,0>((long)N1_, kA,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
cmbv<1,2,1>((kA + 1 > M1) ? kA + 1 : M1, (N1_ - 1 < kB) ? (long)(N1_ - 1) : kB,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
cmbv<1,3,0>((kA + 1 > (long)N1_) ? kA + 1 : (long)N1_, kB,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
cmbv<2,2,1>((kB + 1 > M1) ? kB + 1 : M1, (N1_ - 1 < nmax) ? (long)(N1_ - 1) : nmax,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
cmbv<2,3,0>((kB + 1 > (long)N1_) ? kB + 1 : (long)N1_, nmax,
a, b, Aarr, c, (long)n, (long)m, nmax, SUM, nm_lt);
}
#else
for (long k = 0; k <= nmax; k++) c[k] = c[k];
#endif
_mm_sfence();
(void)E2;
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 3.371 s | 1584 MB + 476 KB | Accepted | Score: 100 | 显示更多 |