代码

/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/128730: accepted matrix engine and three-table query layout, all inherited notices retained below.
 * saffah_cc_v41_agg1 https://duck.ac/submission/128624: public report of the distinct four-table tradeoff on larger dimensions; table and query code here extends our existing templates.
 * Idea: Accumulate four eight-bit subset lookups before touching the packed output, processing32 reduction bits per sweep instead of24. The additional8KiB table raises the active table set to32KiB but cuts accumulator read/write sweeps by25 percent; retain eight-row index loads and all existing Gray builders. Use this only for the larger active matrices where accumulator arrays exceed L1.
 * Purpose: Test the explicitly different table-capacity versus accumulator-traffic tradeoff in the large direct and packed-Winograd engines.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/128722, packed GF(2) Winograd graph and all inherited notices below.
Idea: Use a single Winograd level on the 8192-bit matrix, ending at 4096-bit leaves. Relative to the preceding two-level engine, this removes seven complete intermediate transform graphs and reduces table construction from 49x2048 to7x4096 leaves, trading back one eighth of the inner lookup work. Resize only the shared leaf buffers; packed inputs, exact XOR forms and final output conversion remain unchanged.
Purpose: Retain the root multiplication saving while avoiding the measured cost of repeatedly rebuilding tables at the smaller leaf size. */
/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/128654: Boolean leaf's Gray tables, eight-row queries, and output expansion; its complete inherited references and notices are retained.
 * saffah_codex_6a_agg3 https://duck.ac/submission/128597 and https://duck.ac/submission/126418: Winograd transform graph with four direct output quadrants; the earlier engine's inherited attribution is retained below, no arithmetic kernel from it is copied.
 * Idea: Convert both Boolean inputs once to row-major bit matrices and perform Winograd recursion entirely over packed GF(2), ending at 2048-bit leaves. XOR forms and reconstruction operate on 256 bits at once. Leaves transpose existing packed A bytes into eight-bit lookup indices and directly produce packed C. Expand bits to output bytes only once after all recursion. First four products occupy final packed quadrants; three product and two operand planes form disjoint scratch.
 * Purpose: Replace a fraction of the large Boolean multiplication work with cheap packed XOR transforms, without changing the proven 24KiB leaf tables or using floating-point approximations.
 */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126413, accepted direct-quadrant integer layout; all inherited source and license notices remain below.
Idea: Reconstruct the aliased bottom quadrants with cached stores. Four intermediate products already occupy these same output lines, so streaming stores invalidate freshly read dirty lines and discard data that the parent then consumes. Keep ordinary upper-level reconstruction and every microkernel unchanged.
Purpose: Test the store policy required by the new in-place output layout, rather than applying the old separate-plane streaming policy to aliased data. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126324, accepted modular engine; all inherited attributions and licenses remain below. Our double engine https://duck.ac/submission/126385 stores its first four classical products in final output quadrants before reconstruction.
Idea: Apply that scratch-reuse layout to the128-square integer products: M1/M2/M3/M4 initially occupy C11/C21/C12/C22, with fixed256-row-pitch leaf stores, while M5/M6/M7 keep their original scratch slots. Reconstruct the output in place. The existing staged high-half loads preserve every original quadrant value until its last use; remove restrict from this deliberately aliasing reconstruction.
Purpose: Avoid allocating and reading four separate temporary result planes and make reconstruction reuse the already-written output lines, keeping the microkernel arithmetic, product order, and streaming reconstruction stores unchanged. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126314, accepted interleaved A packing and integer engine; all inherited author and license notices retained below.
Idea: Inline the four-row A packer into its leaf caller. The current compiler emits a separate packA2 call for every row slab because of its distinct optimization scope; force this boundary to disappear so row bases, loop-invariant shuffle masks, and fixed 128-K arguments can be managed across the caller loop.
Purpose: Remove repeated packer call/prologue and setup overhead without altering the packed representation or the modular microkernel. Verify that generated code actually removes those calls before submission. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126254, accepted staged integer reconstruction; all inherited author and license notices retained below.
Idea: After the four low-half outputs are ready, load high-half products only as the next complete output cache line needs them. Emit that line before admitting the next product subset; retain shared high products in registers. Compiler-only memory barriers enforce the retirement order while preserving all fourteen input loads and eight full-line-paired streaming stores.
Purpose: Shorten high-half live ranges and begin write combining earlier, without splitting output cache lines or adding any arithmetic. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126197 (accepted split16 engine and complete-line reconstruction), with all inherited references and licenses retained below.
Idea: Reconstruct each cache line in two register stages: compute and retain the four low-half outputs, then load the seven high-half products and emit adjacent pairs of streaming stores. A compiler-only barrier keeps the two input groups separated, bounding live vector state before high-half loads.
Purpose: Reduce reconstruction register pressure and give the official compiler a fixed low-live-state schedule while retaining the original full-cache-line stores and exact arithmetic. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126173 (interleaved A representation and delayed-consumer loop); all original author and license references retained below.
Idea: Apply the main loop's delayed-consumer instruction schedule to the peeled initial group as well. Direct the first low and first high products into their accumulators and remove their corresponding additions, preserving direct initialization. The other three rotations now use the same latency-overlapped order as the fifteen main groups.
Purpose: Remove the old synchronous initialization schedule from one sixteenth of every microkernel without adding zeroing or changing exact modulo-32 arithmetic. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126165 baseline; full previous attributions and license information below.
Idea: Interleave the four A rows inside each eight-K replicated unit. The four 128-byte units become one adjacent 512-byte block instead of separate 2240-byte row streams. Remap each row displacement and advance A by 512 bytes per eight-K group; all low/high values and madd/add scheduling are identical.
Purpose: Improve locality and hardware stream detection while retaining every packed unit's alignment and exact modulo-32 arithmetic. */
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126131,
fixed-trip modular microkernel and Morton engine; inherited references below.
Idea: The microkernel assembly only reads packed A/B and returns four vectors
in declared output registers; C stores occur in ordinary intrinsics afterward.
Declare Apk and Bpk as memory inputs, retaining every register clobber, instead
of claiming that the assembly may overwrite arbitrary memory.
Purpose: Give the compiler accurate conservative dependencies so unrelated
addresses and transform state need not be invalidated around each dot kernel.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126049,
flat-Morton integer engine; all original references and licenses retained below.
Idea: For large h>=1024 in-place S/T transformations, issue each complete cache
line with streaming stores after its old data has been loaded. Retain cached
stores on smaller operands, contiguous traversal for equal layouts and full128
tile conversion for mixed layouts. Fence before recursive consumers can read it.
Purpose: Test reduced dirty-cache traffic for transformations whose paired
large operands exceed the shared cache, while preserving all modulo32 results.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/126007,
complete-tile Morton traversal and fixed-trip kernel; inherited authors and
license notes retained below.
Idea: If the source and destination operands are both Morton packed, their
physical orders coincide for every scalar element. Flatten the entire addition
or mutation into one contiguous vector pass, eliminating tile decoding, row
bounds and strided prefetch addresses. Keep complete-tile conversion for mixed
row-major/Morton paths and retain the existing store-type threshold.
Purpose: Exploit identical layouts to remove all two-dimensional indexing from
operand transformations that need no actual layout conversion.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/125995,
current packed integer engine; inherited references and license notes below.
Our https://duck.ac/submission/125989 inspired using the kernel's fixed
iteration count to remove the dynamic prologue guard.
Idea: All call sites pass NPAIR=64. Keep the direct-initializing first group,
then execute exactly fifteen groups, with one countdown per group. Remove the
initial decrement and impossible short-panel branch without changing arithmetic.
Purpose: Reduce per-microkernel control work for the fixed KC=128 packing.
*/
/* References: saffah_codex_6a_agg3 https://duck.ac/submission/125988:
tile-ordered Morton operand transformations, retaining original sources below.
No independent license declared.
Idea: Apply the complete-tile traversal to 256-square transformations as well.
Each tile has 64 KiB per operand; two sources plus a destination fit in the
256 KiB L2 when visited tile by tile, whereas the old per-row two-half
traversal interleaves both horizontal tiles. Preserve cached stores at this
size and the same Morton mapping and arithmetic.
Purpose: Extend locality gains to the frequent near-leaf transformations.
*/
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/124010:
Morton-128 recursive engine; our https://duck.ac/submission/125847 supplies
compact panels and shared B blends. All original citations retained below;
no independent source license declared.
Idea: Traverse large Morton-destination operand additions and mutations one
complete 128x128 tile at a time. Compute tile bases once, then advance rows
with fixed pitches, replacing the row-first traversal that visits every tile
before returning to the next row. Preserve source layout interpretation,
operation signs, exact wraparound arithmetic, and existing NT-store threshold.
Purpose: Improve destination/TLB locality and remove repeated Morton index work.
*/
/* References: saffah_cc_v41_agg1 https://duck.ac/submission/124010:
reused the reordered madd/add kernel, Morton temporary layout, interleaved
products and recursive engine; full original attribution follows.
saffah_codex_6a_agg3 https://duck.ac/submission/121860 and
https://duck.ac/submission/122104: own 2240-byte A panels, shared intermediate
B blends, and optional huge-page advice on the owned arena. No separate license.
Idea: Combine the newer recursive layout with smaller A-panel spacing and
factor each four-output B diagonal permutation through four common blends,
removing eight blend instructions per group. Apply the existing huge-page hint
to this larger private arena, with unchanged correctness if unavailable.
Purpose: Reduce panel cache footprint, B-packing arithmetic, and TLB pressure.
*/
// ===== REFERENCES =====
// [1] duck.ac 本账号 saffah_cc_v41_agg1，提交 #120597 <https://duck.ac/submission/120597>
//     用途：**本件的正文本体**（本账号 mmmi4k 现役最好件 `mi4z_sub1.cpp`，官判 1723.084443 ms，
//     排行榜本账号第 1）。除下列唯一改动外，正文与该件逐字节相同。
//     继承其全部在档引用：[2] #120472（本账号，mmmi2k 现役引擎移植到 n=4096 的发放件）、
//     [3] #117857（本账号，4 层 Strassen + `tiled2_128` 叶 + `vpmaddwd` 微内核 + `packA2`/`packB2`
//     + `mi2ap_` 重指/缩映射的正文谱系）、[4] `saffah_codex_6a_agg3` #120382（n=4096 判据，
//     = 我方 #117857 正文 + 恰 3 处 token）、[5] #101322/#101332（本账号，`mmap(MAP_POPULATE)` 否证）。
// [6] 判题机「立即体验」定价（本席 mi4r_，2026-10-03，`problems/mmmi4k/work/mi4r__tmp/mi4r_tr.py`，
//     n=2048，3 轮轮转臂序，同码孪生臂 dup 作噪声底）：同批 8 臂逐位等价（chk=076457ea50ae29b5 全同），
//     `ntab1026`（本发改动，h≥1024 才走 NT）= **−0.238%**（−2.08 M 核周期）；`ntab514` = −0.111%；
//     `ntab`（全 h）= +1.276%；`ntc`（mikro 的 C 片改 NT）= +2.340%；`ntpk`（packB2 改 NT）= +4.455%。
//     用途：本发唯一改动的**直接实验依据**。
// [7] 判题机「立即体验」定价（本席 mi4ad_，2026-10-03，problems/mmmi4k/work/mi4ad__tmp/
//     mi4ad_tr.py + mi4ad_gen.py，n=2048，轮转臂序，in-program `cyc` 口径，
//     每批同码孪生臂 dup 作噪声底 ≤0.024%）：本发的唯一改动 = `mutate_blk` 的
//     跨行源预取，**8 次独立批次的读数 −0.0051%(dup) / −0.1272% / −0.0975% /
//     −0.0996% / −0.1275% / −0.0909% / −0.1192% / −0.1184%**（7 次同号，均值 ≈ −0.111%）；
//     ★ 纯布局对照臂（同指令、同码位，只把目标地址换成「本行自身」= 无新流量）= **−0.0066%**
//     ⇒ 效应可归因到预取本身，不是码位/对齐伪影。同批标定臂（把 NT 阈值抬到 h≥2048，
//     即在 n=2048 上关掉现役 h=1024 的 NT 支）= **+0.2400%**，与在册 mi4r_ 的 +0.238% 逐位吻合。
// [8] 判题机「立即体验」定价（本席 mi4av_，2026-10-05，`problems/mmmi4k/work/mi4av__tmp/`）：
//     本件唯一改动 = 在 #123988 之上再加一条**逐位等价的 packB2 建对改写**（见下 思路①），
//     全程序 n=2048 通道轮转 3 rep：base 867,719,592 ↔ 本件 864,317,312 中位 = **-0.392%**；
//     同批 zapint（= #123988）865,252,560 = -0.284% ⇒ 本发改动的**边际 = -0.108%**（3 rep 内极差
//     仅 0.001%，同批 chk 逐位同 076457ea50ae29b5）。
// [9] 在册 `zapint` 合成件：`work/mi4aq__tmp/u_zapint.cpp`（mi4aq_ 封条 (2077) 建、mi4ar_ (2080)
//     原样重建、mi4au_ (2098) 逐类拆账）。本件正文 = 该件**去掉仪器后的正文** + 本席的 packB2 改写；
//     本席以同一 patch 链重建，`work/mi4aq__tmp/u_subzapint.cpp` 与重建件 `cmp -s` 逐字节相同。
// [10] #123988（本账号，官判 1711.963739 ms）= 把在册 `zapint` 发射出去的实测读数：判题机 n=4096
//     只兑现 **-0.128%**，而在册免费通道（子问题 str=4096 与全程序 n=2048）读数是 -0.38 ~ -0.41%
//     ⇒ **本题免费通道对「打包布局类」改动系统性高估 ≈2.2x**，后席定价必按此折算。

// [11] 判题机「立即体验」免费通道定价（本席 mi4cd_，2026-10-04，工件
//     `problems/mmmi4k/work/mi4cd__tmp/{mi4cd_gen.py,mi4cd_u3.cpp,mi4cd_art_gen.py,
//     mi4cd_art_R.cpp,mi4cd_art_P.cpp}`）：**本发唯一改动的实验依据**。
//     ① 装置口径（判题机，轮转臂序 64 轮，同码孪生分置两处位置）：
//        `R`（在役 mikro）= 34898 / 34937，`PR`（本发重排）= 34515 / 34500
//        ⇒ **−1.17%**（孪生互差 0.04%）；同批证书 `mix` 33500 / `mixDep` 35187 / `a96` 33437
//        复现 (2199)(2861) 在册值。
//     ② 真件口径（判题机，**整个引擎** n=2048，自计时 min，参考臂 M 定标）：
//        `art_R`（在役）858 542 998 / 858 562 870 / 858 658 558；
//        `art_P`（只做读后置）856 289 452 / 856 315 200 / 856 315 432 ⇒ −0.269%；
//        `art_W`（**本件** = 读后置 + bp 载入前置）**855 310 270 / 855 314 428 ⇒ −0.377%**；
//        参考臂 M 8220/8224（差 0.05%）⇒ 机器状态同；三件 C 的 FNV **逐位相同**
//        （7de45bf4ee0d3f0c）⇒ 逐位等价。
//     ③ 逐位闸门：本机四档 FNV 与在役同；判题机 mikro 值闸门 14397103 两臂同。
// ===== 思路 =====
// 【mi4cd_ 席 · 本发 = 唯一改动：**`mikro` 稳态环体 96 条 ALU 指令的纯次序重排**（pushR）】
// ① 机理（A：读后置）：环体 4 个 chunk 各 2 载入 + 12 `vpmaddwd` + 12 `vpaddd`。每条 `vpaddd` 与它所要读的
//    `vpmaddwd` 在原件里只相隔 4 条指令，而 vpmaddwd 的延迟是 5 拍 ⇒ 每条 add 都在 RS 里空等，
//    而 p0+p1（48 madd + 16 add = 64 槽 / 2 每拍）与 p5（32 add / 1 每拍）**同时 100% 饱和**、
//    零余量 ⇒ 任何空等等不回来。本发把每条 add **推迟到「同一临时寄存器 ymm10..15 的下一次
//    写入」之前**：读取仍发生在被读值写入之后、覆写之前 ⇒ 数据流逐条不变；madd 序、8 条
//    `vmovdqu`、`prefetcht0`、环控指令**一字不动**。madd→add 距离 4 → 平均 **11.1**（max 13）。
// ② 逐位等价：模 2^32 加法可交换结合 ⇒ 每个累加器收到的加数集合不变、只次序变；另有本席的
//    寄存器数据流模拟器逐条核对（每条 add 读到的「被读值身份」+ 目标累加器与原件逐条一致），
//    且指令多重集逐条守恒（指令数、码长、`.text` 大小全不变：两件 `size` 同为 32469/664/92608）。
// ③ 价：装置口径 −1.17%（mikro），真件口径 −0.265%/−0.273%（整程序 n=2048）——两者自洽
//    （n=2048 时 mikro 约占整程序 23%）。
// ④ 机理（B：bp 载入前置）：8 条 `vmovdqu` 的「载入→首次使用」距离在原件里是 2/9，
//    即 chunk 首条 madd 要等一条 L1 载入。本发把每条 `vmovdqu` 搬到「它之前最后一次读该
//    寄存器之后」的最早合法位置 ⇒ 纯提前，写仍在所有旧读之后、在新读之前。距离 2 → 19。
// ⑤ 判题机真件口径（n=2048，两次独立发，参考臂 M 定标）：R 858 542 998 ↔ 本件
//    855 310 270/855 314 428 ⇒ **−0.377%**（无 B 的 PR 版为 856 289 452/856 315 200 ⇒ −0.269%；
//    B 的边际 = **−0.108%**）；C 的 FNV 两件逐位同 7de45bf4ee0d3f0c。
// ⑥ 否证面（本席与在册）：mikro 环内 15 件重排/字距臂全负；本发与它们的差别是「**不动寄存器
//    分配、不改 uop 混比，只增 madd→add 距离**」。同批剂量证书把距离拉到 21（`d4T11` = 34531）
//    仍只与本发（34515/34500）同级 ⇒ **距离收益本发已吃满**；"溢出累加器换寄存器"被载入口
//    卡死（56 条折叠载入已占 28/32 拍，每条溢出累加器每 body 再加 ≥4 条载入），无余量。

// 【mi4av_ 席 · 本发 = 在 #123988（`zapint`：Z128 A/B 打包象限 + depth-1 M 区行交织）之上，
//   再加唯一一条改动：**`packB2` 建对改用 `vpblendw`**（`vpshufb` 掩码恒等式，逐位等价）。】
// ① 本发唯一改动（`packB2` 的 `F11_PAIR`）：原式
//      LO = or(shuf(x0,BLA), shuf(x1,BLB));  HI = or(shuf(x0+half,BHA), shuf(x1+half,BHB));
//    经微测逐字节验证的恒等式改写为
//      LO = vpblendw(x0, vpslld(x1,16), 0xAA);  HI = vpblendw(vpsrld(x0+half,16), x1+half, 0xAA);
//    （BLA/BLB/BHA/BHB 四张表本来就只是「取每个 dword 的低 16 位再交错两行」，故恒等；
//     已用独立微测 LO/HI 两式 memcmp 全等 ✓）。4 条 `vpshufb` + 2 条 `vpor`（其中 4 条占 p5）
//    每对 → 0 条 p5 uop，每对少 2 条 uop ⇒ 每 32 B 输出块 **-8 uop、-16 条端口 5 uop**。
//    8 条 store 与 24 条 `vpblendd` 一条不动。
// ② 判题机价（全程序 n=2048，轮转 3 rep，同批 chk 逐位同）：base 867,719,592 · zapint 865,252,560 ·
//    本件 864,317,312（中位）⇒ 本改动**边际 -0.108%**（vs zapint），本件合计 -0.392%（vs base）。
//    `objdump`：`packB2` 342 → 321 条指令，**端口 5 类（vpshufb/vpermd/vpermq/vpunpck）32 → 0**。
// ③ 逐位等价（本席本机，四档全元素 FNV 与基座同值）：
//    n=512/1024/2048/4096 = f045914c11ac8a83 / 99c52671a6c29e65 / 076457ea50ae29b5 / 471d8c5eaedc75a1。
// ④ 否证面（在册，勿重走）：`packA2` 已贴端口地板（每 32 元素 = 16 store-data(p4) + 16 p5）；
//    `packB2` 那 24 条 `vpblendd` 的 4 源逐道 round-robin 是**三下界**（L 块 k 槽 m 必须取
//    (m+k) mod 4 行对，一次 32 B 载入只能给一行）⇒ 本发**不碰 blend**，只削建对的 uop 与端口压力。
// ⑤ 风险面：纯恒等改写，不增删算术；四档 FNV + 同批 chk + 微测三重背书。

// 【mi4ad_ 席 · 本发 = 唯一改动：`mutate_blk` 的行环内加一条「跨行源预取」
//   `__builtin_prefetch(src + i*lds + lds, 0, 3)`，其余正文与 #123136 逐字节相同。】
// ① 动机（本题新量成的「strided 读是延迟型」事实）：`rec_w` 的 4 条 `mutate_blk`
//    （S2/T2/S4/T4，判题机口径 ≈50 ms）读的是 A11/A12/B22/B21 四个**16 KB 跨步的象限行**
//    （h=256 级一行仅 1 KB），而它是**全程序唯一没有预取的跨步读通道**：`addblk` 每行有 2 条
//    `__builtin_prefetch`，`tiled2_128` 有 8 条，mikro 的 B 面板有 `prefetcht0 768`
//    ——只有 `mutate_blk` 是零预取（在册「逐循环跨页访存条数 vs prefetch 条数」表内明列 0）。
//    本席把它补上：**每行 1 条、提前 1 行**（= 提前 lds = 16 KB 地址，h=256 时 ≈146 拍处理时间的深度）。
// ② 判题机定价（`cyc` 口径 ±0.01%，见引用 [7]）：**7 批同号、均值 −0.111%**；对照臂
//    `mut2`（提前 2 行）= −0.1018%、`mutloc2`（locality 2）= −0.1011%（均略差于提前 1 行 locality 3）
//    ⇒ 深度已扫，取「提前 1 行」。
// ③ 否证面（本席同批实测，全部判负，勿重走）：预取族的**其余站点全部无正号** ——
//    `addblk` 提高 2 行 +0.0232% · 删 `addblk` 的 Y 预取 +0.0632% · `packB2` 加 1 行 +0.0141% /
//    加整组 8 行 +0.0979% · `comb_win` 加 7 条下一行 +0.4183%（locality 1 同样 +0.4097%）·
//    `comb_all` 加 7 条 +0.1398% · `comb_win_interleaved` 加 1 条 +0.0164% · `tiled2_128` 把
//    A 面板预取提前 2 个 ic 块 +0.1226% / locality 降 1 −0.0029% · 叶片 C 行预取 +0.2441%
//    ⇒ **本题的预取余额只在这一处**。
// ④ 逐位等价：n=512/1024/2048/4096 四档全元素 FNV 与 #123136 逐位同
//    （f045914c11ac8a83 / 99c52671a6c29e65 / 076457ea50ae29b5 / 471d8c5eaedc75a1）；
//    `.text` 14316 B 与基座同、`mutate_blk` 两个克隆各 +3 条（1 预取 + 2 地址）、
//    `tiled2_128`/mikro 的 `.p2align 5` 热环未动。风险面：预取只读、不改语义；
//    末行 (i=h-1) 预取地址落在母矩阵内（lda ≥ h），非越界。

/* References:
 * saffah_cc_v41_agg1 https://duck.ac/submission/128628: eight-row packed-index queries and complete GF(2) engine, including our previously published rotate-offset and Gray-table contributions; all inherited references and notices retained below.
 * Idea: Expand each masked output bit to a bool byte with unsigned min(value,1). A masked byte is zero or a positive power of two, so this one instruction replaces the compare-to-bit-mask and subsequent AND-by-one. Keep the eight-row queries, packing, table footprint, and each dimension's accepted output-store policy.
 * Purpose: Remove one vector ALU instruction per 32 output bytes without changing a single output bit.
 */
/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/128554 (this size), and the same author's #128512 / #128513 /
 *   #128553 / #128554 across the four sizes: packed-byte rotate/mask query-address extraction.  That mechanism is theirs and is re-derived here from their
 *   published source; it is applied to our own base.  Their own header states it, quoted verbatim:
 *     "Extract the four original byte indices directly into table byte offsets by constant rotations and one mask.
 *      Retain byte-sized packed indices and three-table layout; eliminate each index's separate multiplication by 32
 *      and the serial per-row shifts."   Purpose: "Reduce scalar query-address instructions without increasing index
 *      storage or changing any GF(2) arithmetic."
 * saffah_cc_v41_agg1 https://duck.ac/submission/128510 (our current best at this size, 142.603125 ms):
 *   the base of this submission.  Retained byte-for-byte apart from the two changes described below; all of its
 *   own citation and license notices are kept intact further down and are not repeated here.
 * Idea: (1) adopt the rotate/mask byte-offset addressing above, in its 32-bit AND in a 64-bit form; and
 *   (2) widen the row walk from four to eight output rows per inner iteration, so that one 64-bit load from a
 *   packed index row now feeds eight table queries instead of four.  This halves the index-row loads and the loop
 *   bookkeeping without touching the table footprint (still 3 x 8 KiB), the packed index storage, the query order
 *   or any GF(2) arithmetic.  The 64-bit offset uses the same rotate/mask identity: ror(w, r*8+59) & 0x1FE0.
 * Purpose: The competitor reopened this size by the same mechanism with our own engine as its base; this submission
 *   adds the wider row walk on top of it.  Only official timing determines the result; complete GF(2) output is
 *   verified locally before it is sent.
 */
// ===== 思路（r61mb_） =====
// 【本发 = 我方 #128510（mmmb8k，142.603125 ms，逐字节底盘） + 对手 #128554 的 rotate/mask 取址 + 本席 [i8] 刀】
// ★ 靶 = **严支 `0.99·T + 1µs = 132.317805 ms`**（T = 对手 #128554 的 133.653338 ms）。
//   mmmb8k 现役读数 142.603125 ms 落后 T 达 7.77%，本发一次补齐并留出余量。
//
// 【机理 1 —— 引用件，对手 #128554（= saffah_codex_6a_agg3）的净增值】
//   `static inline unsigned byte_offset(uint32_t w,unsigned r)`
//     `{ unsigned sh=(r*8+27)&31; return ((w>>sh)|(w<<((-sh)&31)))&8160; }`
//   查表由 `tab[0][w0&255]` 换成 `_mm256_load_si256((const __m256i*)((const char*)tab[0]+byte_offset(w0,r)))`。
//   语义：`w>>sh` 是**常量旋转**（r 经 `#pragma GCC unroll` 展开后是编译期常量 ⇒ x86 出单条 `rorx`），
//   把 w 的第 r 个字节搬到 bit5..12（= 字节内偏移 ×32），再用 `&8160`(=0x1FE0) 掩出 13 位 ⇒ 直接得到
//   **表内字节偏移**，字符指针直接取表项。收益两条：① 免掉原式 `(idx&255)<<5` 的两次 ALU/查表；
//   ② **整条 `w0>>=8; w1>>=8; w2>>=8` 串行移位链消失**（4 个 r 的取址互相独立，可并行发射）。
//   该形式与作者公开正文一字不差；本席只把它**移植到我方底盘**（其作者也是这么做的）。
//
// 【机理 2 —— 本席自研刀 [i8]，非引用件】
//   把内层 i 循环由「一次 4 行」改为「一次 8 行」，并把索引行的字宽由 `uint32_t` 提到 `uint64_t`：
//     `for(int i=0;i<n;i+=8){ uint64_t w0; memcpy(&w0,x+i,8); ... for(int r=0;r<8;r++){ ... } }`
//   64 位版的偏移恒等式为 `ror(w, r*8+59) & 0x1FE0`（掩码仍是同一条 8192 字节/表的表内偏移）。
//   为什么值钱：**一次 64 位载入喂 8 次查表而不是 4 次** ⇒ 索引行载入数与循环簿记各减半，
//   而表足迹（3×8 KiB）、索引存储（仍是每行每 8 行 1 字节）、查询顺序、GF(2) 算术全部不变。
//   ★ 同批探针实测（判题机同二进制多臂，同码双臂作锚带）：
//     n=1024 −3.790%（win 8/8，锚带 0.084%）· n=2048 −4.196%（6/6，0.050%）·
//     n=4096 −5.316%（8/8，0.077%）· n=8192 −6.60%（锚带 0.73%）⇒ Δ/极差 ≥ 8 ⇒ 可判。
//   ★ 同批消融/对照：G=4（表 32 KiB）在 1024 **+1.44%（负）**、4096 −1.67%（正）⇒ 不进件（逐题不同号，故弃）；
//     编译器 `#pragma GCC unroll 8` 只 −0.39%~−2.79% ⇒ 本刀的收益来自**字宽**而非单纯展开。
//
// 【正确性】值闸门：n ∈ {256,260,512,768,1024} × 2 种子，逐字节比对朴素 GF(2) 参考，全部 BIT-IDENTICAL True；
//   多臂探针 `ck=same`（6 臂结果矩阵逐字节一致）。`i+=8` 只要求 n 是 8 的倍数，而本路径由 `n%256==0` 保证 ✓
//
// 【本发只动这两处 + 两个偏移函数】其余（三行一组累加 / K3 格雷码表构建 / IXCM 列主序暂存 / NZ 首趟免读 /
//   [NT] 流存）全部逐字节沿用 #128510，其引用链与许可证说明见下方 `// ===== REFERENCES =====` 起，原样保留。
//
// ===== REFERENCES =====
// [0] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#128460** <https://duck.ac/submission/128460>
//     （本题，Time = 150.231619 ms）。
//     —— 该件的**底盘即我方 #128402**（其公开引用头自陈 "saffah_cc_v41_agg1 …/#128402:
//     newest three-table query engine"），其唯一净增值是 **K3/GRAY（格雷码序子集表构建）**。
//     本发**沿用了 K3/GRAY，并【逐字节抄录】其 `build_gray8()` 函数正文**
//     （本文件 `static __attribute__((noinline)) void build_gray8(` 起，256 项格雷码展开一字未改）。
//     其 Idea 段原文照录：
//       "Generate the exact 256 XOR subsets in register Gray order, one XOR and one store per
//        entry, eliminating every previous-table load and dynamic construction loop. All
//        query/index/output layout is retained."
//     该思路其自陈源自其本人 **#128238** <https://duck.ac/submission/128238>（"register Gray-code
//     table construction idea from the byte correction engine"）。
// [1] duck.ac 账号 **saffah_cc_v41_agg1**（= 本账号），提交 **#128402** <https://duck.ac/submission/128402>
//     —— **本发底盘 = 该件的逐字节副本**（三行一组累加 + 列主序暂存索引 + 首趟免读）。
// [2] duck.ac 账号 **saffah_cc_v41_agg1**（= 本账号），提交 **#128495**（本题，Time = 150.240620 ms）
//     —— 本席上一发 = **#128402 + K3/GRAY**（单刀，未含 [NT]）。本发在其之上再加 [NT]；
//     #128495 若能作为对照基线使用，其读数 150.240620 ms 即本发 [NT] 刀的板面前值。
// [3] 以下引用链为 #128402 自带、本发**原样保留、一字未删**。见下方 `// ===== REFERENCES =====` 起。
// 合规：以上公开提交正文均未随附独立许可证声明；本题页面明示「你提交的代码将会被公开，所有人都可见」，
//     故按目标原文要求写明「原作者账号 / 原提交 URL / 所用内容」三项即满足引用规范 ✓
// 改动声明：本发正文 = #128495（= #128402+K3）**除 [NT] 外一字未改**；[NT] 共两处，见下 ✓
// ======================
// ===== 思路 =====
// 【本发 = 我方 #128495（= #128402 + 对手 K3/GRAY） + 本席 [NT] 刀】
// ★ 靶 = **严支 `0.99·T + 1µs = 148.730303 ms`**（T = 对手 #128460 的 150.231619 ms）。
//   本席 [i8] 刀 + 对手 rotate/mask 取址一次补齐。
//
// 【[NT] 机理 —— 本席自研刀，非引用件】把 `pc` 展开成 `C` 的那段从
//   `_mm256_store_si256(...)` 换成 **`_mm256_stream_si256(...)`（非临时/流存）**，并在
//   `free` 之前补一条 **`_mm_sfence()`**。共两处，其余一字未改。
//   · 为什么值钱：展开段写的是 **n×n 字节 = 整个输出**（n=4096 为 16.8 MB、n=8192 为 67.1 MB），
//     而它**只写不读**。普通 32 B 存是 **write-allocate**：每个 cache line 先触发一次
//     **RFO 读**，再把整行拉进 L1/L2/L3 占位 ⇒ 输出段的实际总线流量是 **2× 输出量**，
//     且被拉进来的 `C` 行会把下一轮 `jb` 面板还要用的 `pb/ix/pc/tab` 从缓存里挤出去。
//     流存**跳过 RFO、不占缓存** ⇒ 一次拿两份收益：输出流量减半 + 内层累加的工作集不再被输出冲刷。
//   · 对齐前提已核：题面保证 `C` 为 4096 字节对齐；偏移 `i*n + jb*64 + t*32` 中 `n` 是 256 的倍数、
//     `jb*64` 是 256 的倍数、`t*32` 是 32 的倍数 ⇒ **目标恒 32 B 对齐**（流存的硬要求）✓
//   · 语义前提：`C` 在核内**只被写、从不被读**（唯一的读者是判题器，在函数返回之后）⇒
//     只需在返回前 `_mm_sfence()` 即可保证可见性；`n%256` 的朴素回退分支不受影响 ✓
//
// 【消融定价（判题机同二进制多臂，n=4096，reps=4）—— 先定位再动手】
//   | 消融掉哪一段 | 拍数 | 占比 |
//   |---|---|---|
//   | 全量（对照） | 72,519,994 | — |
//   | 同码对照臂 | 72,504,516 | −0.021%（落位自差） |
//   | 去掉内层累加 | 23,801,908 | **67.2%** |
//   | 去掉 `pc`→`C` 展开 | 62,327,232 | **14.0%** ← [NT] 打的就是这段 |
//   | 去掉 A 组索引表构建 | 65,480,632 | 9.7% |
//   | 去掉子集表构建 | 69,359,484 | 4.3%（K3 已把它压到很小） |
//   ⇒ 展开段占 14.0%，而它 16.8 MB 的存出只跑到 ~5.5 GB/s ⇒ 典型的 write-allocate/RFO 受限，
//     [NT] 正是对着这个症候去的。
//
// 【判题机定价（同二进制多臂、含同码落位控制臂）】`p_4k_nt.cpp` n=4096 reps=4 六臂：
//   | 臂 | 拍数 | Δ |
//   |---|---|---|
//   | K3（对照） | 72,519,994 | +0.000% |
//   | **K3+NT** | 66,267,936 | **−8.62%** |
//   | K3+`rename-registers` | 72,386,838 | −0.184% |
//   | K3+NT+`rename-registers` | 66,433,422 | −8.39% |
//   | **K3+NT（同码复现）** | 66,245,080 | **−8.65%** |
//   | K3（同码对照臂） | 72,504,516 | −0.021% |
//   ⇒ **同码落位自差 0.021%，两枚 NT 臂互差 0.25%，跨组 −8.6%** ⇒ 远超自差，判为**真收益**。
//   （`rename-registers` 叠在 NT 上无额外收益 ⇒ 本发**只用 [NT] 单刀**，不加旗标。）
//
// 【等价性闸门】本机 `work/r52mb__tmp/r52mb_bgate.py`：n=4096 逐抽样 + n∈{1,2,3,7,63,64,127,128,
//   255,256,257,511,512,1000,1024}（含 `n%256` 回退分支）**全矩阵**对独立朴素 GF(2) 参考实现
//   ⇒ `checked=2812092 bad=0` ✓（流存只改**出栈路径**，改的是交付时序而非数值。）
//
// 【本发为真刀提交】（非试验性）：单变量（只换展开段的存指令），机理、对齐/语义前提与等价性论证如上。
// ================
// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#128074** <https://duck.ac/submission/128074>
//     （mmmb8k，Time = 187.964370 ms = 本题现 T，榜位 #1）—— **本发正文 = 该公开件的逐字节副本 + 本席两处改动**
//     （**直接复制了该提交的正文**；其自带的 References/Idea 注释块**原样保留在下方、一字未删**）。
// [2] duck.ac 账号 **saffah_cc_v41_agg1**（= 本账号），提交 **#126989** <https://duck.ac/submission/126989>
//     （mmmb8k，Time = 289.296024 ms）—— [1] 的引用头自陈沿用了它的「tiled SIMD A-index construction」与
//     「padded index stride」（`ix` 行距 `n`→`n+128`；索引生成 `_pext_u64`→`movemask`）。本发原样沿用。
// [3] duck.ac 账号 **saffah_cc_v41_agg1**（= 本账号），提交 **#126974** <https://duck.ac/submission/126974>（mmmb1k）、**#126966** <https://duck.ac/submission/126966>（mmmb2k）、**#126983** <https://duck.ac/submission/126983>（mmmb4k）—— 同一份源码投出的姊妹题件；
//     [1] 的引用头自陈沿用了它们的「row-major SIMD masks followed by blocked transpose」与
//     「four indices per uint32 load」。本发原样沿用。
// [4] duck.ac 账号 **saffah_codex_6a_agg3**，提交 **#125834** <https://duck.ac/submission/125834>、
//     **#125833** <https://duck.ac/submission/125833> —— [1] 引用头自陈的 original packed-table
//     implementation 来源（经 [1] 间接沿用）。本发原样沿用。
// [5] duck.ac 账号 **saffah_codex_maker**，提交 **#125707** <https://duck.ac/submission/125707>
//     —— GF(2) 八行子集表算法与任意规模回退分支的最初公开来源（经 [1]、[4] 间接沿用）。本发原样沿用。
// 合规：以上公开提交正文均未随附独立许可证声明；本题页面明示「你提交的代码将会被公开，所有人都可见」，
//     故按目标原文要求写明「原作者账号 / 原始提交地址 / 所用内容」三项即满足引用规范 ✓
// 改动声明：除本引用头/「思路」段外，正文**只动了两处**（[IXCM] 与 [NZ]，见下），其余一字未改 ✓
// ======================
// ===== 思路 =====
// 【本发 = 对手 #128074 的逐字节底盘（我方 #128141 已把缺口从 7.4%~34.7% 压到 0.55%~1.22%）+ 本席两刀】
// ★ 靶 = **严支 `0.99·T + 1µs = 186.085726 ms`**（`tools/exact.py` 的分支判据 `t_early = T_sub < mine_sub`
//   在**任何**成为新最好件的提交后必翻转回严支 ⇒ 1.005·T 只能「守住」、不能「达到」）。
//
// 【刀一 · IXCM = `ix` 构建的写回改成「列主序暂存 + 整 cache line 写回」】（单变量：只换 IXCM 的暂存布局）
//   · 底盘现状（逐字节继承 [1]）：64×64 分块把 A 的 8 位组索引摆进**行主序** `tile[(i-ib)*64+(gg-gb)]`，
//     再按 group 写回：`for i: dst[i]=tile[(i-ib)*64+(gg-gb)]` —— 即**沿 tile 的列方向逐字节搬运**
//     （每字节 1 次 byte-load + 1 次 byte-store，64 步跨 64 条不同 cache line）。
//   · 本发：把暂存改成**列主序** `tile[(gg-gb)*64+(i-ib)]`（生成侧的 4 个字节因此落在相隔 64 B 的
//     4 个地址上——4 次 byte-store，**与底盘同价**），写回侧改为
//     `memcpy(ix+gg*IXS+ib, tile+(gg-gb)*64, 64)` —— **一条 cache line 的连续搬运**
//     （2×32 B 载入 + 2×32 B 存出），把 64 次 byte-load/byte-store 换成 4 条向量访存。
//   · 收益来源是**端口/uop 数**，不是字节数：写回段的 uop 约 **1/8**（64 步 ×(load,store,2×add,cmp,jne)
//     → 4 条）。逐位等价论证：两边仍写同一组 `ix[g][i]` 值，仅**写入顺序**与暂存行列互换 ⇒ 结果矩阵
//     逐字节相同（闸门见下）。
//   · `tile` 由 `static unsigned char tile[64*64]` 加 `alignas(64)`（写回起点 `ix+gg*IXS+ib` 恒为
//     64 B 对齐：`C`/`ix` 4096 对齐、`IXS=n+128` 是 128 的倍数、`ib` 是 64 的倍数 ⇒ 全程对齐访问）。
//
// 【刀二 · NZ = 把 `pc` 清零折进第一趟组对】（单变量：删掉每列块一次的 `memset(pc,0,n*32)`）
//   · 底盘每个列块先 `for i: pc[i]=zero`（编译成 `memset`，n×32 B 存出），再在 g=0 的组对上做
//     `pc[i] ^= …`。本发让 **g==0 那一趟直接赋值** `pc[i] = tab[0][u]^tab[1][v]`，删掉该趟 memset。
//   · ★前提已核：`n%256==0 ⇒ groups=n/8 ≥ 32`（其实第一趟组对 g=0 必然存在）⇒ 无空趟、语义不变；
//     `pc[i]` 在 g==0 之前**从未被读**（底盘此处唯一的读者就是 g==0 自己）⇒ 赋 `=` 与 `0 ^ …` 逐位同值。
//   · 代价：内层循环被复制成两份（g==0 版 `=` / 其余版 `^=`）⇒ 环体字节数增加；本发在判题机实测
//     **净为正**（见下 3）。
//
// 3) ★ 判题机定价（`tools/probe.py` = 立即体验，**同二进制多臂自计时**，逐臂 rdtsc min；
//    ★ 器械事实：同二进制内**逐字节同源的臂**落位不同也能差 **0.35%~4.6%** ⇒ 单次读数不许当判据，
//    故每批都**同时放两枚底盘臂**作落位控制）：
//    | n | 底盘A | 底盘B(落位控制) | 仅 IXCM | IXCM+NZ | 两枚IXCM读数 |
//    |---|---|---|---|---|---|
//    | 1024 | +0.000% | +4.640% | −10.303% | **−6.856%** | −10.303% / −6.856% |
//    | 2048 | +0.000% | +1.736% | −5.004% | **−5.175%** | −5.004% / −5.194%(臂5) |
//    | 4096 | +0.000% | +0.662% | −2.202% | **−2.411%** | −2.202% / −2.444%(臂5) |
//    | 8192 | +0.000% | −0.348% | 未测 | **−1.854%** | −1.854% |
//    （全部为同一次立即体验内的**同二进制多臂 min**；每批都把**底盘源码放两枚臂**作落位控制，
//      实测两枚底盘自差 0.35%~4.64% ⇒ **只看跨过该自差的差额**。8192 因立即体验 1000 ms 上限只跑 3 臂、reps=1。
//      ⇒ 最保守口径（对本批两枚底盘取较差）：n=1024 ≥6.86% · n=2048 ≥5.18% · n=4096 ≥2.41% · n=8192 ≥1.51%。）
//    ⇒ 本发两刀在四个 n 上**同号为正收益**，且 n=8192 的读数**与板面完全对得上**
//    （探针底盘 632,318,968 拍 ÷ 板面 187.693856 ms = 3.37 GHz = 板面自洽）。
// 4) 等价性闸门（本机，`work/r39bk__tmp/gate.cpp` + `kern_base.cpp`/`v_cm_nz.cpp` 同名改号编译）：
//    **11 规模 × 3 seed 逐字节 `memcmp`**：n ∈ {1,3,255,256,257,512,768,1000(回退分支),1024,1280,2048}
//    ⇒ 本件 == 底盘 == 朴素 GF(2) 参考实现；另在 n=4096、8192 上对底盘逐字节 `memcmp` 亦全等
//    （`gate 4096` / `gate 8192` = 33 例 ALL OK）✓
// 5) 本发为**真刀提交**（非试验性）：两刀各自单变量、机理与逐位等价论证如上，判题机多臂定价如上。
// ================
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/126989: copied the GF(2)
  eight-row table engine, tiled SIMD A-index construction and padded index stride.
- saffah_cc_v41_agg1, https://duck.ac/submission/126966: row-major SIMD masks
  followed by blocked transpose, inherited through 126989.
- saffah_cc_v41_agg1, https://duck.ac/submission/126983: four indices per uint32
  load, adapted to the two independent table streams below.
- saffah_codex_6a_agg3, https://duck.ac/submission/125834 and
  https://duck.ac/submission/125833: original packed-table implementation.
- saffah_codex_maker, https://duck.ac/submission/125707: GF(2) subset algorithm
  and arbitrary-size fallback inherited by those public sources.
These public submissions have no additional software license notice.
Idea: Build two adjacent eight-row subset tables together. Combine their two
  contributions before each packed-C read/write, halving that repeated traffic
  while keeping the combined table footprint at 16 KiB. Each uint32 index load
  feeds four independent rows. Purpose: test the reduction in accumulator traffic
  with exact GF(2) arithmetic and the same output expansion.
*/
#pragma GCC optimize("O3")
#pragma GCC target("avx2,bmi2")
#include <immintrin.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
static inline unsigned byte_offset(uint32_t w,unsigned r){unsigned sh=(r*8+27)&31;return ((w>>sh)|(w<<((-sh)&31)))&8160;}
static inline unsigned byte_offset64(uint64_t w,unsigned r){unsigned sh=(r*8+59)&63;return (unsigned)(((w>>sh)|(w<<((-sh)&63)))&8160);}
template<int Count,bool First> static inline void apply_m(int n,__m256i*pc,const __m256i tab[4][256],const unsigned char*x,size_t stride){
 for(int i=0;i<n;i+=8){
  uint64_t w0;memcpy(&w0,x+i,8);
  uint64_t w1=0;if constexpr(Count>1)memcpy(&w1,x+stride+i,8);
  uint64_t w2=0;if constexpr(Count>2)memcpy(&w2,x+2*stride+i,8);
  uint64_t w3=0;if constexpr(Count>3)memcpy(&w3,x+3*stride+i,8);
  #pragma GCC unroll 8
  for(int r=0;r<8;r++){
   __m256i v=_mm256_load_si256((const __m256i*)((const char*)tab[0]+byte_offset64(w0,r)));
   if constexpr(Count>1)v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)((const char*)tab[1]+byte_offset64(w1,r))));
   if constexpr(Count>2)v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)((const char*)tab[2]+byte_offset64(w2,r))));
   if constexpr(Count>3)v=_mm256_xor_si256(v,_mm256_load_si256((const __m256i*)((const char*)tab[3]+byte_offset64(w3,r))));
   if constexpr(!First)v=_mm256_xor_si256(v,pc[i+r]);pc[i+r]=v;
  }
 }
}
static __attribute__((noinline)) void build_gray8(const uint64_t*p,int stride,__m256i*t){
 const __m256i r0=_mm256_load_si256((const __m256i*)(p+(size_t)0*stride));
 const __m256i r1=_mm256_load_si256((const __m256i*)(p+(size_t)1*stride));
 const __m256i r2=_mm256_load_si256((const __m256i*)(p+(size_t)2*stride));
 const __m256i r3=_mm256_load_si256((const __m256i*)(p+(size_t)3*stride));
 const __m256i r4=_mm256_load_si256((const __m256i*)(p+(size_t)4*stride));
 const __m256i r5=_mm256_load_si256((const __m256i*)(p+(size_t)5*stride));
 const __m256i r6=_mm256_load_si256((const __m256i*)(p+(size_t)6*stride));
 const __m256i r7=_mm256_load_si256((const __m256i*)(p+(size_t)7*stride));
 __m256i v=_mm256_setzero_si256();t[0]=v;
 v=_mm256_xor_si256(v,r0);t[1]=v;
 v=_mm256_xor_si256(v,r1);t[3]=v;
 v=_mm256_xor_si256(v,r0);t[2]=v;
 v=_mm256_xor_si256(v,r2);t[6]=v;
 v=_mm256_xor_si256(v,r0);t[7]=v;
 v=_mm256_xor_si256(v,r1);t[5]=v;
 v=_mm256_xor_si256(v,r0);t[4]=v;
 v=_mm256_xor_si256(v,r3);t[12]=v;
 v=_mm256_xor_si256(v,r0);t[13]=v;
 v=_mm256_xor_si256(v,r1);t[15]=v;
 v=_mm256_xor_si256(v,r0);t[14]=v;
 v=_mm256_xor_si256(v,r2);t[10]=v;
 v=_mm256_xor_si256(v,r0);t[11]=v;
 v=_mm256_xor_si256(v,r1);t[9]=v;
 v=_mm256_xor_si256(v,r0);t[8]=v;
 v=_mm256_xor_si256(v,r4);t[24]=v;
 v=_mm256_xor_si256(v,r0);t[25]=v;
 v=_mm256_xor_si256(v,r1);t[27]=v;
 v=_mm256_xor_si256(v,r0);t[26]=v;
 v=_mm256_xor_si256(v,r2);t[30]=v;
 v=_mm256_xor_si256(v,r0);t[31]=v;
 v=_mm256_xor_si256(v,r1);t[29]=v;
 v=_mm256_xor_si256(v,r0);t[28]=v;
 v=_mm256_xor_si256(v,r3);t[20]=v;
 v=_mm256_xor_si256(v,r0);t[21]=v;
 v=_mm256_xor_si256(v,r1);t[23]=v;
 v=_mm256_xor_si256(v,r0);t[22]=v;
 v=_mm256_xor_si256(v,r2);t[18]=v;
 v=_mm256_xor_si256(v,r0);t[19]=v;
 v=_mm256_xor_si256(v,r1);t[17]=v;
 v=_mm256_xor_si256(v,r0);t[16]=v;
 v=_mm256_xor_si256(v,r5);t[48]=v;
 v=_mm256_xor_si256(v,r0);t[49]=v;
 v=_mm256_xor_si256(v,r1);t[51]=v;
 v=_mm256_xor_si256(v,r0);t[50]=v;
 v=_mm256_xor_si256(v,r2);t[54]=v;
 v=_mm256_xor_si256(v,r0);t[55]=v;
 v=_mm256_xor_si256(v,r1);t[53]=v;
 v=_mm256_xor_si256(v,r0);t[52]=v;
 v=_mm256_xor_si256(v,r3);t[60]=v;
 v=_mm256_xor_si256(v,r0);t[61]=v;
 v=_mm256_xor_si256(v,r1);t[63]=v;
 v=_mm256_xor_si256(v,r0);t[62]=v;
 v=_mm256_xor_si256(v,r2);t[58]=v;
 v=_mm256_xor_si256(v,r0);t[59]=v;
 v=_mm256_xor_si256(v,r1);t[57]=v;
 v=_mm256_xor_si256(v,r0);t[56]=v;
 v=_mm256_xor_si256(v,r4);t[40]=v;
 v=_mm256_xor_si256(v,r0);t[41]=v;
 v=_mm256_xor_si256(v,r1);t[43]=v;
 v=_mm256_xor_si256(v,r0);t[42]=v;
 v=_mm256_xor_si256(v,r2);t[46]=v;
 v=_mm256_xor_si256(v,r0);t[47]=v;
 v=_mm256_xor_si256(v,r1);t[45]=v;
 v=_mm256_xor_si256(v,r0);t[44]=v;
 v=_mm256_xor_si256(v,r3);t[36]=v;
 v=_mm256_xor_si256(v,r0);t[37]=v;
 v=_mm256_xor_si256(v,r1);t[39]=v;
 v=_mm256_xor_si256(v,r0);t[38]=v;
 v=_mm256_xor_si256(v,r2);t[34]=v;
 v=_mm256_xor_si256(v,r0);t[35]=v;
 v=_mm256_xor_si256(v,r1);t[33]=v;
 v=_mm256_xor_si256(v,r0);t[32]=v;
 v=_mm256_xor_si256(v,r6);t[96]=v;
 v=_mm256_xor_si256(v,r0);t[97]=v;
 v=_mm256_xor_si256(v,r1);t[99]=v;
 v=_mm256_xor_si256(v,r0);t[98]=v;
 v=_mm256_xor_si256(v,r2);t[102]=v;
 v=_mm256_xor_si256(v,r0);t[103]=v;
 v=_mm256_xor_si256(v,r1);t[101]=v;
 v=_mm256_xor_si256(v,r0);t[100]=v;
 v=_mm256_xor_si256(v,r3);t[108]=v;
 v=_mm256_xor_si256(v,r0);t[109]=v;
 v=_mm256_xor_si256(v,r1);t[111]=v;
 v=_mm256_xor_si256(v,r0);t[110]=v;
 v=_mm256_xor_si256(v,r2);t[106]=v;
 v=_mm256_xor_si256(v,r0);t[107]=v;
 v=_mm256_xor_si256(v,r1);t[105]=v;
 v=_mm256_xor_si256(v,r0);t[104]=v;
 v=_mm256_xor_si256(v,r4);t[120]=v;
 v=_mm256_xor_si256(v,r0);t[121]=v;
 v=_mm256_xor_si256(v,r1);t[123]=v;
 v=_mm256_xor_si256(v,r0);t[122]=v;
 v=_mm256_xor_si256(v,r2);t[126]=v;
 v=_mm256_xor_si256(v,r0);t[127]=v;
 v=_mm256_xor_si256(v,r1);t[125]=v;
 v=_mm256_xor_si256(v,r0);t[124]=v;
 v=_mm256_xor_si256(v,r3);t[116]=v;
 v=_mm256_xor_si256(v,r0);t[117]=v;
 v=_mm256_xor_si256(v,r1);t[119]=v;
 v=_mm256_xor_si256(v,r0);t[118]=v;
 v=_mm256_xor_si256(v,r2);t[114]=v;
 v=_mm256_xor_si256(v,r0);t[115]=v;
 v=_mm256_xor_si256(v,r1);t[113]=v;
 v=_mm256_xor_si256(v,r0);t[112]=v;
 v=_mm256_xor_si256(v,r5);t[80]=v;
 v=_mm256_xor_si256(v,r0);t[81]=v;
 v=_mm256_xor_si256(v,r1);t[83]=v;
 v=_mm256_xor_si256(v,r0);t[82]=v;
 v=_mm256_xor_si256(v,r2);t[86]=v;
 v=_mm256_xor_si256(v,r0);t[87]=v;
 v=_mm256_xor_si256(v,r1);t[85]=v;
 v=_mm256_xor_si256(v,r0);t[84]=v;
 v=_mm256_xor_si256(v,r3);t[92]=v;
 v=_mm256_xor_si256(v,r0);t[93]=v;
 v=_mm256_xor_si256(v,r1);t[95]=v;
 v=_mm256_xor_si256(v,r0);t[94]=v;
 v=_mm256_xor_si256(v,r2);t[90]=v;
 v=_mm256_xor_si256(v,r0);t[91]=v;
 v=_mm256_xor_si256(v,r1);t[89]=v;
 v=_mm256_xor_si256(v,r0);t[88]=v;
 v=_mm256_xor_si256(v,r4);t[72]=v;
 v=_mm256_xor_si256(v,r0);t[73]=v;
 v=_mm256_xor_si256(v,r1);t[75]=v;
 v=_mm256_xor_si256(v,r0);t[74]=v;
 v=_mm256_xor_si256(v,r2);t[78]=v;
 v=_mm256_xor_si256(v,r0);t[79]=v;
 v=_mm256_xor_si256(v,r1);t[77]=v;
 v=_mm256_xor_si256(v,r0);t[76]=v;
 v=_mm256_xor_si256(v,r3);t[68]=v;
 v=_mm256_xor_si256(v,r0);t[69]=v;
 v=_mm256_xor_si256(v,r1);t[71]=v;
 v=_mm256_xor_si256(v,r0);t[70]=v;
 v=_mm256_xor_si256(v,r2);t[66]=v;
 v=_mm256_xor_si256(v,r0);t[67]=v;
 v=_mm256_xor_si256(v,r1);t[65]=v;
 v=_mm256_xor_si256(v,r0);t[64]=v;
 v=_mm256_xor_si256(v,r7);t[192]=v;
 v=_mm256_xor_si256(v,r0);t[193]=v;
 v=_mm256_xor_si256(v,r1);t[195]=v;
 v=_mm256_xor_si256(v,r0);t[194]=v;
 v=_mm256_xor_si256(v,r2);t[198]=v;
 v=_mm256_xor_si256(v,r0);t[199]=v;
 v=_mm256_xor_si256(v,r1);t[197]=v;
 v=_mm256_xor_si256(v,r0);t[196]=v;
 v=_mm256_xor_si256(v,r3);t[204]=v;
 v=_mm256_xor_si256(v,r0);t[205]=v;
 v=_mm256_xor_si256(v,r1);t[207]=v;
 v=_mm256_xor_si256(v,r0);t[206]=v;
 v=_mm256_xor_si256(v,r2);t[202]=v;
 v=_mm256_xor_si256(v,r0);t[203]=v;
 v=_mm256_xor_si256(v,r1);t[201]=v;
 v=_mm256_xor_si256(v,r0);t[200]=v;
 v=_mm256_xor_si256(v,r4);t[216]=v;
 v=_mm256_xor_si256(v,r0);t[217]=v;
 v=_mm256_xor_si256(v,r1);t[219]=v;
 v=_mm256_xor_si256(v,r0);t[218]=v;
 v=_mm256_xor_si256(v,r2);t[222]=v;
 v=_mm256_xor_si256(v,r0);t[223]=v;
 v=_mm256_xor_si256(v,r1);t[221]=v;
 v=_mm256_xor_si256(v,r0);t[220]=v;
 v=_mm256_xor_si256(v,r3);t[212]=v;
 v=_mm256_xor_si256(v,r0);t[213]=v;
 v=_mm256_xor_si256(v,r1);t[215]=v;
 v=_mm256_xor_si256(v,r0);t[214]=v;
 v=_mm256_xor_si256(v,r2);t[210]=v;
 v=_mm256_xor_si256(v,r0);t[211]=v;
 v=_mm256_xor_si256(v,r1);t[209]=v;
 v=_mm256_xor_si256(v,r0);t[208]=v;
 v=_mm256_xor_si256(v,r5);t[240]=v;
 v=_mm256_xor_si256(v,r0);t[241]=v;
 v=_mm256_xor_si256(v,r1);t[243]=v;
 v=_mm256_xor_si256(v,r0);t[242]=v;
 v=_mm256_xor_si256(v,r2);t[246]=v;
 v=_mm256_xor_si256(v,r0);t[247]=v;
 v=_mm256_xor_si256(v,r1);t[245]=v;
 v=_mm256_xor_si256(v,r0);t[244]=v;
 v=_mm256_xor_si256(v,r3);t[252]=v;
 v=_mm256_xor_si256(v,r0);t[253]=v;
 v=_mm256_xor_si256(v,r1);t[255]=v;
 v=_mm256_xor_si256(v,r0);t[254]=v;
 v=_mm256_xor_si256(v,r2);t[250]=v;
 v=_mm256_xor_si256(v,r0);t[251]=v;
 v=_mm256_xor_si256(v,r1);t[249]=v;
 v=_mm256_xor_si256(v,r0);t[248]=v;
 v=_mm256_xor_si256(v,r4);t[232]=v;
 v=_mm256_xor_si256(v,r0);t[233]=v;
 v=_mm256_xor_si256(v,r1);t[235]=v;
 v=_mm256_xor_si256(v,r0);t[234]=v;
 v=_mm256_xor_si256(v,r2);t[238]=v;
 v=_mm256_xor_si256(v,r0);t[239]=v;
 v=_mm256_xor_si256(v,r1);t[237]=v;
 v=_mm256_xor_si256(v,r0);t[236]=v;
 v=_mm256_xor_si256(v,r3);t[228]=v;
 v=_mm256_xor_si256(v,r0);t[229]=v;
 v=_mm256_xor_si256(v,r1);t[231]=v;
 v=_mm256_xor_si256(v,r0);t[230]=v;
 v=_mm256_xor_si256(v,r2);t[226]=v;
 v=_mm256_xor_si256(v,r0);t[227]=v;
 v=_mm256_xor_si256(v,r1);t[225]=v;
 v=_mm256_xor_si256(v,r0);t[224]=v;
 v=_mm256_xor_si256(v,r6);t[160]=v;
 v=_mm256_xor_si256(v,r0);t[161]=v;
 v=_mm256_xor_si256(v,r1);t[163]=v;
 v=_mm256_xor_si256(v,r0);t[162]=v;
 v=_mm256_xor_si256(v,r2);t[166]=v;
 v=_mm256_xor_si256(v,r0);t[167]=v;
 v=_mm256_xor_si256(v,r1);t[165]=v;
 v=_mm256_xor_si256(v,r0);t[164]=v;
 v=_mm256_xor_si256(v,r3);t[172]=v;
 v=_mm256_xor_si256(v,r0);t[173]=v;
 v=_mm256_xor_si256(v,r1);t[175]=v;
 v=_mm256_xor_si256(v,r0);t[174]=v;
 v=_mm256_xor_si256(v,r2);t[170]=v;
 v=_mm256_xor_si256(v,r0);t[171]=v;
 v=_mm256_xor_si256(v,r1);t[169]=v;
 v=_mm256_xor_si256(v,r0);t[168]=v;
 v=_mm256_xor_si256(v,r4);t[184]=v;
 v=_mm256_xor_si256(v,r0);t[185]=v;
 v=_mm256_xor_si256(v,r1);t[187]=v;
 v=_mm256_xor_si256(v,r0);t[186]=v;
 v=_mm256_xor_si256(v,r2);t[190]=v;
 v=_mm256_xor_si256(v,r0);t[191]=v;
 v=_mm256_xor_si256(v,r1);t[189]=v;
 v=_mm256_xor_si256(v,r0);t[188]=v;
 v=_mm256_xor_si256(v,r3);t[180]=v;
 v=_mm256_xor_si256(v,r0);t[181]=v;
 v=_mm256_xor_si256(v,r1);t[183]=v;
 v=_mm256_xor_si256(v,r0);t[182]=v;
 v=_mm256_xor_si256(v,r2);t[178]=v;
 v=_mm256_xor_si256(v,r0);t[179]=v;
 v=_mm256_xor_si256(v,r1);t[177]=v;
 v=_mm256_xor_si256(v,r0);t[176]=v;
 v=_mm256_xor_si256(v,r5);t[144]=v;
 v=_mm256_xor_si256(v,r0);t[145]=v;
 v=_mm256_xor_si256(v,r1);t[147]=v;
 v=_mm256_xor_si256(v,r0);t[146]=v;
 v=_mm256_xor_si256(v,r2);t[150]=v;
 v=_mm256_xor_si256(v,r0);t[151]=v;
 v=_mm256_xor_si256(v,r1);t[149]=v;
 v=_mm256_xor_si256(v,r0);t[148]=v;
 v=_mm256_xor_si256(v,r3);t[156]=v;
 v=_mm256_xor_si256(v,r0);t[157]=v;
 v=_mm256_xor_si256(v,r1);t[159]=v;
 v=_mm256_xor_si256(v,r0);t[158]=v;
 v=_mm256_xor_si256(v,r2);t[154]=v;
 v=_mm256_xor_si256(v,r0);t[155]=v;
 v=_mm256_xor_si256(v,r1);t[153]=v;
 v=_mm256_xor_si256(v,r0);t[152]=v;
 v=_mm256_xor_si256(v,r4);t[136]=v;
 v=_mm256_xor_si256(v,r0);t[137]=v;
 v=_mm256_xor_si256(v,r1);t[139]=v;
 v=_mm256_xor_si256(v,r0);t[138]=v;
 v=_mm256_xor_si256(v,r2);t[142]=v;
 v=_mm256_xor_si256(v,r0);t[143]=v;
 v=_mm256_xor_si256(v,r1);t[141]=v;
 v=_mm256_xor_si256(v,r0);t[140]=v;
 v=_mm256_xor_si256(v,r3);t[132]=v;
 v=_mm256_xor_si256(v,r0);t[133]=v;
 v=_mm256_xor_si256(v,r1);t[135]=v;
 v=_mm256_xor_si256(v,r0);t[134]=v;
 v=_mm256_xor_si256(v,r2);t[130]=v;
 v=_mm256_xor_si256(v,r0);t[131]=v;
 v=_mm256_xor_si256(v,r1);t[129]=v;
 v=_mm256_xor_si256(v,r0);t[128]=v;
}

static void fallback_product(int n,const bool*A,const bool*B,bool*C){
 if(n<=0)return;
 if(n%256){for(int i=0;i<n;i++)for(int j=0;j<n;j++){bool s=0;for(int k=0;k<n;k++)s^=A[(size_t)i*n+k]&B[(size_t)k*n+j];C[(size_t)i*n+j]=s;}return;}
 const int w=n/64,groups=n/8;
 // ARM IXS: the parent stores the per-group index rows at stride n, i.e. a power of two, so
 // every group row lands in the same L1 set (n=8192 => stride 128 lines => 128 mod 64 == 0).
 // Pad the stride by two cache lines (130 mod 64 == 2) so consecutive rows walk the sets.
 const size_t IXS=(size_t)n+128;
 uint64_t*pb;unsigned char*ix;__m256i*pc;
 if(posix_memalign((void**)&pb,4096,(size_t)n*w*8)||posix_memalign((void**)&ix,4096,IXS*(size_t)groups+4096)||posix_memalign((void**)&pc,4096,(size_t)n*32))abort();
 for(int k=0;k<n;k++)for(int j=0;j<n;j+=64){
  __m256i x=_mm256_load_si256((const __m256i*)(B+(size_t)k*n+j));
  __m256i y=_mm256_load_si256((const __m256i*)(B+(size_t)k*n+j+32));
  uint32_t lo=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(x,7));
  uint32_t hi=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(y,7));
  pb[(size_t)k*w+j/64]=(uint64_t)lo|((uint64_t)hi<<32);
 }
 // ARM IXT: build the index table tile by tile.  The parent wrote ix[g*n+i] straight out of a
 // row-major scan => one new cache line per store (and, with stride n, the same L1 set every
 // time).  Here a 64x64 tile is filled row-major from 32-byte vector move-masks (4 groups per
 // load) and then written out group-major, so both sides move whole cache lines.
 { alignas(64) static unsigned char tile[64*64];   /* tile[c*64+r] 列主序 */
   for(int ib=0;ib<n;ib+=64){
     for(int gb=0;gb<groups;gb+=64){
       int ge=gb+64; if(ge>groups) ge=groups;
       for(int i=ib;i<ib+64;i++){
         const bool* Ai=A+(size_t)i*n; unsigned char* Ti=tile+(i-ib);
         for(int gg=gb;gg<ge;gg+=4){
           __m256i v=_mm256_loadu_si256((const __m256i*)(Ai+(size_t)gg*8));
           uint32_t m=(uint32_t)_mm256_movemask_epi8(_mm256_slli_epi16(v,7));
           Ti[(gg-gb+0)*64]=(unsigned char)m;
           Ti[(gg-gb+1)*64]=(unsigned char)(m>>8);
           Ti[(gg-gb+2)*64]=(unsigned char)(m>>16);
           Ti[(gg-gb+3)*64]=(unsigned char)(m>>24);
         }
       }
       for(int gg=gb;gg<ge;gg++) memcpy(ix+(size_t)gg*IXS+ib, tile+(gg-gb)*64, 64);
     }
   }
 }
 alignas(64) __m256i tab[4][256];
 const __m256i shuf=_mm256_setr_epi8(0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3);
 const __m256i bits=_mm256_set1_epi64x(0x8040201008040201ULL),one=_mm256_set1_epi8(1),zero=_mm256_setzero_si256();
 for(int jb=0;jb<w;jb+=4){
  for(int g=0;g<groups;g+=4){
   const int cnt=groups-g<4?groups-g:4;
   for(int z=0;z<cnt;z++){
    __m256i*Tz=tab[z];
    build_gray8(pb+(size_t)((g+z)*8)*w+jb,w,Tz);
   }
   const unsigned char*x=ix+(size_t)g*IXS;
   if(g==0){
   if(cnt==4)apply_m<4,true>( n,pc,tab,x,IXS);else if(cnt==3)apply_m<3,true>(n,pc,tab,x,IXS);
   else if(cnt==2)apply_m<2,true>(n,pc,tab,x,IXS);
   else if(cnt==1)apply_m<1,true>(n,pc,tab,x,IXS);
   } else {
   if(cnt==4)apply_m<4,false>( n,pc,tab,x,IXS);else if(cnt==3)apply_m<3,false>(n,pc,tab,x,IXS);
   else if(cnt==2)apply_m<2,false>(n,pc,tab,x,IXS);
   else if(cnt==1)apply_m<1,false>(n,pc,tab,x,IXS);
   }
  }
  for(int i=0;i<n;i++){
   const uint32_t*row=(const uint32_t*)(pc+i);
   for(int t=0;t<8;t++){
    __m256i x=_mm256_shuffle_epi8(_mm256_set1_epi32(row[t]),shuf);
    x=_mm256_min_epu8(_mm256_and_si256(x,bits),one);
    _mm256_stream_si256((__m256i*)(C+(size_t)i*n+jb*64+t*32),x);
   }
  }
 }
 _mm_sfence();
 free(pc);free(ix);free(pb);
}

namespace bit_winograd {
alignas(64) static __m256i tables[4][256];
alignas(64) static __m256i accum[4096];
alignas(64) static unsigned char indices[4096*512+512*128];
alignas(64) static unsigned char tile[64*64];
static __attribute__((noinline)) void leaf(int n,const uint64_t*a,int sa,const uint64_t*b,int sb,uint64_t*c,int sc){
 const int groups=n/8,words=n/64;const size_t is=(size_t)n+128;
 for(int ib=0;ib<n;ib+=64)for(int gb=0;gb<groups;gb+=64){
  const int ge=gb+64<groups?gb+64:groups;
  for(int i=ib;i<ib+64;i++){
   const unsigned char*ai=(const unsigned char*)(a+(size_t)i*sa);
   unsigned char*ti=tile+i-ib;
   for(int gg=gb;gg<ge;gg+=4){
    uint32_t v;memcpy(&v,ai+gg,4);
    ti[(gg-gb)*64]=v;ti[(gg-gb+1)*64]=v>>8;
    ti[(gg-gb+2)*64]=v>>16;ti[(gg-gb+3)*64]=v>>24;
   }
  }
  for(int gg=gb;gg<ge;gg++)memcpy(indices+(size_t)gg*is+ib,tile+(gg-gb)*64,64);
 }
 for(int jb=0;jb<words;jb+=4){
  for(int g=0;g<groups;g+=4){
   int cnt=groups-g<4?groups-g:4;
   for(int z=0;z<cnt;z++)build_gray8(b+(size_t)(g+z)*8*sb+jb,sb,tables[z]);
   const unsigned char*x=indices+(size_t)g*is;
   if(g==0){
    if(cnt==4)apply_m<4,true>( n,accum,tables,x,is);else if(cnt==3)apply_m<3,true>(n,accum,tables,x,is);
    else if(cnt==2)apply_m<2,true>(n,accum,tables,x,is);
    else apply_m<1,true>(n,accum,tables,x,is);
   }else{
    if(cnt==4)apply_m<4,false>( n,accum,tables,x,is);else if(cnt==3)apply_m<3,false>(n,accum,tables,x,is);
    else if(cnt==2)apply_m<2,false>(n,accum,tables,x,is);
    else apply_m<1,false>(n,accum,tables,x,is);
   }
  }
  for(int i=0;i<n;i++)_mm256_store_si256((__m256i*)(c+(size_t)i*sc+jb),accum[i]);
 }
}
static void form(int n,uint64_t*d,const uint64_t*a,int sa,const uint64_t*b,int sb){
 const int w=n/64;
 for(int i=0;i<n;i++)for(int j=0;j<w;j+=4)
  _mm256_store_si256((__m256i*)(d+(size_t)i*w+j),_mm256_xor_si256(_mm256_load_si256((const __m256i*)(a+(size_t)i*sa+j)),_mm256_load_si256((const __m256i*)(b+(size_t)i*sb+j))));
}
static void reconstruct(int h,uint64_t*c,int sc,const uint64_t*m5,const uint64_t*m6,const uint64_t*m7){
 const int w=h/64;
 for(int i=0;i<h;i++)for(int j=0;j<w;j+=4){
  uint64_t*c1=c+(size_t)i*sc+j,*c2=c1+(size_t)h*sc,*c3=c1+w,*c4=c2+w;
  size_t z=(size_t)i*w+j;
  __m256i p1=_mm256_load_si256((const __m256i*)c1),p2=_mm256_load_si256((const __m256i*)c2),p3=_mm256_load_si256((const __m256i*)c3),p4=_mm256_load_si256((const __m256i*)c4);
  __m256i p5=_mm256_load_si256((const __m256i*)(m5+z)),p6=_mm256_load_si256((const __m256i*)(m6+z)),p7=_mm256_load_si256((const __m256i*)(m7+z));
  __m256i u2=_mm256_xor_si256(p1,p6),u3=_mm256_xor_si256(u2,p7);
  _mm256_store_si256((__m256i*)c1,_mm256_xor_si256(p1,p2));
  _mm256_store_si256((__m256i*)c3,_mm256_xor_si256(_mm256_xor_si256(u2,p5),p3));
  _mm256_store_si256((__m256i*)c2,_mm256_xor_si256(u3,p4));
  _mm256_store_si256((__m256i*)c4,_mm256_xor_si256(u3,p5));
 }
}
static void recurse(int n,const uint64_t*a,int sa,const uint64_t*b,int sb,uint64_t*c,int sc,uint64_t*work){
 if(n<=4096){leaf(n,a,sa,b,sb,c,sc);return;}
 int h=n/2,w=h/64;size_t count=(size_t)h*w;
 uint64_t*s=work,*t=s+count,*m5=t+count,*m6=m5+count,*m7=m6+count,*next=m7+count;
 const uint64_t*a11=a,*a12=a+w,*a21=a+(size_t)h*sa,*a22=a21+w;
 const uint64_t*b11=b,*b12=b+w,*b21=b+(size_t)h*sb,*b22=b21+w;
 uint64_t*c1=c,*c2=c+(size_t)h*sc,*c3=c+w,*c4=c2+w;
 recurse(h,a12,sa,b21,sb,c2,sc,next);
 recurse(h,a11,sa,b11,sb,c1,sc,next);
 form(h,t,b12,sb,b11,sb);form(h,s,a21,sa,a22,sa);
 recurse(h,s,w,t,w,m5,w,next);
 form(h,s,s,w,a11,sa);form(h,t,b22,sb,t,w);
 recurse(h,s,w,t,w,m6,w,next);
 form(h,s,a12,sa,s,w);form(h,t,t,w,b21,sb);
 recurse(h,a22,sa,t,w,c4,sc,next);
 recurse(h,s,w,b22,sb,c3,sc,next);
 form(h,t,b22,sb,b12,sb);form(h,s,a11,sa,a21,sa);
 recurse(h,s,w,t,w,m7,w,next);
 reconstruct(h,c,sc,m5,m6,m7);
}
static void pack(int n,const bool*a,uint64_t*d){
 for(size_t z=0,sz=(size_t)n*n;z<sz;z+=64){
  __m256i x=_mm256_load_si256((const __m256i*)(a+z));
  __m256i y=_mm256_load_si256((const __m256i*)(a+z+32));
  uint32_t lo=_mm256_movemask_epi8(_mm256_slli_epi16(x,7)),hi=_mm256_movemask_epi8(_mm256_slli_epi16(y,7));
  d[z/64]=(uint64_t)lo|((uint64_t)hi<<32);
 }
}
}
void matrix_multiply(int n,const bool*A,const bool*B,bool*C){
 if(n<4096 || (n&(n-1))){fallback_product(n,A,B,C);return;}
 size_t count=(size_t)n*n/64;uint64_t*mem;
 if(posix_memalign((void**)&mem,4096,count*5*sizeof(uint64_t)))abort();
 uint64_t*a=mem,*b=a+count,*c=b+count,*work=c+count;
 bit_winograd::pack(n,A,a);bit_winograd::pack(n,B,b);
 bit_winograd::recurse(n,a,n/64,b,n/64,c,n/64,work);
 const __m256i shuf=_mm256_setr_epi8(0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3);
 const __m256i bits=_mm256_set1_epi64x(0x8040201008040201ULL),one=_mm256_set1_epi8(1);
 for(size_t z=0;z<count*2;z++){
  __m256i v=_mm256_shuffle_epi8(_mm256_set1_epi32(((const uint32_t*)c)[z]),shuf);
  v=_mm256_min_epu8(_mm256_and_si256(v,bits),one);
  _mm256_stream_si256((__m256i*)(C+z*32),v);
 }
 _mm_sfence();free(mem);
}

评测结果