/* [seat13e 2026-10-31] EDITED copy of the priced 141204 body (sha256 5c7352b98b70...): fused two-stream forward.
* Change: corner_fwd_stream2/fwd2 build each (Ai,p,H) TW table ONCE and apply it to both the A and B SBLK tiles
* (one prep + two corner_tile_avx_to per tile instead of prep+corner twice); the psum accumulators are split
* (g_psumv/g_psumv2 via g_psum_target) so the per-input prefix sums PA/PB are unchanged; colfwd order A then B
* is preserved. The A/B passes' per-tile TW build is the only duplicated work removed; values are unchanged
* (four-mode gate exact). All inherited citation/reference notices are retained verbatim below. */
/* [S8-LC8K] 2026-10-31: 64B line-complete NT stores; rest = #140884 verbatim; the byte-identical wording below is superseded. */
/* Adopted: duck.ac user saffah_codex_6a_agg3, submission #140884 <https://duck.ac/submission/140884> (1002e8, 3.190957753 s Accepted). Body below is that public submission verbatim; all inherited author/license references retained. No separate license declared; site-public submission used by named citation (account+URL+used content). Purpose: fastest-public baseline, byte-exact verified, no new change claimed. */
/* 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.
*/
/* 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.
*/
/* 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.
*/
/* References:
* saffah_codex_6a_agg3 https://duck.ac/submission/139724: positive forward terminal32 fusion, resident subtrees and inherited full reference chain below.
*/
/* 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.
*/
/* 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.
*/
// ===== 引用区 =====
// [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.
* 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.
* 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.
* 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.
*/
/* 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.
* 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。
* 本发目的与改动见下（单变量，试验性；输出逐位不变）。
*/
/* ===== 引用区 =====
* 本件正文 = 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
* 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));
}
/* [S8-LC8K] line-complete NT: one 64B line = its two 32B halves back to back */
#define NT2(p,a,b) {_mm256_stream_si256((__m256i *)(p),(a));_mm256_stream_si256((__m256i *)(p)+1,(b));}
#define NTC(d,s) {_mm256_stream_si256((__m256i *)(d),ldu256((const void *)(s)));_mm256_stream_si256((__m256i *)(d)+1,ldu256((const void *)((s)+8)));}
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], g_psumv2[4]; static __m256i *g_psum_target = g_psumv;
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 u32 psum_take2(){__m256i v=_mm256_add_epi32(_mm256_add_epi32(g_psumv2[0],g_psumv2[1]),_mm256_add_epi32(g_psumv2[2],g_psumv2[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 (u32)_mm256_extract_epi32(v,0);}
static void psum_reset2(){for(int i=0;i<4;i++)g_psumv2[i]=_mm256_setzero_si256();}
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) & 63u) == 0u) {
for (int u = 0; u < TS; u += 16) NTC(d+u,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) & 63u) == 0u) {
for (int u = 0; u < TS; u += 16) NTC(d+u,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_psum_target[k & 3] = _mm256_add_epi32(g_psum_target[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_psum_target[k & 3] = _mm256_add_epi32(g_psum_target[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)) & 15);
for (; i < ka && i < n; i++) dst[i] = src[i];
for (; i + 16 <= n; i += 16) NTC(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 += 16) NTC(d+u,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);
auto hf=[&](int u,__m256i&o0,__m256i&o1,__m256i&o2){
__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);
o0=vv[0];
o1=shoup8(vv[1],ldu256((const void *)(T3F+Ab+u)),ldu256((const void *)(T3FS+Ab+u)));
o2=shoup8(vv[2],ldu256((const void *)(T3F_[1]+Ab+u)),ldu256((const void *)(T3FS_[1]+Ab+u)));
};
for (int u = 0; u < TS; u += 16) {
__m256i p0,p1,p2,q0,q1,q2; hf(u,p0,p1,p2); hf(u+8,q0,q1,q2);
if((((uintptr_t)d0|(uintptr_t)d1|(uintptr_t)d2)&63u)==0u){
NT2(d0+u,p0,q0);NT2(d1+u,p1,q1);NT2(d2+u,p2,q2);
}else{
stu256(d0+u,p0);stu256(d0+u+8,q0);stu256(d1+u,p1);
stu256(d1+u+8,q1);stu256(d2+u,p2);stu256(d2+u+8,q2);
}
}
}
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))), SBLK2[3*N02I_SB] __attribute__((aligned(64)));
static u32 TILE3[3][TS*TS+32] __attribute__((aligned(64))), TILE3B[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);
g_psum_target=g_psumv;
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 corner_fwd_stream2(Eng&e,const u32*srcA,u32 nsrcA,u32*XA,int splitA,int ctrA,const u32*srcB,u32 nsrcB,u32*XB,int splitB,int ctrB){
const int nA=(int)NC_,nB=(int)NR_,strIn=(int)(NR_+PADB_);const u32 root=e.w;int lgB=0;while((1<<lgB)<nB)lgB++;
fwd_prep1(root);const bool hzA=(nB==4096&&(u64)nsrcA<=(u64)(nB/2)*nA),hzB=(nB==4096&&(u64)nsrcB<=(u64)(nB/2)*nA);
const int AM=nA/3,HB=lgB-6;
for(int Ab=0;Ab<AM;Ab+=TS){const int Ai=Ab/TS;
for(int p=0;p<3;p++){const int A0=Ab+p*AM;u32*SA=SBLK+(u64)p*N02I_SB;g_psum_target=g_psumv;
load_block(srcA,nsrcA,SA,A0,hzA?nB/2:nB,nA,strIn,ctrA);
for(int u=0;u<TS;++u){if(hzA)row_dif_avx<true,true>(SA+(u64)u*strIn,nB,e.tfR,e.wsR);else row_dif(SA+(u64)u*strIn,nB,e.tfR,e.wsR);}
u32*SB=SBLK2+(u64)p*N02I_SB;g_psum_target=g_psumv2;
load_block(srcB,nsrcB,SB,A0,hzB?nB/2:nB,nA,strIn,ctrB);
for(int u=0;u<TS;++u){if(hzB)row_dif_avx<true,true>(SB+(u64)u*strIn,nB,e.tfR,e.wsR);else row_dif(SB+(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*sf=GSF[H],*ss=GSS[H];
for(int p=0;p<3;p++){scale_row_mont(TW,GPW2[Ai][p],GWRHM[Ai][p][H]);twiddle_rows_prepared(TW,sf,ss);
corner_tile_avx_to(SBLK+(u64)p*N02I_SB,TILE3[p],B0,strIn);corner_tile_avx_to(SBLK2+(u64)p*N02I_SB,TILE3B[p],B0,strIn);}
corner_head_tile_avx(TILE3[0],TILE3[1],TILE3[2],XA,B0,Ab,splitA);corner_head_tile_avx(TILE3B[0],TILE3B[1],TILE3B[2],XB,B0,Ab,splitB);}}}
static void fwd2(Eng&e,const u32*srcA,u32 nsrcA,u32*XA,int splitA,int ctrA,const u32*srcB,u32 nsrcB,u32*XB,int splitB,int ctrB){
corner_fwd_stream2(e,srcA,nsrcA,XA,splitA,ctrA,srcB,nsrcB,XB,splitB,ctrB);
if(!g_fuse4){for(int k0=0;k0<e.R;k0++)colfwd(e,splitA?BROW((long)k0*e.sx):XA+(u64)k0*e.sx);
for(int k0=0;k0<e.R;k0++)colfwd(e,splitB?BROW((long)k0*e.sx):XB+(u64)k0*e.sx);}}
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();psum_reset2();
fwd2(E1,(const u32*)a,(u32)n+1,A,0,ctr,(const u32*)b,(u32)m+1,B,1,ctr);
u32 PA=psum_take();u32 PB=psum_take2();
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.381 s | 1588 MB + 236 KB | Accepted | Score: 100 | 显示更多 |