/* References:
- saffah_cc_v41_agg1, https://duck.ac/submission/126863: 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 #125857 <https://duck.ac/submission/125857>
// -- rank 1 on `mmmx2k` (judge 6.589 s; 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 125712 <https://duck.ac/submission/125712>
// -- 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 #117857 <https://duck.ac/submission/117857>
// -- 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_ 席 · mmmx2k 第 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 <stdint.h>
#include <stdlib.h>
#include <string.h>
#pragma GCC optimize("O3")
#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 */
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 */
c += lo + ((U)cross << 64);
}
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;
for (; j + 2 <= n; j += 2) {
U x = 0, y = 0;
const U *bp = b + j;
for (int k = 0; k < n; ++k) {
U av = ar[k];
const U *brk = bp + (size_t)k * sb;
mac(x, av, brk[0]);
mac(y, av, brk[1]);
}
cr[j] = x; cr[j + 1] = y;
}
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;
}
}
}
static void combine(int n, const U *a, int sa, const U *b, int sb, U *out, int sign) {
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
out[(size_t)i * n + j] = sign > 0 ? a[(size_t)i * sa + j] + b[(size_t)i * sb + j]
: a[(size_t)i * sa + j] - b[(size_t)i * sb + j];
}
static void addto(int n, U *c, int sc, const U *p, int sign) {
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j) {
if (sign > 0) c[(size_t)i * sc + j] += p[(size_t)i * n + j];
else c[(size_t)i * sc + j] -= p[(size_t)i * n + j];
}
}
#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);
for(int i=0;i<h;i++)for(int j=0;j<h;j++)c11[(size_t)i*sc+j]-=c10[(size_t)i*sc+j];
combine(h, b01, sb, b11, sb, t, -1);
recurse(h, a00, sa, t, h, c01, sc, next);
for(int i=0;i<h;i++)for(int j=0;j<h;j++)c11[(size_t)i*sc+j]+=c01[(size_t)i*sc+j];
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);
}