/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140880 supplies the full engine, forward-factor boundaries and immediate pair stores; https://duck.ac/submission/140869 introduced the middle-factor boundary; https://duck.ac/submission/140846 introduced pair stores. All author/license notices remain; one ancestor experiment paragraph is omitted.
Idea: apply the factor-lifetime boundary to the ordinary inverse radix-eight middle stage. Complete both first-factor uses before loading the second factor triple, enforced by a zero-instruction compiler barrier with the four computed vectors. Preserve lazy/canonical arithmetic, source reads and store order. The fused64 entry remains unchanged.
Purpose: reduce overlapping inverse twiddle lifetimes; validate scalar DIT/lazy bounds and all200000001 convolution coefficients before official measurement. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/140717 supplies the complete current centered one-prime convolution engine, its resident row subtrees, lazy contracts and all inherited source/license notices retained below. No separate license was declared by that public submission.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/139769: full E8 engine and single-use prime residues, with all inherited author/license notices below.
* Idea: Initialize the auxiliary E2 transform engine, its roots, bit reversal tables and inverse scale only when n+m reaches the N1 alias-correction boundary. E1 setup and all arithmetic remain unchanged. Keep a separate persistent initialization flag so a later call requiring E2 remains valid after earlier ordinary calls.
* Purpose: Remove setup and memory writes for the unused secondary transform on the official workload, while preserving the secondary path for larger degrees. Validate the full200000001 coefficients, naive samples, input hashes and guards before one official measurement.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/139746: accepted forward/inverse codelet fusions and inherited latest phase-separated engine.
* saffah_codex_6a_agg3 https://duck.ac/submission/125996: prime-shaped low-word residue identity, formerly evaluated in the ordinary register-factor path.
* saffah_codex_6a_agg3 https://duck.ac/submission/134208: newer single-use memory-operand Shoup helper. All inherited citations/notices remain below.
* Idea: In shoup8_ptr only, compute t-q*P as (t-q)+(q<<27)-(q<<31) modulo2^32, since P=1+2^31-2^27. Remove its dependent packed-low constant multiply while retaining all three memory-operand products and every load address. Shared-register, broadcast and quarter-root multiplication paths remain unchanged.
* Purpose: Test a shorter quotient-to-residue dependency in the newer single-use consumers. Validate full32-bit multiplier inputs against64-bit modular arithmetic and complete200M convolution outputs; official timing decides the ALU-versus-multiply tradeoff.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/139724: positive forward terminal32 fusion, resident subtrees and inherited full reference chain below.
* Idea: Feed the inverse eight-point entry codelet results directly into the same H=8 radix-eight arithmetic, then perform one final store of the complete64-point block. Preserve each original operation and the alternating canonical/lazy64-block contract. Only4096/8192/16384 rows already using dit_pass8<8> select this path; every other row keeps its original behavior. Product entry and later inverse passes remain unchanged.
* Purpose: Eliminate eight intermediate vector stores and reloads per64 coefficients. Check independent scalar DIT residues/ranges, pointwise entry, all200M convolution coefficients and input/output guards; official timing determines whether increased register pressure offsets the deleted traffic.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/139713: positive2048-coefficient forward subtree scheduling on the latest phase-separated tables. All inherited notices remain below.
* Idea: In forward subtrees whose next remaining stage has32 coefficients, feed its four-vector radix-four results directly into the unchanged final eight-point codelets. Both32-coefficient halves retain their original addresses, twiddles and arithmetic. The common64-coefficient iteration eliminates eight vector stores and the following eight reloads; rows without complete64-coefficient blocks keep the original path. No output normalization or inverse path changes.
* Purpose: Remove a concrete intermediate round trip in the newly cache-resident forward path. Validate scalar DIF residues/ranges and all200M output coefficients/guards; reject if added register pressure outweighs deleted work.
*/
/* References:
* saffah_v41f1017_agg4 https://duck.ac/submission/135178: latest padded/phase-separated twiddle tables and complete radix-three/radix-eight NTT, copied with all original author notices below. No separate source license attached.
* saffah_codex_6a_agg3 https://duck.ac/submission/134208: inherited single-use memory twiddle products and preceding arithmetic, credited by that source.
* Idea: For a16384-coefficient forward row, execute the unchanged top three stages once, then complete each independent2048-coefficient subtree before moving to the next. Each subtree executes the exact old radix-eight256/radix-eight32/radix-four32/terminal codelets with the same table offsets and lazy output flag. This reduces the repeatedly scanned data working set from64KiB to8KiB while preserving the mathematical dependency order. Other row lengths and inverse kernels are unchanged; no extra copying or buffer is introduced.
* Purpose: Evaluate cache blocking of existing arithmetic rather than twiddle-copy/layout tuning. Check full-range scalar DIF oracles for both output contracts, full200M coefficients and input/output guards before official timing.
*/
// ===== 引用区 =====
// [1] duck.ac 用户 saffah_v41f1017_agg4（本账号），提交 #135072
// <https://duck.ac/submission/135072>（1002e8，3.257 s Accepted，现役 best；其正文 = #134317 + fwd 预备 hoist）
// —— 本件正文 = 该提交正文的逐字节副本（其自带的完整引用链原样保留在下方），仅加入下述一处值中性改动。
// [2] duck.ac 题目页 <https://duck.ac/problem/1002e8>（站点明示提交正文公开、无独立许可证）。
// LICENSE：均为本账号/站点公开提交，无第三方许可证；按「账号+URL+沿用内容」具名引用。
// ================================================
/* ===== 引用区（REFERENCES，本账号 saffah_v41f1017_agg4，1002e8）=====
* [1] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#134208**
* <https://duck.ac/submission/134208>（1002e8，3.257191855 s = 当前榜面最快件）。
* 用途：**本文件正文 = 该提交全文的逐字节副本**（其自带的 #134179 / #133417 / 我方 #132971 /
* #129561 及更早全部账号、URL、沿用内容与许可说明原样保留在本块之后，未删改一行）。
* 本件沿用其两处改动：① #134179 单次使用 twiddle 的 VPMULLD/VPMULUDQ **内存操作数** Shoup
* （省去 3 条独立载入与临时寄存器）；② #133417 相邻车道蝶形的**补码合并**（偶道 x+y、
* 奇道 (~x)+y+(P+1)=y−x+P，两条 add + 一条 xor 取代 sum/sub/bias/blend）。
* [2] 我方 #132971 <https://duck.ac/submission/132971>（3.275 s）为其引擎正文底盘（其引用块已注明）。
* [3] duck.ac 题目页 <https://duck.ac/problem/1002e8>；站点明示提交正文公开、无独立许可证。
* 按具名引用（账号 + URL + 沿用内容）使用。
* ===== 思路区 =====
* 目标 1002e8：实时 T=3.257191855 s（#134208，05:06:40）。我方原 #132971=3.275001765 s（早于 T）。
* 本发**不引入任何我方新改动**：仅把对手公开的两把刀吸收为新的自方底盘，恢复与 T 同速
* （预期 ~3.257 s），为下一轮"再快 1%（严支 0.99T=3.2244 s）"的结构搜索建立基线。
* 本地闸门：全尺度 xorshift(n=m=1e8) FNV=9128272874767289449 逐位核对。
* 合规：全部网络经 tools/duck.py；单账号；本件非重掷。
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/134179: complemented
* butterflies and complete exact NTT below. The original offset-loaded
* Shoup products and all inherited authors/licenses are preserved below.
* Our earlier129594 ordinary-load experiment changed whole radix kernels;
* its low multiply still forced register operands. This changes only the
* four distinct single-use twiddles, leaving reused factors in registers.
* Idea: Read single-use w,ws,ws+1 directly as VPMULLD/VPMULUDQ memory
* operands. Eliminate their three standalone loads and temporary registers.
* The same eight data lanes, offset quotient lanes, modulo reduction and
* final stores are retained; no demand address or arithmetic bound changes.
* Purpose: Reduce register pressure and instructions at single-use factors
* in the DIF outer level and DIT final level, including half-zero DIF.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/133417: accepted paired
* correction NTT; all inherited authors, source links and notices below.
* Its adjacent-lane butterfly and merged modular correction are retained.
* Idea: Form alternating sum/difference lanes with one complemented input:
* even x+y, odd (~x)+y+(P+1)=y-x+P modulo2^32. Two adds and one XOR
* replace separate sum, subtraction, bias-add and blend. Preserve the exact
* existing canonical/lazy correction at each forward/inverse codelet.
* Purpose: Remove one SIMD operation per adjacent-lane butterfly.
*/
/* References:
* saffah_v41f1017_agg4 https://duck.ac/submission/132971: copied the complete
* radix-three NTT engine and paired correction loop, with all inherited
* references and notices below. No additional external source is used.
* saffah_codex_6a_agg3 https://duck.ac/submission/129561: earlier single-stream
* row cursor experiment; it did not improve that old two-pass correction.
* Idea: The new paired correction has two padded-A streams, the second with
* an unaligned row origin. Split at the earlier of their next row boundaries.
* Inside each long run, advance two contiguous pointers; only the one vector
* crossing the second stream's padded boundary uses scalar gathering. This
* removes both repeated division-based addresses AND the per-vector boundary
* branch newly introduced by paired correction. Prefix, rearrangement and
* all store behavior are preserved. Purpose: test this different fused-loop
* setting; no repeated submission of the earlier single-stream experiment.
*/
/* ===== 引用区 =====
* [1] duck.ac 用户 saffah_cc_v41_agg1，提交 #130739 <https://duck.ac/submission/130739>
* （1002e8，3.313 s，Accepted）—— 本件引擎正文的底盘（两阶段 radix-2/3/8 DIF/DIT、
* Shoup/Montgomery 归约、corner 流、进位与输出合成）。
* [2] 本账号 saffah_v41f1017_agg4，提交 #132721 <https://duck.ac/submission/132721>
* （3.302 s，Accepted，本账号当前最好）—— 沿用其三处改动：两段 cmbv 融合、psum4、
* corner_tile 的 mont8x。以上均为站点公开提交（FAQ 明示公开），无独立许可证，
* 按具名引用（账号 + URL + 沿用内容）。本件未参考任何论文/外部仓库。
* ===== 思路区 =====
* 目标 1002e8：n=m=1e8，T=3.313 s，严支阈值 3.279871 s，我方 3.302 s，缺口 0.675%。
* 判题机分相（本席自旋探针实测，相对底盘）：fwd load 209、fwd rowdif 622、fwd corner_tile 217、
* fwd head 239、inv twcons 155、inv rowdit 354、inv writeback 64、rowconv 1349、cmbv ~105 ms。
* 本发目的与改动见下（单变量，试验性；输出逐位不变）。
*/
/* IDEA: 基座 = 当前最好件 3.299s(sid 132901, ct 预取 +8) + 本发单变量：load_block 的 nsrc 边界行组不再走 64 行标量回退，改为 8x8 瓦片 + 行级判定（整行在界内→向量载入；整行越界→零；唯一跨界行→8 元素标量拼装）。依据：fwd load 段的最后一个 64 行组(r0=1984)跨越 nsrc=100000001，原标量路径对每 stripe/臂/两次 fwd 都执行；本机同构微基准（含边界）3.06M ticks → 修正后 2.19M(-28%)，理想无边界为 1.94M。判题机 load 段 209ms，按 32 组摊薄预期省 14-19ms。psum 累加顺序变化但总和不变；逐位等价(FNV 同底盘)。 */
/* ===== 引用区 =====
* 本件正文 = duck.ac 公开提交 #130739（账号 saffah_cc_v41_agg1，
* <https://duck.ac/submission/130739>，3.313 s，Accepted）的逐字节副本 + 本席单变量（见思路区）。
* 沿用内容 = 完整两阶段多项式乘法引擎（radix-2/3/8 DIF/DIT 内核、Shoup/Montgomery 归约、
* corner 流、进位与输出合成）。上游引用链（全部 duck.ac 公开提交，按站点规则具名引用：
* 账号 + 提交号 + 完整 URL + 沿用内容；下列原文为原链逐代保留的引用/许可文字，仅删叙述行）：
*
* ===== REFERENCES =====
* [1] duck.ac 用户 **saffah_codex_6a_agg3**（对手），提交 **#129718**
* <https://duck.ac/submission/129718>（判题机 3349.773463 ms = 本题 T）。用途 = 本文件正文的**底盘**。
* [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#130671**
* <https://duck.ac/submission/130671>（3349.856085 ms = 我方现役最好）。用途 = 本件的**直接母本**。
* （该件自身 = 本账号 **#130554** <https://duck.ac/submission/130554> 的单变量回退件。）
* [3] 题目页 <https://duck.ac/problem/1002e8>（站点 FAQ：提交正文公开可见）。
* LICENSE：以上均为站点公开提交、按具名引用（账号 + 提交 URL + 沿用内容）；本件未参考任何论文/外部仓库。
* 站点 = `dif_pass8`+`dit_pass8_block`（同 #130711），**外加** `row_dit_avx<Product>` 的 `initial_load` 里 `mont8` 的两个输入移位改偏移载入（`mont8x`）——该站点是 `col_convolve3` 融合点乘的最大消费者。
* 目的（试验性提交必写）= 判题机定价：`#130711` = `−32.471294 ms`，只差 `1.108063 ms` 转绿；`#130719`（5 个 mont8 站点一起换）实测为 **回归 +10.3 ms** ⇒ 本发只换**最热的一个**站点，检验"站点越多越差"是剂量效应还是个别站点反向。本题**红** ⇒ 免费。
* 现读（发前 `tools/exact.py`）：mine = 3317.384791 ms（#130711）· T = 3349.773463 ms（#129718）· **严支**
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129713: retain lazy odd
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129697: retain the complete
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129683: retain the tile-wide
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129663: retain lazy forward
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129659: retain broadcast
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129512: retain the complete
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129498: shared forward
* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/129474: accepted prepared
* References:
* saffah_cc_v41_agg1 https://duck.ac/submission/129436: retain its fixed
* ===== REFERENCES（依目标原文 §"提交的代码中…列出引用的部分"，写在提交代码【开头】）=====
* [1] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#128650** <https://duck.ac/submission/128650>
* （#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 +
* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128584: retain the accepted
* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128547: use the accepted
* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128341: retain complete
* References:
* saffah_codex_6a_agg3, https://duck.ac/submission/128318: retain our
* saffah_cc_v41_agg1, https://duck.ac/submission/127827: retain the exact
* saffah_codex_6a_agg3, https://duck.ac/submission/128249 (1002i): transfer
* References:
* saffah_cc_v41_agg1, https://duck.ac/submission/127827: retain the complete
* saffah_codex_6a_agg3, https://duck.ac/submission/128304: adapt our output
* ===== REFERENCES =====
* [1] duck.ac 用户 **saffah_cc_v41_agg1**，提交 **#127236** <https://duck.ac/submission/127236>
* **逐字节副本**；其自带的 REFERENCES 区与思路段紧随本块之下、**原样保留、未作删改**。
* [2] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#127447**
* <https://duck.ac/submission/127447>（3553.905713 ms = 本账号本题现役最快件 `mine`）
* [3] duck.ac 题目页 <https://duck.ac/problem/1002e8>（提交正文按站点规则公开可见）。
* LICENSE：上述公开提交均**无独立许可证声明**、非来自任何外部仓库 ⇒ 按**具名引用**
* ===== REFERENCES（依目标原文 §"提交的代码中…列出引用的部分"）=====
* [1] duck.ac 用户 **saffah_cc_v41_agg1**（本账号），提交 **#127236**
* <https://duck.ac/submission/127236>（3555.554282 ms = 本题现役最快件 / `mine`）。
* 该件自带的上游引用链（其对 `saffah_codex_6a_agg3` 各代与其自身历史件的引用与许可说明）
* [2] duck.ac 用户 **saffah_codex_6a_agg3**，提交 **#126201** <https://duck.ac/submission/126201>
* LICENSE：二者均为 duck.ac 公开提交（题面明示"你提交的代码将会被公开，所有人都可见"）；
* 无第三方许可证声明、非来自任何外部仓库 ⇒ 按具名引用（账号 + 提交地址 + 沿用内容）使用，
* 遵守其无附加限制的许可。**本件未参考任何已发表论文 / 外部开源仓库。**
* References:
* - saffah_codex_6a_agg3, https://duck.ac/submission/126201 (3579.423873 ms,
* - saffah_cc_v41_agg1, https://duck.ac/submission/127221 (3563.511796 ms,
* (the fused inverse twiddle consumer, submission #125908 of the other
* account, and this account's #123572).
* References:
* - saffah_codex_6a_agg3, https://duck.ac/submission/125996 : use our accepted
* - saffah_cc_v41_agg1, https://duck.ac/submission/123572 : retain the radix-eight
* References:
* - saffah_codex_6a_agg3, https://duck.ac/submission/125941 : retained our
* - saffah_cc_v41_agg1, https://duck.ac/submission/123572 : inherited paired
* References:
* - saffah_codex_6a_agg3, https://duck.ac/submission/125908 : retained our
* - saffah_cc_v41_agg1, https://duck.ac/submission/123953 : adapted its
* - saffah_cc_v41_agg1, https://duck.ac/submission/123572 : inherited all
* References:
* - saffah_cc_v41_agg1, https://duck.ac/submission/123572 : retained its
* - saffah_codex_6a_agg3, https://duck.ac/submission/123490 : inherited our
* ===== REFERENCES =====
* [1] duck.ac user **saffah_codex_6a_agg3**, submission **#123490**
* <https://duck.ac/submission/123490> (3676.709519 ms, Accepted) -- the body
* [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
* duck.ac user **saffah_codex_6a_agg3**:
* #122149 <https://duck.ac/submission/122149> 3745.686476 -- lazy forward
* #122093 <https://duck.ac/submission/122093> -- lazy variance at the
* #122042 <https://duck.ac/submission/122042> 4014.333747 -- exact-parity
* #121791 <https://duck.ac/submission/121791> 4104.879166 -- streamed radix-three
* #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
* #121337 <https://duck.ac/submission/121337> 4254.810910 -- 64-column slab
* #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
* #122215 <https://duck.ac/submission/122215> 3719.856168 -- lazy Shoup operands;
* #122124 <https://duck.ac/submission/122124> 3809.768898 -- paired live-lane row
* #121936 <https://duck.ac/submission/121936> 4038.313189 -- fused forward
* #121886 <https://duck.ac/submission/121886> 4091.567105, #121631
* <https://duck.ac/submission/121631> 4149.798632 -- fused four-sweep row
* #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
* 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,
* - saffah_codex_6a_agg3, https://duck.ac/submission/122149 : inherited our
* [n02q_] ===== REFERENCES =====
* [1] This account **saffah_cc_v41_agg1**, submission **#122215**
* <https://duck.ac/submission/122215> (3719.856168 ms, Accepted) -- the body of
* 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
* [2] duck.ac user **saffah_codex_6a_agg3**, submissions #122149 / #122042 / #122093
* and this account's #122124 / #121936: the row codelets and memory schedule
*
* LICENSE：以上均为 duck.ac 公开提交（站点 FAQ 明示提交正文公开可见），无独立许可证声明、
* 非来自任何外部仓库 ⇒ 按具名引用使用并保留署名。本件未参考任何已发表论文/外部开源仓库。
* ===== 思路区 =====
* 底盘 #130739 的历史改动（继承自本账号前作 #131591/#132721，逐位等值，本机全尺度 FNV 一致）：
* 1) 融合 cmbv：m=n 时第二段输出前缀 SUM2(k+n+1)=F(n)−F(k)，两段 cmbv 合为单遍、a,b 只读一遍；
* 伙伴流 vperm2i128+vpalignr 重排成 32B 对齐后 NT 存（消除 400MB RFO）。
* 2) F(n)=4(PA+PB)−16(n+1)，PA/PB 在前向 load 内核顺带累加（4 条独立链）。
* 3) corner_tile_avx_to 的 mont8 改 mont8x（偏移载入替代 p5 移位）。
* 本席本发 = 单变量试验（见下），目的 = 判题机定价；等值性以双路 oracle + 全尺度 FNV 闸门。
* 本发单变量说明：
* （本地标定件：与 #132721 正文逐字节同码，仅头部压缩；用于取得压缩底盘的本地 FNV 基准。）
* ================================================ */
#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
#define PADC_ 16
#define PADC2_ 16
#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_)
#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 + 8];
static inline u32 mm(u32 a, u32 b) {
u64 z = (u64)a * b;
u64 q = (u64)(((u128)z * IM) >> 64);
u64 r = z - q * (u64)P;
return (u32)(r >= P ? r - P : r);
}
static inline u32 shoup_c(u32 w) { return (u32)(((u64)w << 32) / (u64)P); }
static u32 pw(u32 a, u32 e) {
u32 r = 1;
while (e) { if (e & 1) r = (u32)((u64)r * a % P); a = (u32)((u64)a * a % P); e >>= 1; }
return r;
}
static void build_tab(u32 *tw, int L, u32 w) {
int off = 0;
for (int len = 2; len <= L; len <<= 1) {
u32 wl = pw(w, (u32)(L / len));
u32 c = 1;
int h = len >> 1;
for (int j = 0; j < h; j++) { tw[off + j] = c; c = mm(c, wl); }
off += h;
}
}
static void build_br(u16 *br, int L) {
int bits = 0; while ((1 << bits) < L) bits++;
for (int i = 0; i < L; i++) {
int r = 0;
for (int b = 0; b < bits; b++) if (i & (1 << b)) r |= 1 << (bits - 1 - b);
br[i] = (u16)r;
}
}
#ifdef NO_AVX2
#define HAVE_AVX2 0
#else
#define HAVE_AVX2 1
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __m256i vmulhi8(__m256i a, __m256i b) {
__m256i e = _mm256_mul_epu32(a, b);
__m256i o = _mm256_mul_epu32(_mm256_bsrli_epi128(a, 4), _mm256_bsrli_epi128(b, 4));
e = _mm256_bsrli_epi128(e, 4);
return _mm256_blend_epi32(e, o, 0xAA);
}
static inline __attribute__((always_inline)) __m256i mullo256(__m256i a, __m256i b) {
__m256i r; __asm__("vpmulld %2, %1, %0" : "=x"(r) : "x"(a), "x"(b)); return r;
}
static inline __m256i shoup8(__m256i v, __m256i w, __m256i ws) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i t = mullo256(v, w);
__m256i q = vmulhi8(v, ws);
__m256i qp = mullo256(q, Pv);
__m256i r = _mm256_sub_epi32(t, qp);
return _mm256_min_epu32(r, _mm256_sub_epi32(r, Pv));
}
static inline __m256i vmulhi8x(__m256i a, __m256i b, __m256i bo) {
__m256i e = _mm256_mul_epu32(a, b);
__m256i o = _mm256_mul_epu32(_mm256_bsrli_epi128(a, 4), bo);
return _mm256_blend_epi32(_mm256_bsrli_epi128(e, 4), o, 0xAA);
}
static inline __m256i shoup8x(__m256i v, __m256i w, __m256i ws, __m256i wo) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i t = mullo256(v, w);
__m256i q = vmulhi8x(v, ws, wo);
__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 shoup8_ptr(__m256i v,const u32*w,const u32*ws){
const __m256i pv=_mm256_set1_epi32((int)P);
__m256i t,qe,qo;const __m256i vo=_mm256_bsrli_epi128(v,4);
asm("vpmulld %2,%1,%0":"=x"(t):"x"(v),"m"(*(const __m256i*)w));
asm("vpmuludq %2,%1,%0":"=x"(qe):"x"(v),"m"(*(const __m256i*)ws));
asm("vpmuludq %2,%1,%0":"=x"(qo):"x"(vo),"m"(*(const __m256i*)(ws+1)));
const __m256i q=_mm256_blend_epi32(_mm256_bsrli_epi128(qe,4),qo,0xaa);
const __m256i r=_mm256_add_epi32(_mm256_sub_epi32(t,q),_mm256_sub_epi32(_mm256_slli_epi32(q,27),_mm256_slli_epi32(q,31)));
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 __m256i Pv=_mm256_set1_epi32((int)P);
const u32 *t1=tw+H-1,*s1=tws+H-1;
const u32 *t2=tw+2*H-1,*s2=tws+2*H-1;
const u32 *t3=tw+4*H-1,*s3=tws+4*H-1;
for(int 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 wo,bool lazy=false){
__m256i bm=shoup8x(b,w,ws,wo),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);}
};
auto dit2p=[&](__m256i&a,__m256i&b,const u32*w,const u32*ws,bool lazy){
__m256i bm=shoup8_ptr(b,w,ws),old=a;
if(lazy){a=_mm256_add_epi32(old,bm);b=_mm256_add_epi32(_mm256_sub_epi32(old,bm),Pv);}
else {a=vaddm256(old,bm,Pv);b=vsubm256(old,bm,Pv);}
};
__m256i w1=ldu256(t1+j),ws1=ldu256(s1+j),ws1o=ldu256(s1+j+1);
dit2(v0,v1,w1,ws1,ws1o);dit2(v2,v3,w1,ws1,ws1o,true);
dit2(v4,v5,w1,ws1,ws1o);dit2(v6,v7,w1,ws1,ws1o,true);
{
const __m256i w20=ldu256(t2+j),ws20=ldu256(s2+j),ws20o=ldu256(s2+j+1);
dit2(v0,v2,w20,ws20,ws20o);dit2(v4,v6,w20,ws20,ws20o,true);
}
asm volatile(""::"x"(v0),"x"(v2),"x"(v4),"x"(v6):"memory");
{
const __m256i w21=ldu256(t2+j+H),ws21=ldu256(s2+j+H),ws21o=ldu256(s2+j+H+1);
dit2(v1,v3,w21,ws21,ws21o);dit2(v5,v7,w21,ws21,ws21o,true);
}
dit2p(v0,v4,t3+j,s3+j,LazyOutput);
dit2p(v1,v5,t3+j+H,s3+j+H,LazyOutput);
dit2p(v2,v6,t3+j+2*H,s3+j+2*H,LazyOutput);
dit2p(v3,v7,t3+j+3*H,s3+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){
for(int i=0;i<L;i+=16*H){
dit_pass8_block<H,false>(x+i,tw,tws);
if(i+8*H<L)dit_pass8_block<H,true>(x+i+8*H,tw,tws);
}
}
template<bool LazyOutput> static inline __attribute__((always_inline)) void dit_entry64_registers(
__m256i&v0,__m256i&v1,__m256i&v2,__m256i&v3,__m256i&v4,__m256i&v5,__m256i&v6,__m256i&v7,const u32*tw,const u32*tws){
constexpr int H=8,j=0;
const __m256i Pv=_mm256_set1_epi32((int)P);
const u32*t1=tw+7,*s1=tws+7,*t2=tw+15,*s2=tws+15,*t3=tw+31,*s3=tws+31;
auto dit2=[&](__m256i &a,__m256i &b,__m256i w,__m256i ws,__m256i wo,bool lazy=false){
__m256i bm=shoup8x(b,w,ws,wo),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);}
};
auto dit2p=[&](__m256i&a,__m256i&b,const u32*w,const u32*ws,bool lazy){
__m256i bm=shoup8_ptr(b,w,ws),old=a;
if(lazy){a=_mm256_add_epi32(old,bm);b=_mm256_add_epi32(_mm256_sub_epi32(old,bm),Pv);}
else {a=vaddm256(old,bm,Pv);b=vsubm256(old,bm,Pv);}
};
__m256i w1=ldu256(t1+j),ws1=ldu256(s1+j),ws1o=ldu256(s1+j+1);
dit2(v0,v1,w1,ws1,ws1o);dit2(v2,v3,w1,ws1,ws1o,true);
dit2(v4,v5,w1,ws1,ws1o);dit2(v6,v7,w1,ws1,ws1o,true);
__m256i w20=ldu256(t2+j),ws20=ldu256(s2+j),ws20o=ldu256(s2+j+1);
__m256i w21=ldu256(t2+j+H),ws21=ldu256(s2+j+H),ws21o=ldu256(s2+j+H+1);
dit2(v0,v2,w20,ws20,ws20o);dit2(v4,v6,w20,ws20,ws20o,true);
dit2(v1,v3,w21,ws21,ws21o);dit2(v5,v7,w21,ws21,ws21o,true);
dit2p(v0,v4,t3+j,s3+j,LazyOutput);
dit2p(v1,v5,t3+j+H,s3+j+H,LazyOutput);
dit2p(v2,v6,t3+j+2*H,s3+j+2*H,LazyOutput);
dit2p(v3,v7,t3+j+3*H,s3+j+3*H,LazyOutput);
}
static inline __m256i shoup_of(__m256i w);
static inline __m256i mont8(__m256i a, __m256i b);
static inline __m256i mont8x(__m256i a, __m256i ax, __m256i b, __m256i bx);
template<bool Product=false,bool LazyTail=false>
static void row_dit_avx(u32 *x, int L, const u32 *tw, const u32 *tws, const u32 *other=nullptr) {
auto initial_load=[&](int offset) {
__m256i v=ldu256(x+offset);
if constexpr(Product){__m256i y=ldu256(other+offset);
v=mont8x(v,ldu256(x+offset+1),y,ldu256(other+offset+1));}
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);
vr0=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr0,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr0 = _mm256_min_epu32(vr0,_mm256_sub_epi32(vr0,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr1, 0xB1);
vr1=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr1,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr1 = _mm256_min_epu32(vr1,_mm256_sub_epi32(vr1,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr2, 0xB1);
vr2=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr2,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr2 = _mm256_min_epu32(vr2,_mm256_sub_epi32(vr2,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr3, 0xB1);
vr3=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr3,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr3 = _mm256_min_epu32(vr3,_mm256_sub_epi32(vr3,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr4, 0xB1);
vr4=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr4,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr4 = _mm256_min_epu32(vr4,_mm256_sub_epi32(vr4,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr5, 0xB1);
vr5=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr5,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr5 = _mm256_min_epu32(vr5,_mm256_sub_epi32(vr5,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr6, 0xB1);
vr6=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr6,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
vr6 = _mm256_min_epu32(vr6,_mm256_sub_epi32(vr6,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr7, 0xB1);
vr7=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr7,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
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);
}
}
if(L==4096||L==8192||L==16384){
if(i4&64)dit_entry64_registers<true>(vr0,vr1,vr2,vr3,vr4,vr5,vr6,vr7,tw,tws);
else dit_entry64_registers<false>(vr0,vr1,vr2,vr3,vr4,vr5,vr6,vr7,tw,tws);
}
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);
v=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(v,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
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<64>(x,L,tw,tws);
dit_pass8<512>(x,L,tw,tws);
len=8192;off=4095;
}
for (; (len << 1) <= L; len <<= 2) {
int h = len >> 1;
const u32 *t1 = tw + off, *ts1 = tws + off;
const u32 *t2 = t1 + h, *ts2 = ts1 + h;
const u32 *t3 = t2 + h, *ts3 = ts2 + h;
off += 3 * h;
for (int i = 0; i < L; i += (len << 1)) {
u32 *p = x + i;
for (int j = 0; j < h; j += 8) {
__m256i WV1 = ldu256((const void *)((t1 + j))), WS1 = ldu256((const void *)((ts1 + j)));
__m256i WV2 = ldu256((const void *)((t2 + j))), WS2 = ldu256((const void *)((ts2 + j)));
__m256i WV3 = ldu256((const void *)((t3 + j))), WS3 = ldu256((const void *)((ts3 + j)));
__m256i A = ldu256((const void *)(p + j));
__m256i B = ldu256((const void *)(p + j + h));
__m256i C = ldu256((const void *)(p + j + 2 * h));
__m256i D = ldu256((const void *)(p + j + 3 * h));
__m256i Bm = shoup8(B, WV1, WS1), Dm = shoup8(D, WV1, WS1);
__m256i Ap = vaddm256(A, Bm, Pv), Bp = vsubm256(A, Bm, Pv);
__m256i Cp = _mm256_add_epi32(C, Dm), Dp = _mm256_add_epi32(_mm256_sub_epi32(C, Dm), Pv);
__m256i Cm = shoup8(Cp, WV2, WS2), Dn = shoup8(Dp, WV3, WS3);
if(LazyTail && (len<<1)==L){
stu256(p+j,_mm256_add_epi32(Ap,Cm));
stu256(p+j+2*h,_mm256_add_epi32(_mm256_sub_epi32(Ap,Cm),Pv));
stu256(p+j+h,_mm256_add_epi32(Bp,Dn));
stu256(p+j+3*h,_mm256_add_epi32(_mm256_sub_epi32(Bp,Dn),Pv));
}else{
stu256((void *)(p + j), vaddm256(Ap, Cm, Pv));
stu256((void *)(p + j + 2 * h), vsubm256(Ap, Cm, Pv));
stu256((void *)(p + j + h), vaddm256(Bp, Dn, Pv));
stu256((void *)(p + j + 3 * h), vsubm256(Bp, Dn, Pv));
}
}
}
}
if (len <= L) {
int h = len >> 1;
const u32 *t = tw + off, *ts = tws + off; off += h;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j += 8) {
__m256i u = ldu256((const void *)(p + j));
__m256i v = shoup8(ldu256((const void *)(q + j)),
ldu256((const void *)(t + j)), ldu256((const void *)(ts + j)));
stu256((void *)(p + j), vaddm256(u, v, Pv));
stu256((void *)(q + j), vsubm256(u, v, Pv));
}
}
}
} else {
for (int len = 2; len <= L; len <<= 1) {
int h = len >> 1;
const u32 *t = tw + off; off += h;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = mm(q[j], t[j]);
u32 sa = u + v; if (sa >= P) sa -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = sa; q[j] = d;
}
}
}
}
}
template<int H3, bool HalfZero=false> static void dif_pass8(u32 *x, int L, const u32 *tw, const u32 *tws) {
const __m256i Pv = _mm256_set1_epi32((int)P);
const int len = H3 << 3;
const u32 *t1 = tw + (len >> 1) - 1, *s1 = tws + (len >> 1) - 1;
const u32 *t2 = tw + (len >> 2) - 1, *s2 = tws + (len >> 2) - 1;
const u32 *t3 = tw + H3 - 1, *s3 = tws + H3 - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i;
for (int j = 0; j < H3; j += 8) {
__m256i v0=ldu256(p+j),v1=ldu256(p+j+H3),
v2=ldu256(p+j+2*H3),v3=ldu256(p+j+3*H3);
__m256i v4,v5,v6,v7;
auto dif2=[&](__m256i &a,__m256i &b,__m256i w,__m256i ws,__m256i wo) {
__m256i sm=vaddm256(a,b,Pv),df=_mm256_add_epi32(_mm256_sub_epi32(a,b),Pv);
a=sm;b=shoup8x(df,w,ws,wo);
};
auto dif2p=[&](__m256i&a,__m256i&b,const u32*w,const u32*ws){
__m256i sm=vaddm256(a,b,Pv),df=_mm256_add_epi32(_mm256_sub_epi32(a,b),Pv);
a=sm;b=shoup8_ptr(df,w,ws);
};
if constexpr(HalfZero) {
v4=shoup8_ptr(v0,t1+j,s1+j);
v5=shoup8_ptr(v1,t1+j+H3,s1+j+H3);
v6=shoup8_ptr(v2,t1+j+2*H3,s1+j+2*H3);
v7=shoup8_ptr(v3,t1+j+3*H3,s1+j+3*H3);
} else {
v4=ldu256(p+j+4*H3);v5=ldu256(p+j+5*H3);
v6=ldu256(p+j+6*H3);v7=ldu256(p+j+7*H3);
dif2p(v0,v4,t1+j,s1+j);
dif2p(v1,v5,t1+j+H3,s1+j+H3);
dif2p(v2,v6,t1+j+2*H3,s1+j+2*H3);
dif2p(v3,v7,t1+j+3*H3,s1+j+3*H3);
}
{
const __m256i w20=ldu256(t2+j),s20=ldu256(s2+j),s20o=ldu256(s2+j+1);
dif2(v0,v2,w20,s20,s20o);dif2(v4,v6,w20,s20,s20o);
}
asm volatile(""::"x"(v0),"x"(v2),"x"(v4),"x"(v6):"memory");
{
const __m256i w21=ldu256(t2+j+H3),s21=ldu256(s2+j+H3),s21o=ldu256(s2+j+H3+1);
dif2(v1,v3,w21,s21,s21o);dif2(v5,v7,w21,s21,s21o);
}
asm volatile(""::"x"(v1),"x"(v3),"x"(v5),"x"(v7):"memory");
__m256i w3=ldu256(t3+j),ws3=ldu256(s3+j),ws3o=ldu256(s3+j+1);
dif2(v0,v1,w3,ws3,ws3o);
stu256(p+j,v0);stu256(p+j+H3,v1);
dif2(v2,v3,w3,ws3,ws3o);
stu256(p+j+2*H3,v2);stu256(p+j+3*H3,v3);
dif2(v4,v5,w3,ws3,ws3o);
stu256(p+j+4*H3,v4);stu256(p+j+5*H3,v5);
dif2(v6,v7,w3,ws3,ws3o);
stu256(p+j+6*H3,v6);stu256(p+j+7*H3,v7);
}
}
}
template<bool HalfZero=false, bool LazyOut=false, bool ResidentInner=false>
static void row_dif_avx(u32 *x, int L, const u32 *tw, const u32 *tws) {
__m256i Pv = _mm256_set1_epi32((int)P), M1 = _mm256_set1_epi32(-1), Z = _mm256_setzero_si256();
if (L >= 8) {
int len = L;
if constexpr(ResidentInner){
dif_pass8<256>(x,L,tw,tws);dif_pass8<32>(x,L,tw,tws);len=32;
}else{
if(L==16384){
dif_pass8<2048>(x,L,tw,tws);
for(int block=0;block<L;block+=2048)row_dif_avx<false,LazyOut,true>(x+block,2048,tw,tws);
return;
}
if (L == 4096) {
dif_pass8<512,HalfZero>(x, L, tw, tws); dif_pass8<64>(x, L, tw, tws); dif_pass8<8>(x, L, tw, tws);
len = 8;
} else if (L == 8192) {
dif_pass8<1024>(x, L, tw, tws); dif_pass8<128>(x, L, tw, tws); dif_pass8<16>(x, L, tw, tws);
len = 16;
}
if (L == 16384) {
dif_pass8<2048>(x,L,tw,tws); dif_pass8<256>(x,L,tw,tws); dif_pass8<32>(x,L,tw,tws);
len=32;
}
}
const bool Fuse32=(len==32 && (L&63)==0);
if (!Fuse32 && L != 4096 && L != 8192)
for (; len >= 32; len >>= 2) {
int h = len >> 1, h2 = len >> 2;
const u32 *t1 = tw + h - 1, *ts1 = tws + h - 1;
const u32 *t2 = t1 + h2, *ts2 = ts1 + h2;
const u32 *t3 = tw + h2 - 1, *ts3 = tws + h2 - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i;
for (int j = 0; j < h2; j += 8) {
__m256i WV1 = ldu256((const void *)((t1 + j))), WS1 = ldu256((const void *)((ts1 + j)));
__m256i WV2 = ldu256((const void *)((t2 + j))), WS2 = ldu256((const void *)((ts2 + j)));
__m256i WV3 = ldu256((const void *)((t3 + j))), WS3 = ldu256((const void *)((ts3 + j)));
__m256i v0 = ldu256((const void *)((p + j)));
__m256i v1 = ldu256((const void *)((p + j + h2)));
__m256i v2 = ldu256((const void *)((p + j + h)));
__m256i v3 = ldu256((const void *)((p + j + h + h2)));
__m256i s02 = _mm256_min_epu32(_mm256_add_epi32(v0, v2), _mm256_sub_epi32(_mm256_add_epi32(v0, v2), Pv));
__m256i d02 = _mm256_sub_epi32(v0, v2);
__m256i s13 = _mm256_min_epu32(_mm256_add_epi32(v1, v3), _mm256_sub_epi32(_mm256_add_epi32(v1, v3), Pv));
__m256i d13 = _mm256_sub_epi32(v1, v3);
d02 = _mm256_add_epi32(d02, Pv);
d13 = _mm256_add_epi32(d13, Pv);
d02 = shoup8(d02, WV1, WS1);
d13 = shoup8(d13, WV2, WS2);
__m256i ra = _mm256_min_epu32(_mm256_add_epi32(s02, s13), _mm256_sub_epi32(_mm256_add_epi32(s02, s13), Pv));
__m256i rb = _mm256_add_epi32(_mm256_sub_epi32(s02, s13), Pv);
rb = shoup8(rb, WV3, WS3);
__m256i rc = _mm256_min_epu32(_mm256_add_epi32(d02, d13), _mm256_sub_epi32(_mm256_add_epi32(d02, d13), Pv));
__m256i rd = _mm256_add_epi32(_mm256_sub_epi32(d02, d13), Pv);
rd = shoup8(rd, WV3, WS3);
stu256((void *)((p + j)), ra);
stu256((void *)((p + j + h2)), rb);
stu256((void *)((p + j + h)), rc);
stu256((void *)((p + j + h + h2)), rd);
}
}
}
if(Fuse32)len=8;
if (len >= 16) {
int h = len >> 1;
const u32 *t = tw + h - 1, *ts = tws + h - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j += 8) {
__m256i u = ldu256((const void *)((p + j)));
__m256i v = ldu256((const void *)((q + j)));
__m256i w = ldu256((const void *)((t + j)));
__m256i ws = ldu256((const void *)((ts + j)));
__m256i sa = _mm256_min_epu32(_mm256_add_epi32(u, v), _mm256_sub_epi32(_mm256_add_epi32(u, v), Pv));
__m256i d = _mm256_add_epi32(_mm256_sub_epi32(u, v), Pv);
d = shoup8(d, w, ws);
stu256((void *)((p + j)), sa);
stu256((void *)((q + j)), d);
}
}
}
const __m256i V8 = _mm256_setr_epi32(1, 1, 1, 1, (int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6]);
const __m256i V8s = _mm256_setr_epi32(2, 2, 2, 2, (int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6]);
const __m256i V8N = _mm256_setr_epi32((int)tw[3], (int)tw[4], (int)tw[5], (int)tw[6],
(int)(P - tw[3]), (int)(P - tw[4]), (int)(P - tw[5]), (int)(P - tw[6]));
const __m256i V8Ns = _mm256_setr_epi32((int)tws[3], (int)tws[4], (int)tws[5], (int)tws[6],
(int)shoup_c(P - tw[3]), (int)shoup_c(P - tw[4]), (int)shoup_c(P - tw[5]), (int)shoup_c(P - tw[6]));
const __m256i V4 = _mm256_setr_epi32(1, 1, (int)tw[1], (int)tw[2], 1, 1, (int)tw[1], (int)tw[2]);
const __m256i V4s = _mm256_setr_epi32(2, 2, (int)tws[1], (int)tws[2], 2, 2, (int)tws[1], (int)tws[2]);
const __m256i V4P=_mm256_set1_epi64x((long long)tw[2]);
const __m256i V4Ps=_mm256_set1_epi64x((long long)tws[2]);
const __m256i IDX8 = _mm256_setr_epi32(4, 5, 6, 7, 0, 1, 2, 3);
auto finish32=[&](__m256i&v0,__m256i&v1,__m256i&v2,__m256i&v3){
__m256i s02=vaddm256(v0,v2,Pv),s13=vaddm256(v1,v3,Pv);
__m256i d02=shoup8_ptr(_mm256_add_epi32(_mm256_sub_epi32(v0,v2),Pv),tw+15,tws+15);
__m256i d13=shoup8_ptr(_mm256_add_epi32(_mm256_sub_epi32(v1,v3),Pv),tw+23,tws+23);
__m256i ra=vaddm256(s02,s13,Pv),rc=vaddm256(d02,d13,Pv);
__m256i rb=shoup8_ptr(_mm256_add_epi32(_mm256_sub_epi32(s02,s13),Pv),tw+7,tws+7);
__m256i rd=shoup8_ptr(_mm256_add_epi32(_mm256_sub_epi32(d02,d13),Pv),tw+7,tws+7);
v0=ra;v1=rb;v2=rc;v3=rd;
};
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)));
if(Fuse32){finish32(vr0,vr1,vr2,vr3);finish32(vr4,vr5,vr6,vr7);}
{
{
__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);
vr0=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr0,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr0=_mm256_min_epu32(vr0,_mm256_sub_epi32(vr0,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr1, 0xB1);
vr1=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr1,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr1=_mm256_min_epu32(vr1,_mm256_sub_epi32(vr1,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr2, 0xB1);
vr2=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr2,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr2=_mm256_min_epu32(vr2,_mm256_sub_epi32(vr2,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr3, 0xB1);
vr3=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr3,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr3=_mm256_min_epu32(vr3,_mm256_sub_epi32(vr3,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr4, 0xB1);
vr4=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr4,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr4=_mm256_min_epu32(vr4,_mm256_sub_epi32(vr4,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr5, 0xB1);
vr5=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr5,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr5=_mm256_min_epu32(vr5,_mm256_sub_epi32(vr5,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr6, 0xB1);
vr6=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr6,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
if constexpr(!LazyOut)vr6=_mm256_min_epu32(vr6,_mm256_sub_epi32(vr6,Pv));
}
}
{
{
__m256i y = _mm256_shuffle_epi32(vr7, 0xB1);
vr7=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(vr7,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
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);
v=_mm256_add_epi32(_mm256_add_epi32(_mm256_xor_si256(v,_mm256_set1_epi64x((long long)0xffffffff00000000ull)),y),
_mm256_set1_epi64x((long long)((u64)(P+1)<<32)));
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 __m256i g_psumv[4];
static u32 g_psums;
static u32 psum_take() {
__m256i v = _mm256_add_epi32(_mm256_add_epi32(g_psumv[0], g_psumv[1]), _mm256_add_epi32(g_psumv[2], g_psumv[3]));
v = _mm256_add_epi32(v, _mm256_permute2x128_si256(v, v, 0x01));
v = _mm256_add_epi32(v, _mm256_shuffle_epi32(v, 0x4E));
v = _mm256_add_epi32(v, _mm256_shuffle_epi32(v, 0xB1));
return g_psums + (u32)_mm256_extract_epi32(v, 0);
}
static void psum_reset() { for (int i = 0; i < 4; i++) g_psumv[i] = _mm256_setzero_si256(); g_psums = 0; }
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 mont8x(__m256i a, __m256i ax, __m256i b, __m256i bx) {
__m256i P64 = _mm256_set1_epi64x((long long)P), NP = _mm256_set1_epi64x((long long)NPI);
__m256i Pv = _mm256_set1_epi32((int)P);
__m256i tlo = _mm256_mul_epu32(a, b), thi = _mm256_mul_epu32(ax, bx);
__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 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 *pf = in + B0 + v0 + (u64)(u0 + 8) * SRF_;
for (int k = 0; k < 8; k++) _mm_prefetch((const char *)(pf + (u64)k * SRF_), _MM_HINT_T0); }
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] = mont8x(ldu256((const void *)(sb + (u64)k * SRF_)), ldu256((const void *)(sb + (u64)k * SRF_ + 1)),
ldu256((const void *)(tb + (long)k * TS)), ldu256((const void *)(tb + (long)k * TS + 1)));
}
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_));
g_psumv[k & 3] = _mm256_add_epi32(g_psumv[k & 3], x);
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; }
const __m256i c4 = _mm256_set1_epi32(4), c4p = _mm256_set1_epi32((int)(P - 4));
for (int rb = 0; rb < TS; rb += 8) {
for (int cb0 = 0; cb0 < TS; cb0 += 8) {
__m256i r[8];
for (int k = 0; k < 8; k++) {
u64 idx = (u64)(r0 + rb + k) * C + c0 + cb0;
if (idx + 8 <= (u64)nsrc) {
__m256i x = ldu256((const void *)(src + idx));
g_psumv[k & 3] = _mm256_add_epi32(g_psumv[k & 3], x);
if (ctr) x = _mm256_min_epu32(_mm256_sub_epi32(x, c4), _mm256_add_epi32(x, c4p));
r[k] = x;
} else if (idx >= (u64)nsrc) {
r[k] = _mm256_setzero_si256();
} else {
u32 tmp[8];
for (int q = 0; q < 8; q++) {
u64 id = idx + q; u32 v = 0;
if (id < (u64)nsrc) { v = src[id]; g_psums += v; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
tmp[q] = v;
}
r[k] = ldu256((const void *)tmp);
}
}
t8x8(r);
const u32 *db = S + (u64)r0 + rb;
for (int j = 0; j < 8; j++) stu256((void *)(db + (long)(cb0 + j) * dstr), r[j]);
}
}
}
}
static void copy_u32_avx(u32 *dst, const u32 *src, long n) {
long i = 0;
for (; i + 8 <= n; i += 8) stu256((void *)(dst + i), ldu256((const void *)(src + i)));
for (; i < n; i++) dst[i] = src[i];
}
static void copy_u32_nt(u32 *dst, const u32 *src, long n) {
long i = 0;
long ka = ((-(long)(((uintptr_t)dst) >> 2)) & 7);
for (; i < ka && i < n; i++) dst[i] = src[i];
for (; i + 8 <= n; i += 8) {
__m256i v = ldu256((const void *)(src + i));
_mm256_stream_si256((__m256i *)(dst + i), v);
}
for (; i < n; i++) dst[i] = src[i];
_mm_sfence();
}
static void untr_tile_avx(const u32 *in, u32 *out, int A0, int B0, int strIn, int strOut) {
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 16) {
__m256i r[8];
for (int k = 0; k < 8; k++)
r[k] = ldu256((const void *)(BROW((long)(A0 + u0 + k) * strIn + B0 + v0)));
t8x8(r);
for (int k = 0; k < 8; k++)
stu256((void *)(TBUF + (v0 + k) * TS + u0), r[k]);
{
int v0b = v0 + 8;
__m256i r2[8];
for (int k = 0; k < 8; k++)
r2[k] = ldu256((const void *)(BROW((long)(A0 + u0 + k) * strIn + B0 + v0b)));
t8x8(r2);
for (int k = 0; k < 8; k++)
stu256((void *)(TBUF + (v0b + k) * TS + u0), r2[k]);
}
}
}
for (int v = 0; v < TS; v++) {
u32 *d = out + (u64)(B0 + v) * strOut + A0;
const u32 *s = TBUF + v * TS;
for (int u = 0; u < TS; u += 8)
_mm256_stream_si256((__m256i *)(d + u), ldu256((const void *)(s + u)));
}
}
static void build_tab2(u32 *tw, u32 *tws, int L, u32 w) {
int off = 0;
for (int len = 2; len <= L; len <<= 1) {
u32 wl = pw(w, (u32)(L / len));
u32 c = 1;
int h = len >> 1;
for (int j = 0; j < h; j++) { tw[off + j] = c; tws[off + j] = shoup_c(c); c = mm(c, wl); }
off += h;
}
}
#pragma GCC pop_options
#endif
static u32 R3F[2], R3FS[2];
static u32 R3I[2], R3IS[2];
static u32 T3F_[2][M3_ + 16], T3FS_[2][M3_ + 16];
#define T3F T3F_[0]
#define T3FS T3FS_[0]
static u32 T3I[2 * M3_], T3IS[2 * M3_];
static u32 tsubFw[M3_ + 8], tsubFsw[M3_ + 8], tsubIw[M3_ + 8], tsubIsw[M3_ + 8];
static u32 TSUB[4][M3_ + 64] __attribute__((aligned(64)));
static u32 TT3I[2*M3_ + 128] __attribute__((aligned(64))), TT3IS[2*M3_ + 128] __attribute__((aligned(64)));
static u32 RTF[NR_ + 64] __attribute__((aligned(64))), RTI[NR_ + 64] __attribute__((aligned(64))), RWS[2*NR_ + 64] __attribute__((aligned(64)));
static void build_radix3(u32 root) {
const u32 H2 = 1006632961u;
u32 w3 = pw(root, (u32)(NC_ / 3));
u32 w32 = mm(w3, w3);
u32 cA = (u32)(P - H2);
u32 cP = mm((u32)(((u64)w3 + P - w32) % P), H2);
R3F[0] = cA; R3F[1] = cP;
R3I[0] = cA; R3I[1] = (cP ? P - cP : 0);
for (int i = 0; i < 2; i++) { R3FS[i] = shoup_c(R3F[i]); R3IS[i] = shoup_c(R3I[i]); }
for (int s = 1; s < 3; s++) {
u32 g = pw(root, (u32)s), gi = pw(g, P - 2u), c = 1, ci = 1;
u32 *pf = T3F_[s-1], *pf2 = T3FS_[s-1];
u32 *pi = T3I + (s - 1) * M3_, *pi2 = T3IS + (s - 1) * M3_;
for (int j = 0; j < (int)M3_; j++) { pf[j] = c; pf2[j] = shoup_c(c); c = mm(c, g);
pi[j] = ci; pi2[j] = shoup_c(ci); ci = mm(ci, gi); }
}
build_tab2(tsubFw, tsubFsw, (int)M3_, pw(root, 3));
build_tab2(tsubIw, tsubIsw, (int)M3_, pw(pw(root, 3), P - 2u));
}
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __attribute__((always_inline)) __m256i vaddm(__m256i a, __m256i b, __m256i Pv) {
__m256i s = _mm256_add_epi32(a, b);
return _mm256_min_epu32(s, _mm256_sub_epi32(s, Pv));
}
static inline __attribute__((always_inline)) __m256i vsubm(__m256i a, __m256i b, __m256i Pv) {
__m256i d = _mm256_sub_epi32(a, b);
return _mm256_min_epu32(d, _mm256_add_epi32(d, Pv));
}
static inline __attribute__((always_inline)) void dft3(__m256i *v, const __m256i *c,
const __m256i *cs, __m256i Pv) {
__m256i A = vaddm(v[1], v[2], Pv);
__m256i B = _mm256_add_epi32(_mm256_sub_epi32(v[1],v[2]),Pv);
const __m256i parity=_mm256_slli_epi32(A,31);
const __m256i base=_mm256_castps_si256(_mm256_blendv_ps(_mm256_castsi256_ps(Pv),
_mm256_castsi256_ps(_mm256_set1_epi32((int)((P-1u)/2u))),_mm256_castsi256_ps(parity)));
__m256i T=_mm256_sub_epi32(base,_mm256_srli_epi32(A,1));
__m256i U = shoup8_broadcast(B, c[1], cs[1]);
__m256i X0 = vaddm(v[0], A, Pv);
__m256i Z = vaddm(v[0], T, Pv);
v[0] = X0; v[1] = _mm256_add_epi32(Z,U); v[2] = _mm256_add_epi32(_mm256_sub_epi32(Z,U),Pv);
}
static void col_dif3_head(u32 *x) {
const __m256i Pv = _mm256_set1_epi32((int)P);
__m256i c[2], cs[2];
for (int i = 0; i < 2; i++) { c[i] = _mm256_set1_epi32((int)R3F[i]); cs[i] = _mm256_set1_epi32((int)R3FS[i]); }
for (int j = 0; j < (int)M3_; j += 8) {
__m256i v[3];
for (int i = 0; i < 3; i++) v[i] = ldu256((const void *)(x + (u64)i * M3_ + j));
dft3(v, c, cs, Pv);
stu256((void *)(x + j), v[0]);
for (int s = 1; s < 3; s++) {
__m256i w = ldu256((const void *)(T3F_[s-1] + j));
__m256i ws = ldu256((const void *)(T3FS_[s-1] + j));
stu256((void *)(x + (u64)s * M3_ + j), shoup8(v[s], w, ws));
}
}
}
static void col_dif3(u32 *x) {
col_dif3_head(x);
for (int s = 0; s < 3; s++) row_dif_avx(x + (u64)s * M3_, (int)M3_, TSUB[0]+16, TSUB[1]+16);
}
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 *)(TT3I + 16 + (s - 1) * M3_ + (J))); \
ws = ldu256((const void *)(TT3IS + 16 + (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_, TSUB[2]+16, TSUB[3]+16);
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_,TSUB[0]+16,TSUB[1]+16);
row_dif_avx(pb,(int)M3_,TSUB[0]+16,TSUB[1]+16);
if(arm==0)row_dit_avx<true>(pa,(int)M3_,TSUB[2]+16,TSUB[3]+16,pb);
else row_dit_avx<true,true>(pa,(int)M3_,TSUB[2]+16,TSUB[3]+16,pb);
}
col_dit3_tail(a);
}
#pragma GCC pop_options
static void row_dit(u32 *x, int L, const u32 *tw, const u32 *tws) {
#if HAVE_AVX2
row_dit_avx(x, L, tw, tws);
return;
#else
(void)tws;
#endif
int off = 0;
for (int len = 2; len <= L; len <<= 1) {
int h = len >> 1;
const u32 *t = tw + off; off += h;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = mm(q[j], t[j]);
u32 s = u + v; if (s >= P) s -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = s; q[j] = d;
}
}
}
}
static void pointwise(u32 *A, const u32 *B, u32 n, u32 sc, u32 sx) {
#if HAVE_AVX2
pointwise_mul_avx(A, B, n, sc, sx);
return;
#else
for (long i = 0; i < (long)n; i++) { long o = (i / (long)(sx - PADC_)) * sx + (i % (long)(sx - PADC_)); A[o] = mm(mm(A[o], B[o]), sc); }
#endif
}
static void row_dif(u32 *x, int L, const u32 *tw, const u32 *tws) {
#if HAVE_AVX2
row_dif_avx(x, L, tw, tws);
return;
#else
(void)tws;
#endif
for (int len = L; len >= 2; len >>= 1) {
int h = len >> 1;
const u32 *t = tw + h - 1;
for (int i = 0; i < L; i += len) {
u32 *p = x + i, *q = x + i + h;
for (int j = 0; j < h; j++) {
u32 u = p[j], v = q[j];
u32 s = u + v; if (s >= P) s -= P;
u32 d = u + P - v; if (d >= P) d -= P;
p[j] = s; q[j] = mm(d, t[j]);
}
}
}
}
static inline void bitrev_pass(u32 *x, int L, const u16 *br) {
for (int i = 0; i < L; i++) {
int k = br[i];
if (i < k) { u32 t = x[i]; x[i] = x[k]; x[k] = t; }
}
}
struct Eng {
int R, C, N, sx;
int blkr;
int r5;
u32 w;
u32 *tfR; u32 *tiR; u32 *tfC; u32 *tiC;
u32 *wsR; u32 *wsC;
u16 *brR; u16 *brC;
};
static void col_dif3_arms(u32 *x) {
for (int s = 0; s < 3; s++) row_dif_avx(x + (u64)s * M3_, (int)M3_, TSUB[0]+16, TSUB[1]+16);
}
static void colfwd(Eng &e, u32 *row) {
#if HAVE_AVX2
if (e.r5) { col_dif3_arms(row); return; }
#endif
row_dif(row, e.C, e.tfC, e.wsC);
}
static void colinv(Eng &e, u32 *row) {
#if HAVE_AVX2
if (e.r5) { col_dit3(row); return; }
#endif
row_dit(row, e.C, e.tiC, e.wsC + e.C);
}
static void load_T(Eng &e, const u32 *src, u32 nsrc, u32 *S, int ctr) {
const int R = e.R, C = e.C;
#if HAVE_AVX2
if (nsrc <= (u32)R * (u32)C) {
for (int r0 = 0; r0 < R; r0 += TS) {
u64 tilerow = (u64)r0 * C;
for (int c0 = 0; c0 < C; c0 += TS) {
u64 firstidx = tilerow + c0;
u64 lastidx = (u64)(r0 + TS - 1) * C + c0 + TS - 1;
if (lastidx < nsrc) { load_tile_avx(src, S, r0, c0, R, C, ctr, 0); continue; }
if (firstidx >= nsrc) { load_tile_avx(src, S, r0, c0, R, C, ctr, 1); continue; }
for (int rb = 0; rb < TS; rb++) {
u64 base = (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u64 id = base + cb;
u32 v = 0;
if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
TBUF[cb * TS + rb] = v;
}
}
for (int cb = 0; cb < TS; cb++) {
u32 *d = S + (u64)(c0 + cb) * R + r0;
const u32 *s = TBUF + cb * TS;
for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
}
}
}
return;
}
#endif
for (int r0 = 0; r0 < R; r0 += TS) {
for (int c0 = 0; c0 < C; c0 += TS) {
u64 lastidx = (u64)(r0 + TS - 1) * C + (c0 + TS - 1);
if (lastidx < nsrc) {
for (int rb = 0; rb < TS; rb++) {
const u32 *s = src + (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u32 v = s[cb];
if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u));
TBUF[cb * TS + rb] = v;
}
}
} else {
for (int rb = 0; rb < TS; rb++) {
u64 base = (u64)(r0 + rb) * C + c0;
for (int cb = 0; cb < TS; cb++) {
u64 id = base + cb;
u32 v = 0;
if (id < nsrc) { v = src[id]; if (ctr) v = (v >= 4u) ? (v - 4u) : (v + (P - 4u)); }
TBUF[cb * TS + rb] = v;
}
}
}
for (int cb = 0; cb < TS; cb++) {
u32 *d = S + (u64)(c0 + cb) * R + r0;
const u32 *s = TBUF + cb * TS;
for (int rb = 0; rb < TS; rb++) d[rb] = s[rb];
}
}
}
}
static u32 REV6[TS];
static u32 g_inv_scale = 1u;
static u32 POWT[1024], POWT64[1024];
static void build_rev6(void) {
for (int i = 0; i < TS; i++) {
int r = 0;
for (int b = 0; b < 6; b++) if (i & (1 << b)) r |= 1 << (5 - b);
REV6[i] = (u32)r;
}
}
static inline u32 revbits(u32 x, int lg) {
u32 r = 0;
for (int b = 0; b < lg; b++) if (x & (1u << b)) r |= 1u << (lg - 1 - b);
return r;
}
template<int perm>
static void corner_generic(const u32 *in, u32 *out, int nA, int nB, int strIn, int strOut, u32 root, int split) {
u32 wv[TS];
{ u32 c = 1; for (int v = 0; v < TS; v++) { wv[v] = c; c = mm(c, root); } }
const u32 w64 = pw(root, TS);
int lgB = 0, lgA = 0;
while ((1 << lgB) < nB) lgB++;
while ((1 << lgA) < nA) lgA++;
const int S1 = nB / TS, S2 = nA / TS;
const int HB = lgB - 6, HA = lgA - 6;
u32 RV2[TS], RV2b[TS], PW2[TS], cu[TS], cu64[TS], base[TS];
u32 cu64_plain[TS],cu64_shoup[TS];
if (perm == 1) {
u32 g = pw(root, (u32)S1);
{ u32 gM = mm(g, R32), c = R32; for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; c = mont1(c, gM); } }
{ u32 rootM = mm(root, R32), p = R32; for (int t = 0; t < S1; t++) { POWT[t] = p; p = mont1(p, rootM); } }
} else if (perm == 2) {
u32 g = pw(root, (u32)S2), g64 = pw(g, 64);
const u32 rootM2 = mm(root, R32), gM = mm(g, R32), g64M = mm(g64, R32);
{ u32 c = R32, c64 = R32;
for (int t = 0; t < TS; t++) { RV2[REV6[t]] = c; RV2b[REV6[t]] = c64; c = mont1(c, gM); c64 = mont1(c64, g64M); } }
{ u32 p = R32; for (int t = 0; t < S2; t++) { POWT[t] = p; p = mont1(p, rootM2); } }
{ u32 q = pw(root, 64), qM = mm(q, R32), p = R32; for (int t = 0; t < S2; t++) { POWT64[t] = p; p = mont1(p, qM); } }
}
u32 wr = 1;
const u32 INVSCALE = g_inv_scale;
for (int A0 = 0; A0 < nA; A0 += TS) {
u32 wr64 = pw(wr, TS);
if (perm == 1) {
u32 g = pw(wr, (u32)S1), gM2 = mm(g, R32), c = R32;
for (int t = 0; t < TS; t++) { PW2[REV6[t]] = c; c = mont1(c, gM2); }
} else if (perm == 2) {
int hi = A0 / TS;
u32 rH = POWT[revbits((u32)hi, HA)], r64H = POWT64[revbits((u32)hi, HA)];
const u32 INVSCALE_M = mm(INVSCALE, R32);
for (int u = 0; u < TS; u++) {
cu[u] = mont1(RV2[u], rH); cu64[u] = mont1(RV2b[u], r64H); base[u] = INVSCALE_M;
cu64_plain[u]=mont1(cu64[u],1u); cu64_shoup[u]=shoup_c(cu64_plain[u]);
}
} else {
u32 c = 1; for (int v = 0; v < TS; v++) { PW2[v] = c; c = mm(c, wr); }
}
u32 Vloc = 1;
u32 Vr = 1;
for (int B0 = 0; B0 < nB; B0 += TS) {
if (perm == 1) {
u32 H = revbits((u32)(B0 / TS), HB);
u32 wrH = pw(wr, H), rootH = POWT[H];
u32 RV[TS];
#if HAVE_AVX2
scale_row_mont(RV, RV2, rootH);
scale_row_mont(TW, PW2, mm(wrH, R32));
twiddle_rows_mont(TW, RV);
#else
for (int v = 0; v < TS; v++) RV[v] = mont1(RV2[v], rootH);
for (int v = 0; v < TS; v++) TW[v] = mont1(PW2[v], mm(wrH, R32));
for (int u = 1; u < TS; u++) {
const u32 *prv = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v++) cur[v] = mont1(prv[v], RV[v]);
}
#endif
} else if (perm == 2) {
if (B0 == 0) {
for (int u = 0; u < TS; u++) {
u32 *cur = TW + u * TS;
cur[0] = base[u];
for (int v = 1; v < TS; v++) cur[v] = mont1(cur[v - 1], cu[u]);
}
} else {
for (int u = 0; u < TS; u++) {
u32 *cur = TW + u * TS;
#if HAVE_AVX2
scale_row_prepared(cur, cur, cu64_plain[u], cu64_shoup[u]);
#else
u32 d = cu64[u];
for (int v = 0; v < TS; v++) cur[v] = mont1(cur[v], d);
#endif
}
}
} else {
u32 RV[TS];
for (int v = 0; v < TS; v++) RV[v] = mm(Vloc, wv[v]);
for (int v = 0; v < TS; v++) TW[v] = mm(Vr, PW2[v]);
for (int u = 1; u < TS; u++) {
const u32 *prv = TW + (u - 1) * TS;
u32 *cur = TW + u * TS;
for (int v = 0; v < TS; v++) cur[v] = mm(prv[v], RV[v]);
}
}
#if HAVE_AVX2
corner_tile_avx(in, out, A0, B0, strIn, strOut, split);
#else
for (int u = 0; u < TS; u++) {
const u32 *s = in + (u64)(A0 + u) * strIn + B0;
const u32 *tw = TW + u * TS;
for (int v = 0; v < TS; v++) TBUF[u * TS + v] = mont1(s[v], tw[v]);
}
for (int v = 0; v < TS; v++) {
u32 *d = out + (u64)(B0 + v) * strOut + A0;
for (int u = 0; u < TS; u++) d[u] = TBUF[u * TS + v];
}
#endif
Vr = mm(Vr, wr64);
Vloc = mm(Vloc, w64);
}
wr = mm(wr, w64);
}
}
static inline void corner_fwd(Eng &e, const u32 *S, u32 *X, u32 root) {
corner_generic<1>(S, X, e.C, e.R, e.R, e.sx, root, 0);
}
static inline void corner_inv(Eng &e, const u32 *X, u32 *S, u32 root) {
corner_generic<2>(X, S, e.R, e.C, e.sx, e.R, root, 1);
}
static void untranspose(Eng &e, const u32 *S, u32 *dst) {
const int R = e.R, C = e.C;
#if HAVE_AVX2
for (int c0 = 0; c0 < C; c0 += TS)
for (int r0 = 0; r0 < R; r0 += TS)
untr_tile_avx(S, dst, c0, r0, R, e.sx);
return;
#endif
for (int c0 = 0; c0 < C; c0 += TS) {
for (int r0 = 0; r0 < R; r0 += TS) {
for (int u = 0; u < TS; u++) {
const u32 *s = S + (u64)(c0 + u) * R + r0;
u32 *d = TBUF + u * TS;
for (int v = 0; v < TS; v++) d[v] = s[v];
}
for (int v = 0; v < TS; v++) {
u32 *d = dst + (u64)(r0 + v) * e.sx + c0;
for (int u = 0; u < TS; u++) d[u] = TBUF[u * TS + v];
}
}
}
}
#define N02I_SB (TS * (NRMAX_ + PADB_))
static u32 SBLK[3 * N02I_SB] __attribute__((aligned(64)));
static u32 TILE3[3][TS * TS + 32] __attribute__((aligned(64)));
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static u32 GRV2[TS], GPOWT[NR_/TS];
static u32 GRV[TS][TS];
static __m256i GSF[TS][TS/8], GSS[TS][TS/8];
static u32 GWRHM[256][3][NR_/TS], GPW2[256][3][TS];
static int g_prep_done=0; static u32 g_prep_root=0;
static void fwd_prep1(u32 root){
if(g_prep_done && g_prep_root==root) return;
g_prep_done=1; g_prep_root=root;
const int S1=NR_/TS;
{ u32 g=pw(root,(u32)S1), gM=mm(g,R32), c=R32;
for(int t=0;t<TS;t++){ GRV2[REV6[t]]=c; c=mont1(c,gM); } }
{ u32 rootM=mm(root,R32), p=R32;
for(int t=0;t<S1;t++){ GPOWT[t]=p; p=mont1(p,rootM); } }
for(int H=0;H<S1;++H){ scale_row_mont(GRV[H],GRV2,GPOWT[H]); prepare_shared_row_factors(GRV[H],GSF[H],GSS[H]); }
const int AM=NC_/3;
for(int Ai=0,Ab=0;Ab<AM;Ab+=TS,++Ai){
for(int p=0;p<3;p++){
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){GWRHM[Ai][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++){GPW2[Ai][p][REV6[t]]=c; c=mont1(c,gM2);} }
}
}
}
#pragma GCC pop_options
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;
fwd_prep1(root);
const bool halfzero=(nB==4096 && (u64)nsrc<=(u64)(nB/2)*nA);
const int AM = nA / 3;
for (int Ab = 0; Ab < AM; Ab += TS) {
const int Ai = Ab / TS;
for (int p = 0; p < 3; p++) {
u32 *S = SBLK + (u64)p * N02I_SB;
const int A0 = Ab + p * AM;
load_block(src, nsrc, S, A0, halfzero?nB/2:nB, nA, strIn, ctr);
for (int u=0;u<TS;++u){
if(halfzero)row_dif_avx<true,true>(S+(u64)u*strIn,nB,e.tfR,e.wsR);
else row_dif(S+(u64)u*strIn,nB,e.tfR,e.wsR);
}
}
for (int B0 = 0; B0 < nB; B0 += TS) {
u32 H = revbits((u32)(B0 / TS), HB);
const __m256i *shared_factor=GSF[H]; const __m256i *shared_shoup=GSS[H];
for (int p = 0; p < 3; p++) {
scale_row_mont(TW, GPW2[Ai][p], GWRHM[Ai][p][H]);
twiddle_rows_prepared(TW, shared_factor, shared_shoup);
corner_tile_avx_to(SBLK + (u64)p * N02I_SB, TILE3[p], B0, strIn);
}
corner_head_tile_avx(TILE3[0],TILE3[1],TILE3[2],X,B0,Ab,split);
}
}
}
static void fwd(Eng &e, const u32 *src, u32 nsrc, u32 *S, u32 *X, int split, int ctr) {
#ifdef OLD_FWD
load_T(e, src, nsrc, S, ctr);
for (int j0 = 0; j0 < e.C; j0++) {
u32 *row = S + (u64)j0 * e.R;
row_dif(row, e.R, e.tfR, e.wsR);
}
corner_fwd(e, S, X, e.w);
#else
(void)S;
{ corner_fwd_stream(e, src, nsrc, X, split, ctr); }
#endif
if (!g_fuse4) { for (int k0 = 0; k0 < e.R; k0++)
colfwd(e, split ? BROW((long)k0 * e.sx) : X + (u64)k0 * e.sx); }
}
static void inv(Eng &e, u32 *X, u32 *S) { for (int k0 = 0; k0 < e.R; k0++) colinv(e, X + (u64)k0 * e.sx);
corner_inv(e, X, S, pw(e.w, P - 2u));
for (int j0 = 0; j0 < e.C; j0++) {
u32 *row = BROW((long)j0 * e.R);
row_dit(row, e.R, e.tiR, e.wsR + e.R);
}
}
static u32 INV_TW[NRMAX_ * TS] __attribute__((aligned(64)));
static u32 INV_STEP[NRMAX_], INV_STEP_SHOUP[NRMAX_];
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
template<bool Advance>
static inline __attribute__((always_inline)) void inverse_twiddle_consumer(Eng &e, u32 *X, int r0, int c0) {
(void)e;
const long SXF_ = (long)CSX_, BLF_ = (long)(NR_ + PADB_);
for (int u0 = 0; u0 < TS; u0 += 8) {
for (int v0 = 0; v0 < TS; v0 += 8) {
{ const u32 *pf = X + (u64)(r0 + u0 + 8) * SXF_ + c0 + v0;
for (int k = 0; k < 8; k++) _mm_prefetch((const char *)(pf + (u64)k * SXF_), _MM_HINT_T0); }
__m256i rows[8];
{
const u32 *sbase = X + (u64)(r0 + u0) * SXF_ + c0 + v0;
u32 *tbase = INV_TW + (r0 + u0) * (long)TS + v0;
u32 *sblkbase = SBLK + (u64)v0 * BLF_ + r0 + u0;
for (int k = 0; k < 8; ++k) {
const u32 *src = sbase + (u64)k * SXF_;
u32 *tw = tbase + (long)k * TS;
__m256i w = ldu256(tw);
if constexpr(Advance) {
const int j = r0 + u0 + k;
stu256(tw, shoup8_broadcast(w, _mm256_set1_epi32((int)INV_STEP[j]),
_mm256_set1_epi32((int)INV_STEP_SHOUP[j])));
}
rows[k] = mont8(ldu256(src), w);
}
t8x8(rows);
for (int k = 0; k < 8; ++k)
stu256(sblkbase + (u64)k * BLF_, rows[k]);
}
}
}
}
static void inv_natural_stream(Eng &e, u32 *X) {
if (!g_fuse4)
for (int k = 0; k < e.R; ++k) colinv(e, X + (u64)k * e.sx);
const u32 root = pw(e.w, P - 2u);
const u32 seed = mm(mm(g_inv_scale, R32), R32);
int lg = 0; while ((1 << lg) < e.R) ++lg;
for (int k = 0; k < e.R; ++k) {
const u32 factor = pw(root, revbits(k, lg));
const u32 factorM = mm(factor, R32);
INV_STEP[k] = pw(factor, TS);
INV_STEP_SHOUP[k] = shoup_c(INV_STEP[k]);
u32 *tw = INV_TW + k * TS;
tw[0] = seed;
for (int v = 1; v < TS; ++v) tw[v] = mont1(tw[v - 1], factorM);
}
for (int c0 = 0; c0 < e.C; c0 += TS) {
for (int r0 = 0; r0 < e.R; r0 += TS) {
inverse_twiddle_consumer<true>(e,X,r0,c0);
}
for (int v = 0; v < TS; ++v)
row_dit(SBLK + (u64)v * e.blkr, e.R, e.tiR, e.wsR + e.R);
{
const long SXF_ = (long)CSX_, BLF_ = (long)(NR_ + PADB_);
for (int r0 = 0; r0 < (int)NR_; r0 += TS) {
for (int v0 = 0; v0 < TS; v0 += 8) {
{
const u32 *sbase = SBLK + (u64)(r0 + v0);
u32 *dbase = X + (u64)(r0 + v0) * SXF_ + c0;
for (int u0 = 0; u0 < TS; u0 += 16) {
const u32 *sp = sbase + (u64)u0 * BLF_;
__m256i rows[8];
for (int k = 0; k < 8; ++k)
rows[k] = ldu256(sp + (u64)k * BLF_);
t8x8(rows);
for (int k = 0; k < 8; ++k) stu256(TBUF+k*8,rows[k]);
for (int k = 0; k < 8; ++k)
rows[k] = ldu256(sp + (u64)(8 + k) * BLF_);
t8x8(rows);
for (int k = 0; k < 8; ++k) {
u32 *dst = dbase + (u64)k * SXF_ + u0;
_mm256_stream_si256((__m256i *)dst, ldu256(TBUF+k*8));
_mm256_stream_si256((__m256i *)(dst+8), rows[k]);
}
}
}
}
}
}
}
_mm_sfence();
}
static void fuse4_row_avx(u32 *pa, const u32 *pb, int n) {
int i = 0;
for (; i + 8 <= n; i += 8) {
__m256i x = ldu256((const void *)(pa + i));
__m256i y = ldu256((const void *)(pb + i));
stu256((void *)(pa + i), shoup8(x, y, shoup_of(y)));
}
for (; i < n; i++) pa[i] = mm(pa[i], pb[i]);
}
#pragma GCC pop_options
static void eng_init(Eng &e, int lgR, int nC, u32 w, int sx,
u32 *tfR, u32 *tiR, u32 *tfC, u32 *tiC, u16 *brR, u16 *brC,
u32 *wsR, u32 *wsC) {
e.R = 1 << lgR; e.C = nC; e.N = e.R * e.C; e.w = w; e.sx = sx;
e.blkr = e.R + PADB_;
e.r5 = (nC == (int)NC_);
e.tfR = tfR; e.tiR = tiR; e.tfC = tfC; e.tiC = tiC; e.brR = brR; e.brC = brC;
e.wsR = wsR; e.wsC = wsC;
u32 wR = pw(w, (u32)e.C), wC = pw(w, (u32)e.R);
#if HAVE_AVX2
build_tab2(tfR, e.wsR, e.R, wR); build_tab2(tiR, e.wsR + e.R, e.R, pw(wR, P - 2u));
if (!e.r5) { build_tab2(tfC, e.wsC, e.C, wC); build_tab2(tiC, e.wsC + e.C, e.C, pw(wC, P - 2u)); }
#else
build_tab(tfR, e.R, wR); build_tab(tiR, e.R, pw(wR, P - 2u));
build_tab(tfC, e.C, wC); build_tab(tiC, e.C, pw(wC, P - 2u));
#endif
build_br(brR, e.R);
if (!e.r5) build_br(brC, e.C);
build_rev6();
}
#pragma GCC push_options
#pragma GCC target("avx2,tune=skylake")
static inline __attribute__((always_inline)) __m256i scan8(__m256i v, __m256i &carry) {
v = _mm256_add_epi32(v, _mm256_slli_si256(v, 4));
v = _mm256_add_epi32(v, _mm256_slli_si256(v, 8));
__m256i t = _mm256_permute2x128_si256(_mm256_setzero_si256(), v, 0x20);
t = _mm256_shuffle_epi32(t, 0xFF);
v = _mm256_add_epi32(v, t);
v = _mm256_add_epi32(v, carry);
carry = _mm256_permutevar8x32_epi32(v, _mm256_set1_epi32(7));
return v;
}
#define CMBSTEP() \
{ u32 d; \
if (PH == 0) d = 4u * (a[k] + b[k]) - 16u; \
else if (PH == 1) { if (nm_lt) d = 4u * (b[k] - b[k - 1 - n]); \
else d = 4u * (a[k] - a[k - 1 - m]); } \
else d = 16u - 4u * (a[k - 1 - m] + b[k - 1 - n]); \
SUMr += d; \
u32 tab; \
if (TB == 1) { u32 x = A[AOFF(k)], y = c[N1_ + k]; tab = (x >= y) ? (x - y) : (x + (P - y)); } \
else if (TB == 2) tab = A[AOFF(k)]; \
else tab = c[k]; \
i32 ts = (i32)tab; \
if (ts > (i32)(P >> 1)) ts -= (i32)P; \
c[k] = (u32)((u32)ts + SUMr); \
}
template<int PH, int TB, int NT>
static void cmbv(long k0, long k1, const u32 *a, const u32 *b, u32 *A, u32 *c,
long n, long m, long nmax, u32 &SUMr, int nm_lt)
{
if (k0 > k1) return;
long k = k0;
{
long ka = k + ((-(long)k) & 7);
long kend = (ka < k1 + 1) ? ka : k1 + 1;
if (ka > k) { const __m256i Pv0 = _mm256_set1_epi32((int)P); (void)Pv0; }
for (; k < kend; k++) { CMBSTEP(); }
}
const __m256i Pv = _mm256_set1_epi32((int)P);
const __m256i Phv = _mm256_set1_epi32((int)(P >> 1));
#if HAVE_AVX2
const __m256i c16 = _mm256_set1_epi32(16);
const __m256i idx7 = _mm256_set1_epi32(7);
__m256i carry = _mm256_set1_epi32((int)SUMr);
if (NT) {
carry = _mm256_set1_epi32((int)SUMr);
}
for (; k + 7 <= k1; k += 8) {
__m256i d;
if (PH == 0) {
__m256i va = ldu256(a + k), vb = ldu256(b + k);
d = _mm256_sub_epi32(_mm256_slli_epi32(_mm256_add_epi32(va, vb), 2), c16);
} else if (PH == 1) {
if (nm_lt) { __m256i v1 = ldu256(b + k), v0 = ldu256(b + k - 1 - n);
d = _mm256_slli_epi32(_mm256_sub_epi32(v1, v0), 2); }
else { __m256i v1 = ldu256(a + k), v0 = ldu256(a + k - 1 - m);
d = _mm256_slli_epi32(_mm256_sub_epi32(v1, v0), 2); }
} else {
__m256i va = ldu256(a + k - 1 - m), vb = ldu256(b + k - 1 - n);
d = _mm256_sub_epi32(c16, _mm256_slli_epi32(_mm256_add_epi32(va, vb), 2));
}
d = scan8(d, carry);
__m256i tab;
if (TB == 1) {
__m256i x = ldu256(A + AOFF(k)), y = ldu256(c + N1_ + k);
__m256i dd = _mm256_sub_epi32(x, y);
__m256i mk = _mm256_cmpgt_epi32(y, x);
tab = _mm256_add_epi32(dd, _mm256_and_si256(mk, Pv));
} else if (TB == 2) {
tab = ldu256(A + AOFF(k));
} else {
tab = ldu256(c + k);
}
__m256i ts = _mm256_sub_epi32(tab, _mm256_and_si256(_mm256_cmpgt_epi32(tab, Phv), Pv));
__m256i out = _mm256_add_epi32(ts, d);
if (NT && (((uintptr_t)(c + k)) & 31) == 0) _mm256_stream_si256((__m256i *)(c + k), out);
else stu256(c + k, out);
}
SUMr = (u32)_mm256_extract_epi32(carry, 0);
#endif
for (; k <= k1; k++) { CMBSTEP(); }
}
static void cmbv_pair_exact(long n, const u32 *a, const u32 *b, u32 *A, u32 *c, u32 S1T) {
const __m256i Pv = _mm256_set1_epi32((int)P);
const __m256i c16 = _mm256_set1_epi32(16);
const __m256i Phv = _mm256_set1_epi32((int)(P >> 1));
const __m256i sp = _mm256_set1_epi32((int)S1T);
__m256i carry = _mm256_setzero_si256(), prev = _mm256_setzero_si256();
long k = 0; int have = 0;
auto emit=[&](__m256i tab,__m256i tab2) {
__m256i va = ldu256(a + k), vb = ldu256(b + k);
__m256i d = _mm256_sub_epi32(_mm256_slli_epi32(_mm256_add_epi32(va, vb), 2), c16);
d = scan8(d, carry);
__m256i ts = _mm256_sub_epi32(tab, _mm256_and_si256(_mm256_cmpgt_epi32(tab, Phv), Pv));
__m256i out = _mm256_add_epi32(ts, d);
if ((((uintptr_t)(c + k)) & 31) == 0) _mm256_stream_si256((__m256i *)(c + k), out); else stu256(c + k, out);
long k2 = k + n + 1;
__m256i ts2 = _mm256_sub_epi32(tab2, _mm256_and_si256(_mm256_cmpgt_epi32(tab2, Phv), Pv));
__m256i out2 = _mm256_add_epi32(ts2, _mm256_sub_epi32(sp, d));
if (!have) { u32 tmp[8]; _mm256_storeu_si256((__m256i *)tmp, out2); for (int q = 0; q < 7; q++) c[k2 + q] = tmp[q]; have = 1; }
else {
__m256i t = _mm256_permute2x128_si256(prev, out2, 0x21);
__m256i u = _mm256_alignr_epi8(out2, t, 12);
long idx = k2 - 1;
if ((((uintptr_t)(c + idx)) & 31) == 0) _mm256_stream_si256((__m256i *)(c + idx), u); else stu256(c + idx, u);
}
prev = out2;
};
while(k+8<=n) {
const long k2=k+n+1;
const long col0=k%(long)NC_,col1=k2%(long)NC_;
long runs=((long)NC_-col0)/8;
const long runs1=((long)NC_-col1)/8;
if(runs1<runs)runs=runs1;
if((n-k)/8<runs)runs=(n-k)/8;
const u32 *p0=A+AOFF(k),*p1=A+AOFF(k2);
if(runs) {
const long end=k+8*runs;
for(;k<end;k+=8,p0+=8,p1+=8)emit(ldu256(p0),ldu256(p1));
} else {
u32 tmp[8];for(int q=0;q<8;++q)tmp[q]=A[AOFF(k2+q)];
emit(ldu256(p0),ldu256(tmp));k+=8;
}
}
u32 SUM = (u32)_mm256_extract_epi32(carry, 0);
if (have) { u32 tmp[8]; _mm256_storeu_si256((__m256i *)tmp, prev); c[2 * n] = tmp[7]; }
for (; k <= n; k++) {
u32 tab = A[AOFF(k)];
i32 ts = (i32)tab; if (ts > (i32)(P >> 1)) ts -= (i32)P;
SUM += 4u * (a[k] + b[k]) - 16u;
c[k] = (u32)((u32)ts + SUM);
}
}
#pragma GCC pop_options
static u32 Aarr[NR_ * CSX_] __attribute__((aligned(64)));
#if (NR_ * CSX_) > BSTAIL_SPLIT
#define BTAILN_ ((int)((NR_ * CSX_) - BSTAIL_SPLIT))
#else
#define BTAILN_ 1
#endif
static u32 BTAIL[BTAILN_] __attribute__((aligned(64)));
static Eng E1, E2;
static u32 t1fR[NR_], t1iR[NR_], t1fC[1], t1iC[1];
static u16 b1R[NR_], b1C[1];
static u32 w1R[2 * NR_], w1C[2];
static u32 t2fR[NR2_], t2iR[NR2_], t2fC[NC2_], t2iC[NC2_];
static u16 b2R[NR2_], b2C[NC2_];
static u32 w2R[2 * NR2_], w2C[2 * NC2_];
static int inited = 0;
void poly_multiply(unsigned *a, int n, unsigned *b, int m, unsigned *c)
{
IM = ~0ull / (u64)P;
if (!inited) {
inited = 1;
{ u64 inv = 1; for (int i = 0; i < 5; i++) inv = inv * (2 - (u64)P * inv); NPI = (u32)(-(u64)inv); }
u32 w1 = pw(G, (u32)((P - 1) / N1_));
eng_init(E1, LGR_, (int)NC_, w1, (int)CSX_, t1fR, t1iR, t1fC, t1iC, b1R, b1C, w1R, w1C);
memcpy(RTF+16, t1fR, (size_t)NR_*4); memcpy(RTI+16, t1iR, (size_t)NR_*4); memcpy(RWS+16, w1R, (size_t)(2*NR_)*4);
E1.tfR=RTF+16; E1.tiR=RTI+16; E1.wsR=RWS+16;
E1.r5 = 1;
#if HAVE_AVX2
build_radix3(pw(w1, (u32)E1.R));
memcpy(TT3I+16, T3I, (size_t)(2*M3_)*4); memcpy(TT3IS+16, T3IS, (size_t)(2*M3_)*4);
memcpy(TSUB[0]+16, tsubFw, (size_t)M3_*4); memcpy(TSUB[1]+16, tsubFsw, (size_t)M3_*4);
memcpy(TSUB[2]+16, tsubIw, (size_t)M3_*4); memcpy(TSUB[3]+16, tsubIsw, (size_t)M3_*4);
#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 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) {
static bool secondary_ready=false;
if(!secondary_ready){
const u32 w2=pw(G,(u32)((P-1)/N2_));
eng_init(E2,LGR2_,(1<<LGC2_),w2,(int)CSX2_,t2fR,t2iR,t2fC,t2iC,b2R,b2C,w2R,w2C);
E2.r5=0;secondary_ready=true;
}
const u32 sc2=pw((u32)(N2_%P),P-2u);
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;
psum_reset();
fwd(E1, (const u32 *)a, (u32)n + 1, S, A, 0, ctr);
u32 PA = psum_take();
psum_reset();
fwd(E1, (const u32 *)b, (u32)m + 1, S, B, 1, ctr);
u32 PB = psum_take();
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;
if (n == m && n >= 0 && ((n & 7) == 0) && M < 0) {
cmbv_pair_exact((long)n, a, b, Aarr, c, 4u * (PA + PB) - 16u * (u32)((long)n + 1));
} else {
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<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);
}
}
#else
for (long k = 0; k <= nmax; k++) c[k] = c[k];
#endif
_mm_sfence();
(void)E2;
}
| Compilation | N/A | N/A | Compile OK | Score: N/A | 显示更多 |
| Testcase #1 | 3.191 s | 1585 MB + 164 KB | Accepted | Score: 100 | 显示更多 |