代码

// ===== REFERENCES（引用段）=====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#129718** <https://duck.ac/submission/129718>
//     （判题机 3349.773463 ms = 取底时本题排行榜第一 = 实时 T）。
//     用途：**本文件正文 = 该件正文的逐字节底盘**（经我方在册件 #130554 的逐字吸收件转承，
//     其自带的全上游引用链 / 署名 / 许可声明原样保留在下方），本发在其上只加**一处单变量**。
// [2] 本队内部模型账 `problems/_coord/1002e8_depchain_forms.md`（席 `r136a_`，2026-10-16 19:19 CST）
//     —— 本发单变量的**唯一依据**：该账对 `row_dif_avx<true,true>` 热环（判题同款 gcc 9.3 件）
//     实测 `sim 101.10` / `max 单口 82.00(p5)` / 残余 `19.10`，并给出敏感度矩阵。
// 许可遵守：duck.ac 站内公开提交（站内明示提交正文公开可见），未见单独许可证声明
//     ⇒ 按「站内公开提交 + 作者账号 + 提交 URL + 所用内容」引用。
// ===== 思路 =====
// 【本发性质】**单变量 · 试验性提交（目的 = 定价，不作「转绿」声称）**
//   单变量 = **常数寄存器化**：模板 `row_dif_avx<bool,bool>` 内
//   `__m256i Pv = _mm256_set1_epi32((int)P);` 之后加**一行** `asm("" : "+x"(Pv));`
//   （空 asm + `"+x"` 输出约束 ⇒ 禁止 GCC 把 P 常数再物化进循环体；无指令发出、无副作用）。
//   **机器码级自证**（判题同款 gcc 9.3 · `-S`）：热函数 `row_dif_avx<true,true>` 内
//   折叠常数 **36 → 16**，其中 `vpaddd .LC(%rip)` 折叠 **28 → 8**（该热环里 20 处全消）；
//   `.text` 其余差异逐节未变（本发只改这一处源码行）。
//   ★ 诚实登记：**同源第二处**（`shoup8h` 内的 `const Pv`）本席**试过并撤回** ——
//   加上后热函数折叠反而回到 34（寄存器压力外溢）⇒ **本发只含一处**，余 4 处 `vpmuludq .LC4` 折叠**未覆盖**。
// 【预测区间（写死，发前定标）】**`[−14.8%, 0]`**（模型口径；`−14.8%` = 24 处全寄存器化的模型值，
//   本发只覆盖 20/24 ⇒ 实际预期更靠近区间下沿的 **`≈−12%`~`0`**）。**达标所需 Δ = −1.095%**（缺口 36.723795 ms）。
//   ★ **【模型·未标定】**（`律 1612`）：本环**已在溢出**（`(%rsp)/32/64/96(%rsp)` 读写）⇒ 多 2 个活跃寄存器
//   可能引发新溢出 ⇒ 模型值是**上界**，**判题机单价无实测 ⇒ 本预测只作相对排序**。
// 【内存轴三数（`律 1626`/`律 1631`）】① 表体积变化 = **0**（常数向量本就在 `.rodata`，2×32 B）
//   ② 每迭代 L1 装载数 = **−24**（去掉 24 次折叠 load-op 的装载；★方向 = **减装载**，与在册判负的
//   `[SHQH]`「ALU→load」**反号**）③ 驻留档 = **L1 不变**（2×32 B）。
// 【值闸门（`律 1625`：只声称闸门能给的）】同一 driver（`problems/1002e8/work/r135a__tmp/b_iso.cpp` 的
//   判题形状 driver）跑底盘与本发：
//   ① **输出 `chk`/`fnv` 逐位同**（`n=1e8 pat=0` 判题形状 + `n=2e6 pat=1/3` 形状）；
//   ② **独立 oracle**（driver 自带：奇点 Horner 环同态 `C(X)=A(X)B(X)` 6 点 + 小 n 逐系数朴素卷积）
//      `bad=0` · **`GATE=on`**；
//   ③ **正控有牙**：把常数改成 `P+1` 的同形件 ⇒ `SPOTCHECK bad=6 FAIL` · `SMALL bad=8193 FAIL` · **`GATE=off`** ✓
// ======================
/* ===== REFERENCES =====
 * [1] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#129718**
 *     <https://duck.ac/submission/129718> —— 本题榜面 rank 1（判题机 3349.773463 ms = 本题现役 T）。
 *     **本文件正文 = 该提交正文的逐字节副本**（含其全部上游引用链、署名与许可声明，原样保留在下方），
 *     仅本席标注的一处改动（见下）。本席未新增任何上游引用。
 * [2] duck.ac 题目页 <https://duck.ac/problem/1002e8>（提交正文按站点规则公开可见）。
 * LICENSE：上述为站点公开提交，无独立许可证声明、非来自任何外部仓库 ⇒ 按具名引用
 *     （账号 + 提交地址 + 沿用内容）使用。按 /home/yjp/duck.ac/RULES.md 第 3 节执行。
 * ===== 思路 =====
 * 性质 = **试验性提交（定价）**：本发只做一件事 —— 量【端口再平衡】在本题的兑现。
 * 【单变量 `[SHQH]`】把 Shoup 商常数表的**高半拆位**搬出内环：原状每次 `shoup8` 都要对**常数**做一条
 *   `_mm256_bsrli_epi128(ws,4)`（SKL = `vpsrldq`，**p5**）；本发在建表期另存拆位结果（`*_h` 表），
 *   内环改为**多一条 load**（**p2/p3**）。**逐位等价**：`wsh[i]=(i%4==3)?0:ws[i+1]` 恰等于该位移的逐元素值。
 * 【前置读数（本席 `work/r135a__tmp/`，`llvm-mca -mcpu=skylake`）】`dif_pass8<64>` 环 76.05 拍/迭代、
 *   `p0=72.9 p1=72.0 p5=68.1` 而 `p2=p3=18`（取数口空着）⇒ 12 条 p5 位移换 12 条 load，模型 **−6.4%**；
 *   受影响环 = `dif_pass8` 族 28.3% + `dit_pass8` 族 16.6% Ir ⇒ 名义 ≈ 90 ms（在册 `律 1537` 保守 ≈ 67 ms）。
 * 【预测区间（发前写死；`律 1602` ×0.91）】`Δ ∈ [−84, −33] ms`；**达标所需 Δ ≤ −36.724 ms**。
 *   **不作转绿声称**，只报 Δ 与兑现率。
 * 【值闸门】独立 oracle 双路 + 与底盘**逐位同**（含 `n=m=1e8` 生产档）+ 有牙正控。
 * 行外声明：为不超 102400 字符上限，压缩了四处继承散文（只删叙述，账号/URL/许可全保留）。102315
 * 落款 r135a_。
 * ======================================================================== */
/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/129713: retain lazy odd
 * entry blocks and all original NTT code; inherited authors/licenses below.
 * DIT pass. Its odd-numbered8H blocks feed exclusively into the next pass's
 * first-level Shoup multipliers, so leave their four final butterflies in
 * [0,2P). Even blocks remain canonical. Process one even and one odd block
 * using separate compile-time specializations to keep this choice outside
 * every vector loop. For the4096-point final pass only the even block exists;
 * the16384-point final radix-four pass also first multiplies both odd blocks.
 * twiddle values, transform order, direct-add bounds or final output bounds.
 */
/* 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.
 */
/* 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;
 */
/* 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.
 */
/* 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.
 */
/* 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.
 * 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.
 * 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.
 * keeping the outer scalar caller and passing only pointers across that
 * boundary. Validate with judge-default ISA flags, without global -mavx2.
 * 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
 * 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.
 * 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.
 * 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.
 * 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.
 * 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.
 * 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.
 * 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.
 * 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.
 * 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.
 * 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.
 * quotient products use that same broadcast directly, removing its otherwise
 * redundant byte shift. Retain canonical reduction and all twiddle values.
 */
/* 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.
 * 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.
 */
// ===== 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`）
//     —— 本发的比较基准；本件正文与其同码（差异仅本块纯注释 + 上述单变量）。
// LICENSE：上述公开提交均**无独立许可证声明**、非来自任何外部仓库 ⇒ 按**具名引用**
//     （账号 + 提交地址 + 沿用内容）使用并保留全部上游署名。★本件**未参考任何已发表
//     论文 / 外部开源仓库**。
// ================================================
// ===== 思路（r32gs_ 席 · 同码锚【测量器械】；本席 r135a_ 只删叙述）=====
//   单变量 = 一处保值 token 重拼写（真代码位、非注释）。用途 = 给同批剂量提供紧邻同窗口
//   的前/后标尺。完整预登记判据/同码自证/读数见 notes.md。
// ================================================================
//
//
// ===== 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_ 席；本席 r135a_ 只删叙述压缩）=====
//   单变量 = `shoup8` 的 Montgomery 归约一次**链重结合**（uop 不变、只缩链）。
//   完整读数/相位账/闸门记录见 problems/1002e8/notes.md。
/**
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)
* Eight `_mm_prefetch(_MM_HINT_T0)` per sub-tile aimed at the NEXT 8-row group,
* in `load_tile_avx_d` and in the `mont8` twiddle pass of `inv_natural_stream`
* (strided gather the HW prefetcher cannot track).  No load/store/address or
* arithmetic change: whole-output checksum 7586622636533660264 at n = m = 1e8,
* identical to the base.  Local judge-flag A/B: loader alone -2.3%, loader +
* inverse twiddle -4.0%.  Full prose: problems/1002e8/notes.md.
* 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 (n02q_ seat; 本席 r135a_ 只删叙述) =====
* Two bit-exact changes (chk 9478245449648222414, GI = 0/4/5): two independent
* 8-element groups per `col_dit3_tail` iteration; shortened row-constant live
* ranges in `dft3`/row codelets.  Full prose: problems/1002e8/notes.md.
*/
#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 __attribute__((always_inline)) __m256i vmulhi8h(__m256i a,__m256i b,__m256i bh){
__m256i e = _mm256_mul_epu32(a, b);
__m256i o = _mm256_mul_epu32(_mm256_bsrli_epi128(a, 4), bh);
e = _mm256_bsrli_epi128(e, 4);
return _mm256_blend_epi32(e, o, 0xAA);
}
static inline __attribute__((always_inline)) __m256i shoup8h(__m256i v,__m256i w,__m256i ws,__m256i wsh){
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i t = mullo256(v, w);
__m256i q = vmulhi8h(v, ws, wsh);
__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,bool LazyOutput> static inline __attribute__((always_inline)) void dit_pass8_block(u32 *p,const u32 *tw,const u32 *tws,const u32 *twsh){// [r135a_ SHQH]
const __m256i Pv=_mm256_set1_epi32((int)P);
const u32 *t1=tw+H-1,*s1=tws+H-1;
const u32 *u1=twsh+H-1;
const u32 *t2=tw+2*H-1,*s2=tws+2*H-1;
const u32 *u2=twsh+2*H-1;
const u32 *t3=tw+4*H-1,*s3=tws+4*H-1;
const u32 *u3=twsh+4*H-1;
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,__m256i wsh,bool lazy=false){// [r135a_ SHQH]
__m256i bm=shoup8h(b,w,ws,wsh),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),u_1=ldu256(u1+j);
dit2(v0,v1,w1,ws1,u_1);dit2(v2,v3,w1,ws1,u_1,true);
dit2(v4,v5,w1,ws1,u_1);dit2(v6,v7,w1,ws1,u_1,true);
__m256i w20=ldu256(t2+j),ws20=ldu256(s2+j),u_20=ldu256(u2+j);
__m256i w21=ldu256(t2+j+H),ws21=ldu256(s2+j+H),u_21=ldu256(u2+j+H);
dit2(v0,v2,w20,ws20,u_20);dit2(v4,v6,w20,ws20,u_20,true);
dit2(v1,v3,w21,ws21,u_21);dit2(v5,v7,w21,ws21,u_21,true);
dit2(v0,v4,ldu256(t3+j),ldu256(s3+j),ldu256(u3+j),LazyOutput);
dit2(v1,v5,ldu256(t3+j+H),ldu256(s3+j+H),ldu256(u3+j+H),LazyOutput);
dit2(v2,v6,ldu256(t3+j+2*H),ldu256(s3+j+2*H),ldu256(u3+j+2*H),LazyOutput);
dit2(v3,v7,ldu256(t3+j+3*H),ldu256(s3+j+3*H),ldu256(u3+j+3*H),LazyOutput);
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);
}
}
template<int H> static void dit_pass8(u32*x,int L,const u32*tw,const u32*tws,const u32*twsh){// [r135a_ SHQH]
for(int i=0;i<L;i+=16*H){
dit_pass8_block<H,false>(x+i,tw,tws,twsh);
if(i+8*H<L)dit_pass8_block<H,true>(x+i+8*H,tw,tws,twsh);
}
}
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, const u32 *twsh=nullptr) {// [r135a_ SHQH]
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,twsh);
dit_pass8<64>(x,L,tw,tws,twsh);
dit_pass8<512>(x,L,tw,tws,twsh);
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;
const u32 *us1 = twsh + off, *us2 = us1 + h, *us3 = us2 + h;// [r135a_ SHQH]
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))), US1 = ldu256((const void *)((us1 + j)));
__m256i WV2 = ldu256((const void *)((t2 + j))), WS2 = ldu256((const void *)((ts2 + j))), US2 = ldu256((const void *)((us2 + j)));
__m256i WV3 = ldu256((const void *)((t3 + j))), WS3 = ldu256((const void *)((ts3 + j))), US3 = ldu256((const void *)((us3 + 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 = shoup8h(B, WV1, WS1, US1), Dm = shoup8h(D, WV1, WS1, US1);
__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 = shoup8h(Cp, WV2, WS2, US2), Dn = shoup8h(Dp, WV3, WS3, US3);
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 u32 *twsh) {// [r135a_ SHQH]
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 *u1 = twsh + (len >> 1) - 1;// [r135a_ SHQH]
const u32 *t2 = tw + (len >> 2) - 1, *s2 = tws + (len >> 2) - 1;
const u32 *u2 = twsh + (len >> 2) - 1;
const u32 *t3 = tw + H3 - 1, *s3 = tws + H3 - 1;
const u32 *u3 = twsh + 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 wsh) {// [r135a_ SHQH]
__m256i sm=vaddm256(a,b,Pv),df=_mm256_add_epi32(_mm256_sub_epi32(a,b),Pv);
a=sm;b=shoup8h(df,w,ws,wsh);
};
if constexpr(HalfZero) {
v4=shoup8h(v0,ldu256(t1+j),ldu256(s1+j),ldu256(u1+j));
v5=shoup8h(v1,ldu256(t1+j+H3),ldu256(s1+j+H3),ldu256(u1+j+H3));
v6=shoup8h(v2,ldu256(t1+j+2*H3),ldu256(s1+j+2*H3),ldu256(u1+j+2*H3));
v7=shoup8h(v3,ldu256(t1+j+3*H3),ldu256(s1+j+3*H3),ldu256(u1+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),ldu256(u1+j));
dif2(v1,v5,ldu256(t1+j+H3),ldu256(s1+j+H3),ldu256(u1+j+H3));
dif2(v2,v6,ldu256(t1+j+2*H3),ldu256(s1+j+2*H3),ldu256(u1+j+2*H3));
dif2(v3,v7,ldu256(t1+j+3*H3),ldu256(s1+j+3*H3),ldu256(u1+j+3*H3));
}
__m256i w20=ldu256(t2+j),s20=ldu256(s2+j);
__m256i u20=ldu256(u2+j);
__m256i w21=ldu256(t2+j+H3),s21=ldu256(s2+j+H3);
__m256i u21=ldu256(u2+j+H3);
dif2(v0,v2,w20,s20,u20);dif2(v4,v6,w20,s20,u20);
dif2(v1,v3,w21,s21,u21);dif2(v5,v7,w21,s21,u21);
__m256i w3=ldu256(t3+j),ws3=ldu256(s3+j);
__m256i u3v=ldu256(u3+j);
dif2(v0,v1,w3,ws3,u3v);dif2(v2,v3,w3,ws3,u3v);
dif2(v4,v5,w3,ws3,u3v);dif2(v6,v7,w3,ws3,u3v);
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, const u32 *twsh) {// [r135a_ SHQH]
__m256i Pv = _mm256_set1_epi32((int)P), M1 = _mm256_set1_epi32(-1), Z = _mm256_setzero_si256();
asm("" : "+x"(Pv));   /* [r136a_-F1] 常数寄存器化：禁止再物化 .LC4/.LC6 */
if (L >= 8) {
int len = L;
if (L == 4096) {
dif_pass8<512,HalfZero>(x, L, tw, tws, twsh); dif_pass8<64>(x, L, tw, tws, twsh); dif_pass8<8>(x, L, tw, tws, twsh);
len = 8;
} else if (L == 8192) {
dif_pass8<1024>(x, L, tw, tws, twsh); dif_pass8<128>(x, L, tw, tws, twsh); dif_pass8<16>(x, L, tw, tws, twsh);
len = 16;
}
if (L == 16384) {
dif_pass8<2048>(x,L,tw,tws,twsh); dif_pass8<256>(x,L,tw,tws,twsh); dif_pass8<32>(x,L,tw,tws,twsh);
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;
const u32 *us1 = twsh + h - 1, *us2 = us1 + h2, *us3 = twsh + h2 - 1;// [r135a_ SHQH]
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))), US1 = ldu256((const void *)((us1 + j)));
__m256i WV2 = ldu256((const void *)((t2 + j))), WS2 = ldu256((const void *)((ts2 + j))), US2 = ldu256((const void *)((us2 + j)));
__m256i WV3 = ldu256((const void *)((t3 + j))), WS3 = ldu256((const void *)((ts3 + j))), US3 = ldu256((const void *)((us3 + 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 = shoup8h(d02, WV1, WS1, US1);
d13 = shoup8h(d13, WV2, WS2, US2);
__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 = shoup8h(rb, WV3, WS3, US3);
__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 = shoup8h(rd, WV3, WS3, US3);
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;
const u32 *us = twsh + h - 1;// [r135a_ SHQH]
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 usv = ldu256((const void *)((us + 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 = shoup8h(d, w, ws, usv);
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;
}}

static void build_tab2_h(const u32 *ws, u32 *wsh, int L){
// [r135a_ SHQH] all DIF/DIT load sites sit at element offset == 3 (mod 4), so the
// 128-bit lane boundary is i == 2 (mod 4): wsh[i] = ws[i+1] except at those boundaries.
for (int i = 0; i < L; i++) wsh[i] = (((i + 1) & 3) == 3 || i + 1 >= L) ? 0u : ws[i + 1];
}
#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 u32 tsubFsw_h[M3_], tsubIsw_h[M3_];// [r135a_ SHQH] pre-split Shoup quotients
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_h(tsubFsw, tsubFsw_h, (int)M3_);// [r135a_ SHQH]
build_tab2(tsubIw, tsubIsw, (int)M3_, pw(pw(root, 3), P - 2u));
build_tab2_h(tsubIsw, tsubIsw_h, (int)M3_);// [r135a_ SHQH]
}
#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, tsubFsw_h);
}
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, nullptr, tsubIsw_h);
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,tsubFsw_h);
row_dif_avx(pb,(int)M3_,tsubFw,tsubFsw,tsubFsw_h);
if(arm==0)row_dit_avx<true>(pa,(int)M3_,tsubIw,tsubIsw,pb,tsubIsw_h);
else row_dit_avx<true,true>(pa,(int)M3_,tsubIw,tsubIsw,pb,tsubIsw_h);
}
col_dit3_tail(a);
}
#pragma GCC pop_options
static void row_dit(u32 *x, int L, const u32 *tw, const u32 *tws, const u32 *twsh) {// [r135a_ SHQH]
#if HAVE_AVX2
row_dit_avx(x, L, tw, tws, nullptr, twsh);
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, const u32 *twsh) {// [r135a_ SHQH]
#if HAVE_AVX2
row_dif_avx(x, L, tw, tws, twsh);
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;
u32 *wsRh; u32 *wsC_h;// [r135a_ SHQH]
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, tsubFsw_h);
}
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, e.wsC_h);
}
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, e.wsC_h + 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,e.wsRh);
else row_dif(S+(u64)u*strIn,nB,e.tfR,e.wsR,e.wsRh);
}
}
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, e.wsRh);
}
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, e.wsRh);
}
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, e.wsRh + 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, e.wsRh + 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, u32 *wsRh, u32 *wsC_h) {// [r135a_ SHQH]
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; e.wsRh = wsRh; e.wsC_h = wsC_h;// [r135a_ SHQH]
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));
build_tab2_h(e.wsR, wsRh, e.R); build_tab2_h(e.wsR + e.R, wsRh + e.R, e.R);// [r135a_ SHQH]
if (!e.r5) { build_tab2(tfC, e.wsC, e.C, wC); build_tab2(tiC, e.wsC + e.C, e.C, pw(wC, P - 2u)); build_tab2_h(e.wsC, wsC_h, e.C); build_tab2_h(e.wsC + e.C, wsC_h + e.C, e.C); }// [r135a_ SHQH]
#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 u32 w1Rh[2 * NR_];// [r135a_ SHQH] pre-split Shoup quotients (E1)
static u32 wsRh[NR_];// [r135a_ SHQH] pre-split Shoup quotients (row table)
static u16 b1R[NR_], b1C[1];
static u32 w1R[2 * NR_], w1C[2];
static u32 t2fR[NR2_], t2iR[NR2_], t2fC[NC2_], t2iC[NC2_];
static u32 w2Rh[2 * NR2_], w2Ch[2 * NC2_];// [r135a_ SHQH]
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, w1Rh, w1Rh);
eng_init(E2, LGR2_, (1 << LGC2_), w2, (int)CSX2_, t2fR, t2iR, t2fC, t2iC, b2R, b2C, w2R, w2C, w2Rh, w2Ch);
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;
}

评测结果