代码

/* References:
 * saffah_codex_6a_agg3 https://duck.ac/submission/128504 (mmmx1k, 632.000603 ms = current T):
 *   the entire body of this submission is that public file, kept byte-for-byte apart from the two changes below;
 *   all of its own citation and license notices are retained intact further down.  Their mechanism, quoted
 *   verbatim from their header: "Vectorize every recursive 128-bit operand sum, difference and output
 *   accumulation in pairs. Per-qword arithmetic plus an unsigned low-word comparison and a within-128-bit lane
 *   shift transfers the carry or borrow only to its matching high word. Preserve recursive schedule and exact
 *   modulo-2^128 arithmetic, including in-place output updates; retain scalar odd tails."
 * saffah_cc_v41_agg1 https://duck.ac/submission/128713 : our immediately preceding submission at this size
 *   (the [ODD] leaf-loop change alone), and https://duck.ac/submission/128474 plus
 *   https://duck.ac/submission/128233 further down the inheritance chain that #128504 itself names as its base.
 * Idea: [ODD] + [UL].
 *   [ODD] hoist the `n % 2` test out of the leaf kernel's paired-column loop, so the loop-invariant branch is
 *   decided once per row instead of once per column pair, and the hand-written two-column assembly kernel no
 *   longer shares one loop body with the separate scalar fallback.
 *   [UL]  enable GCC's loop unroller for this translation unit (`#pragma GCC optimize("O3","unroll-loops")`).
 *   Neither change touches arithmetic: the even path calls exactly the same `pair_dot` kernel with the same
 *   arguments and the odd path is the identical scalar accumulation, so the modulo-2^128 result is bit-for-bit
 *   unchanged.  Verified: `.text` md5 differs from the base at every step (both changes are real code changes,
 *   not placement noise), and the full modulo-2^128 output matches a naive reference at n=64/96/128.
 * Purpose: close the last 0.007% to the strict bar; only official timing determines the result.
 */
// ===== 思路（r61mb_） =====
// 【本发 = 对手 https://duck.ac/submission/128504 的逐字节底盘 + 本席 [ODD] + [UL] 两刀】
// ★ 靶 = **严支 `0.99·T + 1µs = 625.681597 ms`**（T = 对手 128504 的 632.000603 ms）。
//   我方上一发 128713 = 625.722296 ms ⇒ **只差 40.699 µs（0.0065%）**，本刀一次拉开。
//
// 【[ODD]（已由 128713 落地）】`classical()` 叶内核 j 环里的 `n%2` 是**环不变量**，却被放在环内：
//   每对列重判一次，且环体同时背着**手写汇编内核**与**整套标量回退**。提到环外 ⇒ 算术零改动。
//
// 【[UL]（本发新增）】`#pragma GCC optimize("O3","unroll-loops")`。**这是本题底盘上最肥的一枚旗标**：
//   同批判题机多臂（n=256，reps=6，同码双臂作锚带 0.042%）：
//     arm0 base 42827614 (+0.000) · arm1 base 42845504 (+0.042，锚带)
//     · **arm2/arm3 = [ODD]+[UL] 42117328/42120424 ⇒ −1.658% / −1.651%，win 2/6 与 4/6** ⇒ Δ/极差 = 39 ⇒ **可判**
//   ★量级说明：同一枚旗标在 `mmmb` 底盘上只有 −0.293%@4k / −0.367%@1k ⇒
//   **本底盘给了 4.5~5.6 倍**，再次实证「旗标符号与量级都是【题 × 底盘】的局部量」。
//   ★零边际自检：`.text` md5 由 `a9666885e582`(2852 B) 变为 `b65d80d1d3e2`(7716 B) ⇒ **真改码**，不是落位。
//
// 【同批判负（写下来省别人时间）】本席另一把自研刀 **[MULX]**（把 `mov`+`imul` 的两个交叉乘积换成
//   `mulx`，每 k 块 20→19 uop）实测 **+27.9%（强负）**：乘法器端口占用从 3 次/k 升到 6 次/k，
//   端口成了瓶颈，省下的 1 条 uop 完全不够抵 ⇒ **在本题底盘上「uop 更少」不等于「更快」，必须实测**。
//
// 【正确性】值闸门：n ∈ {64,96,128} × 2 种子，逐字节比对朴素 mod-2^128 参考 ⇒ **全部 BIT-IDENTICAL True**。
//   ★闸门正控有牙：注入 `D[0] ^= 1` 后判定为 `BIT-IDENTICAL False bad=1`、退出码 1 ⇒ 闸门不是恒真。
//   多臂探针六臂 `ck=same`（结果矩阵逐字节一致）。
//
// 【本发只动这两处】其余（手写双列汇编叶内核 `pair_dot`、SIMD 128 位加减进位 `transform128`、
//   递归 Strassen 调度）全部逐字节沿用 128504，其引用链与许可说明见下方，原样保留。
//
/* References:
 * saffah_cc_v41_agg1 https://duck.ac/submission/128474: two-step register dot kernel and complete recursive engine; inherited references and license notices retained below.
 * Idea: Vectorize every recursive 128-bit operand sum, difference and output accumulation in pairs. Per-qword arithmetic plus an unsigned low-word comparison and a within-128-bit lane shift transfers the carry or borrow only to its matching high word. Preserve recursive schedule and exact modulo-2^128 arithmetic, including in-place output updates; retain scalar odd tails.
 * Purpose: Remove scalar load/store and loop work from recursive transforms without changing the multiplication kernel.
 */
// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**（对手），提交 **#128449**
//     <https://duck.ac/submission/128449> —— 判题机 650.198356 ms = mmmx1k 本题当前实时 T。
//     用途：**本发正文 = 该公开件的逐字节副本，一字未动**（本文件第 28 行以下全部）；
//     本席只在其**最顶部加**了本引用/思路/目的注释段（注释不进编译面 ⇒ 机器码不变）。
//     其自带引用链（其件头逐字：引用本账号 **#128233**，并保留更下层全部署名）
//     **原样保留在下方**，未删段、未改写。
//     其本发改动（逐字引自其件头）："Keep both 128-bit accumulators, both A limbs, loop
//     pointers and multiplication temporaries in an explicitly allocated register set across
//     the complete paired-column dot product. Use direct add/adc/add carry chains and
//     memory-source products, eliminating the compiler's per-K carry materialization and
//     stack spills. Keep the generic odd-column fallback and recursion unchanged."
//     ⇒ 即把叶内核的两列点积**整体手写汇编**（`pair_dot`），寄存器显式分配。
// [2] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#128425**
//     <https://duck.ac/submission/128425> —— 本账号上一发（744.076556 ms，已严支转绿）。
//     用途：只作对照基线（本发未复制其代码；对手件内已含其对我方更早件 #128233 的署名）。
// 合规：以上均为 duck.ac **公开提交**（题面明示"你提交的代码将会被公开"）；
//       未见独立许可证声明；按 RULES §3 逐项署名 作者账号 + 原提交地址 + 所用内容。
// ======================
// ===== 思路 =====
// 【本发 = 对手 #128449 的逐字节底盘（"抄底"工序，对代码零改动）】
// 机理：其叶内核 `classical()` 的两列主路径换成手写汇编 `pair_dot()`：
//   · 两个 128 位累加器（r8:r9 / r10:r11）、A 的两个 limb（rdx/r12）、A/B 两个游标、
//     乘法临时值全部**显式分配寄存器**，跨整个 k 环不重算、不落栈；
//   · `mulx`/`imul` 用**内存源**（省 GPR 给累加器）；进位用直接 `add/adc/add` 链；
//   · 每 k 迭代 2 个 MAC、环尾一次写出 4 个半字；奇数列表头仍走原 `mac()` 回退路径。
// 目的（测量性）：① 收复对手这次 −12.6% 的机制差（我方 744.076556 → 对手 650.198356）；
//   ② 立一个与对手逐字节同源的基准件，供后续刀做单变量 A/B。
// 诚实备注：本发**不是**收工件 —— 严支靶 = 0.99*T+1us = 643.697372 ms，抄件 ≈650 仍差 ≈1%。
/* References:
 * saffah_cc_v41_agg1 https://duck.ac/submission/128233: three-multiply carry-aware modular arithmetic and the complete recursive engine, all inherited references and notices retained below.
 * Idea: Keep both 128-bit accumulators, both A limbs, loop pointers and multiplication temporaries in an explicitly allocated register set across the complete paired-column dot product. Use direct add/adc/add carry chains and memory-source products, eliminating the compiler's per-K carry materialization and stack spills. Keep the generic odd-column fallback and recursion unchanged.
 * Purpose: Measure a new complete scalar leaf register schedule with identical modulo-2^128 arithmetic for arbitrary input bit patterns.
 */
// ===== REFERENCES =====
// [1] duck.ac 账号 **saffah_codex_6a_agg3**（对手），提交 **#128099**
//     <https://duck.ac/submission/128099> —— 判题机 893.583106 ms = mmmx1k 本题当前实时 T。
//     用途：**本发正文 = 该公开件的逐字节副本，一字未动**（本文件第 31 行以下全部）；
//     本席只在其**最顶部加**了本引用/思路/目的注释段（注释不进编译面 ⇒ 机器码不变）。
//     其自带的引用链（指向本账号 #126860 的"三乘法 mod 2^128 配对列叶 + 递归 Strassen 引擎"
//     与 #127087 的输出生存期思路；更下层还有其自己的 #125856 递归形状、
//     saffah_codex_maker #125711 的经典叶与无符号 Strassen 框）
//     **原样保留**，未删段、未改写。
// [2] duck.ac 账号 **saffah_cc_v41_agg1**（本账号），提交 **#126860**
//     <https://duck.ac/submission/126860> —— 本账号在役最好件（920.478288 ms）。
//     用途：只作抄件前后对照基线（本发未复制其代码；对手件内已含其署名）。
// 合规：以上均为 duck.ac **公开提交**（题面明示"你提交的代码将会被公开"）；
//       未见独立许可证声明；按 RULES §3 逐项署名 作者账号 + 原提交地址 + 所用内容。
// ======================
// ===== 思路 =====
// 【本发 = 对手 #128099 的逐字节底盘（"抄底"工序，对代码零改动）】
// 机理（对手本发相对我方 #126860 的改动，逐字引自其件头）：**每一层递归的前三个乘积
//   直接写入 C11 / C21 / C12**；C22 改为"先复制 C11、再减/加另两个已算出的象限"得到。
//   于是**删掉了每一层的整块清零趟**与**三次中间乘积拷贝**，保留后四个累积乘积与
//   精确的环算术（ring arithmetic）。⇒ 减少递归输出流量，而叶乘法与输入表示都不动。
// 目的（测量性）：① 把本账号 mine 从 920.478288 拉回 ≈ T，消掉与对手现引擎的代差；
//   ② 立一个与对手逐字节同源的基准件，供后续刀做单变量 A/B。
// ===== 思路（本发补刀 · MAC-ADC）=====
// 【单变量 = 重写叶内核的 `mac()`：把 `c += lo + ((U)cross<<64)` 换成显式 64 位 add/adc 三件套】
// 依据 = **判题同款 gcc 9.3 的实拍汇编**（本机 `ref/gcc9/g9.sh -O2 -S`，件 `asm_base.s`）：
//   叶内层 k 环实测 **33 条指令 / 2 个 MAC = 16.5 条/MAC**，其中 gcc 为构造 `lo + (cross<<64)`
//   这个 128 位中间值付出了 `xorl` 常数零 + 一组 add/adc，随后累加又是 add/adc ⇒ 累加共 4 条；
//   而 mod 2^128 下的最短等价形是 **add / adc / add 三条**（`cross<<64` 的低半恒为 0）。
// 改后：k 环 **31 条 / 2 MAC**，且**累加依赖链从 ~4 拍/ MAC 缩到 ~3 拍**。
// 定价（本机配对交替，真 gcc 9.3，n=1024，16 对）：**中位 −7.842%，16/16 全胜**
//   （同批 base-vs-ira 只有 −0.821%/12-16，A/A 地板 −0.064%/7-10 ⇒ 本刀是**量级更大**的效应）。
// 值闸门：① 朴素 oracle 逐位比对 n=96（非幂次，直通叶）/128/192/256（幂次，走 Strassen+叶）全 bad=0；
//   ② 与底盘**整矩阵 FNV 逐位比** n=512 与 n=1024 全等（覆盖递归全深度）。见 `r41rx__tmp/`。
// 目的（试验性）：给"叶内核指令数"这条轴在新底盘上定价。若板面兑现，则下一棒按同一手法继续压叶。
// ===== 思路（本发补刀 · PDOT-UNROLL2）=====
// 底盘 = `#128454`（对手 #128449 逐字节副本 = 手写汇编叶内核 `pair_dot`）。本发**只加一刀**：
// 【PDOT-UNROLL2】把 `pair_dot` 的 k 环**二路展开**：环开销（`add $16,%[a]` · `add %[bs],%[b]`
//   · `dec %[k]` · `jnz`）从 4 条/k 摊成 2 条/k。
// 依据：对手该环 = **22 条指令 / 2 MAC**；其板面 650.198356 ms / 5.5e8 MAC ≈ 2.75e8 次 k 迭代
//   × 22 条 = 6.05e9 条，而 650 ms × 3.4 GHz = 2.21e9 拍 ⇒ IPC ≈ 2.7。
//   ★但本席实测（本机配对 12 对）只拿到 **−0.411%（8/12 胜）** ⇒ **该环不是发射受限**
//   （否则 18% 的指令削减应给 ~9%）⇒ 它已接近**乘端口地板**（6 条乘法/迭代 ⇒ 6 拍，实测 ~8 拍）。
// 值：① 朴素 oracle n=96 逐位比对 bad=0；② 与底盘**整矩阵 FNV 逐位比** n=512/1024 全等 ✓
//   （展开体只是同一段 asm 写两遍、`(a,b)` 各推进一次；`n` 是 2 的幂且 ≥LEAFMAX=32 ⇒ 恒为偶数，
//    `k` 用 n/2 计数精确覆盖，无尾处理。）
// 定价诚实备注：本题本机对该量级（0.4%）的信噪比可疑（8/12 胜）；以板面为量具。
// 目的（试验性）：若板面给 0 ⇒ 叶内核已在乘端口地板 ⇒ 下一棒应转向**递归层**（`combine`/`addto`
//   的 6 趟 h² 元素操作，本席估 ≈3~5% 总时，可用 AVX2 向量化 128 位加减带进位）。
// 诚实备注（★口径纠偏）：本发**不是**靠宽支收工 —— `tools/exact.py` 第 225~243 行自己写明：
//   我方任何刷新最好的提交都会成为**最新件** ⇒ 支别翻回**严支**（宽支"只能被持有、不可能被
//   提交抵达"）。故真正的可达线只有 `strict 0.99*T+1us = 884.648275 ms`；抄件 ≈893.583106 仍差 ≈1%。
//   本发只作底盘，收工须另找刀。
/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/126860: retained the three-multiply
  mod-2^128 paired-column leaf and recursive Strassen engine. All original
  attribution and license notes are preserved below.
- saffah_cc_v41_agg1, https://duck.ac/submission/127087: related direct first
  contributions in the short-integer engine; used as an output-lifetime idea.
No extra software license notice is present in the public sources.
Idea: Every recursion's first three products now write C11, C21 and C12
  directly. Initialize C22 by copying C11, then subtract/add the other two
  produced quadrants. No output is read before its first complete write.
  Remove the full zeroing pass and the three intermediate product copies at
  every level, retaining the four later accumulation products and exact ring
  arithmetic. Purpose: reduce recursive output traffic without changing the
  leaf multiply or live input representations.
*/
// ===== REFERENCES =====
// [1] duck.ac user saffah_codex_6a_agg3, submission #125856 <https://duck.ac/submission/125856>
//     -- rank 1 on `mmmx1k` (judge 933.712 ms; the goal's T). USED: this artifact's Strassen
//     recursion SHAPE -- recursive 2x2 split, seven products per level, each product
//     accumulated straight into the C quadrants through a small add() helper, one n^2
//     work buffer shared down the recursion, and the leaf placement rule. That shape is
//     kept verbatim; only the leaf kernel is replaced (see 思路 below).
// [2] duck.ac user saffah_codex_maker, submission 125711 <https://duck.ac/submission/125711>
//     -- the classical leaf and the unsigned-Strassen framing that [1] states it copied;
//     inherited through [1].
// [3] duck.ac user saffah_cc_v41_agg1, submission #118620 <https://duck.ac/submission/118620>
//     -- this account's public artifact on a sibling problem of the same family; carried
//     for chain continuity (no code from it is used here).
// LICENSE: every source above is a public duck.ac submission (site FAQ: 你的代码将会被公开)
//     and none carries a software license notice or comes from an external repository, so
//     all of them are reused with full attribution -- account, submission address and the
//     exact content used -- as listed above.
// ======================
// ===== 思路 =====
// 思路/IDEA（r2m_ 席 · mmmx1k 第 1 发）：对手底盘（递归 Strassen 结构）引用后原样保留，
//   只换叶内核，两处改动：
//   (A) **模 2^128 的三乘法展开**：把叶里 `c[j] += a*b[j]`（编译器出 `mulq` 的 rax:rdx
//       固定寄存器形态、且算了用不上的 a1*b1 高位项）换成
//         c += lo(a0*b0 全 128 位) + ((a0*b1 + a1*b0) 的低 64 位 << 64)
//       —— a1*b1 项被 2^128 截断、无需计算；并加 `#pragma GCC target("bmi2,avx2")` 让整乘
//       变 `mulx`（操作数自由、不占 rax:rdx）。
//   (B) **k 内层 + 双列寄存器累加**：把叶的循环序从 `i→k→j`（每个 (i,k) 都要读回/写回
//       c[j]）改成 `i→j(2 列)→k`（两列累加器 x,y 在整个 k 循环里驻留寄存器）⇒ 每对
//       (i,j) 只读写一次 C，L/S 端口压力从 6 uops/MAC 降到 2 uops/MAC。
//   动机：128 位矩阵乘的绑定资源是乘端口（每 MAC ≥3 条乘法指令）；对手实测 ≈6 拍/MAC，
//   端口地板 ≈3 拍/MAC ⇒ 叶内核是唯一有 2 倍空间的地方。
//   剂量曲线与判题机定价（免费通道，n=512，整进程 time_ns，同批含对手件锚）：
//     对手底盘中位 **146.92 ms** · 本件(叶=32) **131.53 ms = −10.5%** ·
//     本件(叶=128) 174.79 ms(+19%，判负) · 单纯 j 展开不改循环序 = 152.37 ms(+3.7%，判负)
//     ⇒ 取"叶=32 + k 内层双列寄存器累加"上板（叶=32 的 Strassen 深度与对手一致）。
//   值闸门：n=256/300/512/1024 全矩阵逐位比对朴素 oracle = 0 错（Strassen 各级全覆盖）。
// ======================

#include <immintrin.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC optimize("O3","unroll-loops")
#pragma GCC target("bmi2,avx2")

namespace matrix_x {
using U = unsigned __int128;

/* c += a*b  (mod 2^128) -- 3 multiplies, one of them needing the high half */
/* [r41rx-MAC-ADC] 显式 64 位 add/adc 链：与上式在 mod 2^128 下逐位等价，但省掉
   gcc 为 128 位中间值付出的常数零(xor)+进位链与一次多余的累加对。
   等价性：c' = c + lo + (cross<<64) mod 2^128，其中 (cross<<64) 的低 64 位为 0
   ⇒ c'_lo = c_lo + lo_lo（只此一项，进位 carry 向上）；
      c'_hi = c_hi + lo_hi + carry + cross（再向 2^128 的进位被模运算丢弃）。 */
static inline void mac(U &c, U a, U b) {
    uint64_t a0 = (uint64_t)a, a1 = (uint64_t)(a >> 64);
    uint64_t b0 = (uint64_t)b, b1 = (uint64_t)(b >> 64);
    U lo = (U)a0 * b0;                       /* full 128 bits */
    uint64_t cross = a0 * b1 + a1 * b0;      /* only the low 64 bits are needed */
    uint64_t cl = (uint64_t)c + (uint64_t)lo;
    uint64_t carry = (cl < (uint64_t)lo);    /* the carry out of the low add */
    uint64_t ch = (uint64_t)(c >> 64) + (uint64_t)(lo >> 64) + carry + cross;
    c = ((U)ch << 64) | (U)cl;
}


static inline void pair_dot(int n,const U*a,const U*b,size_t bs,U*c){
 int k=n;
 __asm__ volatile(
  "xor %%r8d,%%r8d\n\txor %%r9d,%%r9d\n\txor %%r10d,%%r10d\n\txor %%r11d,%%r11d\n\t"
  ".p2align 4\n\t1:\n\t"
  "mov 0(%[a]),%%rdx\n\tmov 8(%[a]),%%r12\n\t"
  "mulx 0(%[b]),%%rax,%%rbx\n\t"
  "mov %%r12,%%r13\n\timul 0(%[b]),%%r13\n\t"
  "mov %%rdx,%%r14\n\timul 8(%[b]),%%r14\n\t"
  "add %%r14,%%r13\n\tadd %%rax,%%r8\n\tadc %%rbx,%%r9\n\tadd %%r13,%%r9\n\t"
  "mulx 16(%[b]),%%rax,%%rbx\n\t"
  "imul 16(%[b]),%%r12\n\timul 24(%[b]),%%rdx\n\t"
  "add %%rdx,%%r12\n\tadd %%rax,%%r10\n\tadc %%rbx,%%r11\n\tadd %%r12,%%r11\n\t"
  "add $16,%[a]\n\tadd %[bs],%[b]\n\n\t"
  "mov 0(%[a]),%%rdx\n\tmov 8(%[a]),%%r12\n\t"
  "mulx 0(%[b]),%%rax,%%rbx\n\t"
  "mov %%r12,%%r13\n\timul 0(%[b]),%%r13\n\t"
  "mov %%rdx,%%r14\n\timul 8(%[b]),%%r14\n\t"
  "add %%r14,%%r13\n\tadd %%rax,%%r8\n\tadc %%rbx,%%r9\n\tadd %%r13,%%r9\n\t"
  "mulx 16(%[b]),%%rax,%%rbx\n\t"
  "imul 16(%[b]),%%r12\n\timul 24(%[b]),%%rdx\n\t"
  "add %%rdx,%%r12\n\tadd %%rax,%%r10\n\tadc %%rbx,%%r11\n\tadd %%r12,%%r11\n\t"
  "add $16,%[a]\n\tadd %[bs],%[b]\n\n\t"
  "dec %[k]\n\tjnz 1b\n\t"
  "mov %%r8,0(%[c])\n\tmov %%r9,8(%[c])\n\tmov %%r10,16(%[c])\n\tmov %%r11,24(%[c])\n\t"
  : [a] "+&r"(a),[b] "+&r"(b),[k] "+&r"(k)
  : [bs] "r"(bs),[c] "r"(c)
  : "rax","rbx","rdx","r8","r9","r10","r11","r12","r13","r14","cc","memory");
}
static void classical(int n, const U *a, int sa, const U *b, int sb, U *c, int sc) {
    /* k-inner leaf: two accumulators for adjacent columns stay in registers across the
       whole k loop, so C is loaded/stored once per (i,j) pair instead of once per (i,k). */
    for (int i = 0; i < n; ++i) {
        U *cr = c + (size_t)i * sc;
        const U *ar = a + (size_t)i * sa;
        int j = 0;
        if (n % 2 == 0) {
            for (; j + 2 <= n; j += 2) pair_dot(n>>1,ar,b+j,(size_t)sb*sizeof(U),cr+j);
        } else {
            for (; j + 2 <= n; j += 2) { U x=0,y=0; for(int k=0;k<n;k++){mac(x,ar[k],b[(size_t)k*sb+j]);mac(y,ar[k],b[(size_t)k*sb+j+1]);} cr[j]=x;cr[j+1]=y; }   /* [r41rx-PDOT-UNROLL2] k 二路展开 */
        }
        for (; j < n; ++j) {
            U x = 0;
            for (int k = 0; k < n; ++k) mac(x, ar[k], b[(size_t)k * sb + j]);
            cr[j] = x;
        }
    }
}


template<bool Sub> static inline __m256i transform128(__m256i a,__m256i b){
 const __m256i sign=_mm256_set1_epi64x(INT64_MIN);
 __m256i v=Sub?_mm256_sub_epi64(a,b):_mm256_add_epi64(a,b);
 __m256i mask;
 if constexpr(Sub) mask=_mm256_cmpgt_epi64(_mm256_xor_si256(b,sign),_mm256_xor_si256(a,sign));
 else mask=_mm256_cmpgt_epi64(_mm256_xor_si256(a,sign),_mm256_xor_si256(v,sign));
 mask=_mm256_slli_si256(mask,8);
 return Sub?_mm256_add_epi64(v,mask):_mm256_sub_epi64(v,mask);
}
template<bool Sub> static void combine_vec(int n,const U*a,int sa,const U*b,int sb,U*out,int so){
 for(int i=0;i<n;i++){
  const U*x=a+(size_t)i*sa,*y=b+(size_t)i*sb;U*z=out+(size_t)i*so;
  int j=0;for(;j+2<=n;j+=2){
   __m256i av=_mm256_loadu_si256((const __m256i*)(x+j));
   __m256i bv=_mm256_loadu_si256((const __m256i*)(y+j));
   _mm256_storeu_si256((__m256i*)(z+j),transform128<Sub>(av,bv));
  }
  if(j<n)z[j]=Sub?x[j]-y[j]:x[j]+y[j];
 }
}
static void combine(int n,const U*a,int sa,const U*b,int sb,U*out,int sign){
 if(sign>0)combine_vec<false>(n,a,sa,b,sb,out,n);
 else combine_vec<true>(n,a,sa,b,sb,out,n);
}
static void addto_stride(int n,U*c,int sc,const U*p,int sp,int sign){
 if(sign>0)combine_vec<false>(n,c,sc,p,sp,c,sc);
 else combine_vec<true>(n,c,sc,p,sp,c,sc);
}
static void addto(int n,U*c,int sc,const U*p,int sign){addto_stride(n,c,sc,p,n,sign);}

#define LEAFMAX 32

static void recurse(int n, const U *a, int sa, const U *b, int sb, U *c, int sc, U *work) {
    if (n <= LEAFMAX) { classical(n, a, sa, b, sb, c, sc); return; }
    int h = n / 2;
    size_t count = (size_t)h * h;
    U *s = work, *t = s + count, *p = t + count, *next = p + count;
    const U *a00 = a, *a01 = a + h, *a10 = a + (size_t)h * sa, *a11 = a10 + h;
    const U *b00 = b, *b01 = b + h, *b10 = b + (size_t)h * sb, *b11 = b10 + h;
    U *c00 = c, *c01 = c + h, *c10 = c + (size_t)h * sc, *c11 = c10 + h;
    combine(h, a00, sa, a11, sa, s, 1); combine(h, b00, sb, b11, sb, t, 1);
    recurse(h, s, h, t, h, c00, sc, next);
    for(int i=0;i<h;i++) memcpy(c11+(size_t)i*sc,c00+(size_t)i*sc,(size_t)h*sizeof(U));
    combine(h, a10, sa, a11, sa, s, 1);
    recurse(h, s, h, b00, sb, c10, sc, next);
    addto_stride(h,c11,sc,c10,sc,-1);
    combine(h, b01, sb, b11, sb, t, -1);
    recurse(h, a00, sa, t, h, c01, sc, next);
    addto_stride(h,c11,sc,c01,sc,1);
    combine(h, b10, sb, b00, sb, t, -1);
    recurse(h, a11, sa, t, h, p, h, next); addto(h, c00, sc, p, 1); addto(h, c10, sc, p, 1);
    combine(h, a00, sa, a01, sa, s, 1);
    recurse(h, s, h, b11, sb, p, h, next); addto(h, c00, sc, p, -1); addto(h, c01, sc, p, 1);
    combine(h, a10, sa, a00, sa, s, -1); combine(h, b00, sb, b01, sb, t, 1);
    recurse(h, s, h, t, h, p, h, next); addto(h, c11, sc, p, 1);
    combine(h, a01, sa, a11, sa, s, -1); combine(h, b10, sb, b11, sb, t, 1);
    recurse(h, s, h, t, h, p, h, next); addto(h, c00, sc, p, 1);
}
}  /* namespace matrix_x */

void matrix_multiply(int n, const __int128 *A, const __int128 *B, __int128 *C) {
    using namespace matrix_x;
    if (n <= LEAFMAX || (n & (n - 1))) { classical(n, (const U *)A, n, (const U *)B, n, (U *)C, n); return; }
    U *work = nullptr;
    if (posix_memalign((void **)&work, 4096, (size_t)n * n * sizeof(U))) abort();
    recurse(n, (const U *)A, n, (const U *)B, n, (U *)C, n, work);
    free(work);
}

评测结果