代码

/* lane_ch5k ARM bshuf  (base md5 de925e97c9e3)
   EDIT:  */
/* lane_ch4k ARM blk256  (base md5 f3826e5b825e, arm md5 below)
   EDIT: Block the n=256 materialised buffers too (gate 512 -> 256): makes the LEAF's
       packA/packB (and the level-256 dest) dense with ld == 128 == chunk, instead of
       ld == 256 (span 2x).  Record named this 'a pure repair ... price it first'. */
#pragma GCC optimize("O3","unroll-loops")
/* mmml: C = A*B mod 2^64, long long, FUNCTION-style C++ linkage.
   5-ALU-uop mikro (vpmuludq + vpaddq + vpmulld + vpaddd per k x row x 4 cols, no p5 uop)
   plus Strassen with fused operand combination and dual-destination stores
   (no extra memory).  MR=6 KC=512 MC=48 NC=512 CUTL=512 U=1 BSW=0 EXP=0 */
#include <immintrin.h>

/* Single-instruction 256-bit unaligned access, NO alignment precondition.
   g++-9 -O2 with no -march expands the unaligned 256-bit intrinsics into 2-3 uops, one
   of them a p5-only vinsertf128/vextractf128; these emit the single instruction the
   intrinsic is named for.  AT&T order is `vmovdqu src,dst`: the STORE is %0,%1 (ymm
   first, memory second).  A REVERSED store operand silently emits a LOAD instead --
   that is how this was first written and it produced wrong output with no diagnostic.
   volatile + "memory" is REQUIRED: a non-volatile asm with only an "m" input is pure
   and gets eliminated. */
#define LDQ(p) ({ __m256i _ldq_v; __asm__("vmovdqu %1,%0" : "=x"(_ldq_v) : "m"(*(const __m256i *)(p)) : "memory"); _ldq_v; })
#define STQ(p, v) __asm__ volatile("vmovdqu %0,%1" :: "x"(v), "m"(*(__m256i *)(p)) : "memory")
#include <stddef.h>
#pragma GCC push_options
#pragma GCC target("avx2")

typedef long long TE;

#define MR 6
#define KC 128
#define MC 24
#define NC 512
#define UU 1
#define CUTL 128

/* A-panel row stride PADDED by 8 uint64 (64 B).  With the natural stride KC*8 = 4096 B
   every one of the 6 rows of a micro-tile maps to the SAME L1 set (4096/64 = 64 sets
   apart, 64 sets in the cache) -- a textbook conflict.  +8 uint64 makes the stride
   4160 B = 65 lines, so consecutive rows land in consecutive sets. */
#define KCP (KC + 0)
static TE Apanel[(size_t)MC * KCP + 64] __attribute__((aligned(64)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(64)));
/* lane_ch5k: BpanelS is the SAME panel with each 64-bit lane's two 32-bit halves
   swapped -- byte-for-byte what mikro's per-k `vpshufd $0xB1` produced.  Building it
   once per packB call (4096 shuffles/leaf) replaces 86016 per-k shuffles/leaf, and the
   k-loop's binding p015 port count drops from 31 to 30.  No new memory: the Bpanel
   reservation is 2*(NC/4+1)*KC*8 elements and only (nc/4)*KC*4 are used. */
static TE BpanelS[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(64)));
#define BSOFF ((const char *)BpanelS - (const char *)Bpanel)
static __m256i AC[12] __attribute__((aligned(64)));
static volatile int g_blk = 1;
#define BLKSEL 1   /* 0 = frame control (original path) */


/* op: 0 = store, 1 = add, 2 = sub;  c2 may be 0 (single destination) */
static inline void mikro(long cnt, const long long *ap, const long long *bp,
                        long long *c1, int ldc1, int op1,
                        long long *c2, int ldc2, int op2, int rows) {
  long idx = -(cnt << 3);
  long long *apE = (long long *)ap + cnt;
  long long *bpE = (long long *)bp + (cnt << 2);
    __asm__ volatile(
    "vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
    "vmovdqa %%ymm0, %%ymm1\n\t"
    "vmovdqa %%ymm0, %%ymm2\n\t"
    "vmovdqa %%ymm0, %%ymm3\n\t"
    "vmovdqa %%ymm0, %%ymm4\n\t"
    "vmovdqa %%ymm0, %%ymm5\n\t"
    "vmovdqa %%ymm0, %%ymm6\n\t"
    "vmovdqa %%ymm0, %%ymm7\n\t"
    "vmovdqa %%ymm0, %%ymm8\n\t"
    "vmovdqa %%ymm0, %%ymm9\n\t"
    "vmovdqa %%ymm0, %%ymm10\n\t"
    "vmovdqa %%ymm0, %%ymm11\n\t"
    "test %[idx], %[idx]\n\t"
    "jz 2f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu (%[bs],%[idx],4), %%ymm13\n\t"
    ".balign 4096\n\t"
    "nopl 0x0\n\t"
    "nop\n\t"
    "1:\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 1024(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
    "vbroadcastsd 2048(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
    "vbroadcastsd 3072(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
    "vbroadcastsd 4096(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm9, %%ymm9\n\t"
    "vbroadcastsd 5120(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
    ".byte 0x0f,0x1f,0x80,0x00,0x00,0x00,0x00\n\t"
    "vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 32(%[bs],%[idx],4), %%ymm13\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "2:\n\t"
    "vmovdqa %%ymm0, 0(%[ac])\n\t"
    "vmovdqa %%ymm1, 32(%[ac])\n\t"
    "vmovdqa %%ymm2, 64(%[ac])\n\t"
    "vmovdqa %%ymm3, 96(%[ac])\n\t"
    "vmovdqa %%ymm4, 128(%[ac])\n\t"
    "vmovdqa %%ymm5, 160(%[ac])\n\t"
    "vmovdqa %%ymm6, 192(%[ac])\n\t"
    "vmovdqa %%ymm7, 224(%[ac])\n\t"
    "vmovdqa %%ymm8, 256(%[ac])\n\t"
    "vmovdqa %%ymm9, 288(%[ac])\n\t"
    "vmovdqa %%ymm10, 320(%[ac])\n\t"
    "vmovdqa %%ymm11, 352(%[ac])\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
    : [ac] "r"(&AC[0]), [bs] "r"((const long long *)((const char *)bpE + BSOFF))
    : "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
  if (rows > 0) {
    __m256i lo = AC[0], hi = AC[1];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)0 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)0 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
  if (rows > 1) {
    __m256i lo = AC[2], hi = AC[3];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)1 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)1 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
  if (rows > 2) {
    __m256i lo = AC[4], hi = AC[5];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)2 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)2 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
  if (rows > 3) {
    __m256i lo = AC[6], hi = AC[7];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)3 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)3 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
  if (rows > 4) {
    __m256i lo = AC[8], hi = AC[9];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)4 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)4 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
  if (rows > 5) {
    __m256i lo = AC[10], hi = AC[11];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)5 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)5 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
}

/* rows==2 TAIL kernel: the same instruction stream as the 6-row body, truncated to
   two rows (4 accumulators) and the two A-row displacements that group needs. */
static inline void mikro2(long cnt, const long long *ap, const long long *bp,
                        long long *c1, int ldc1, int op1,
                        long long *c2, int ldc2, int op2, int rows) {
  long idx = -(cnt << 3);
  long long *apE = (long long *)ap + cnt;
  long long *bpE = (long long *)bp + (cnt << 2);
    __asm__ volatile(
    "vpxor %%ymm0, %%ymm0, %%ymm0\n\t"
    "vmovdqa %%ymm0, %%ymm1\n\t"
    "vmovdqa %%ymm0, %%ymm2\n\t"
    "vmovdqa %%ymm0, %%ymm3\n\t"
    "test %[idx], %[idx]\n\t"
    "jz 2f\n\t"
    "vmovdqu (%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu (%[bs],%[idx],4), %%ymm13\n\t"
    ".balign 4096\n\t"
    "nopl 0x0\n\t"
    "nop\n\t"
    "1:\n\t"
    "vbroadcastsd 0(%[ap],%[idx],1), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 1024(%[ap],%[idx],1), %%ymm14\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm14\n\t"
    "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vmovdqu 32(%[bp],%[idx],4), %%ymm12\n\t"
    "vmovdqu 32(%[bs],%[idx],4), %%ymm13\n\t"
    "add $8, %[idx]\n\t"
    "jnz 1b\n\t"
    "2:\n\t"
    "vmovdqa %%ymm0, 0(%[ac])\n\t"
    "vmovdqa %%ymm1, 32(%[ac])\n\t"
    "vmovdqa %%ymm2, 64(%[ac])\n\t"
    "vmovdqa %%ymm3, 96(%[ac])\n\t"
    ""
    : [ap] "+r"(apE), [bp] "+r"(bpE), [idx] "+r"(idx)
    : [ac] "r"(&AC[0]), [bs] "r"((const long long *)((const char *)bpE + BSOFF))
    : "ymm0", "ymm1", "ymm2", "ymm3", "ymm12", "ymm13", "ymm14", "ymm15", "cc", "memory");
  if (rows > 0) {
    __m256i lo = AC[0], hi = AC[1];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)0 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)0 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
  if (rows > 1) {
    __m256i lo = AC[2], hi = AC[3];
    __m256i y = _mm256_add_epi32(hi, _mm256_srli_epi64(hi, 32));
    __m256i res = _mm256_add_epi64(lo, _mm256_slli_epi64(y, 32));
    long long *cp = c1 + (size_t)1 * ldc1;
    if (op1 == 0) STQ((__m256i *)cp, res);
    else if (op1 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
    else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)1 * ldc2;
      if (op2 == 0) STQ((__m256i *)cp, res);
      else if (op2 == 1) STQ((__m256i *)cp, _mm256_add_epi64(LDQ((const __m256i *)cp), res));
      else STQ((__m256i *)cp, _mm256_sub_epi64(LDQ((const __m256i *)cp), res)); } }
}


struct Opd {                 /* value(i,j) = sum_t s[t]*p[(r[t]+i)*ld + c[t]+j]  (nt <= 8) */
  const long long *p; int ld; int nt; int blk; int r[16], c[16], s[16];
};
/* m155zb: a quadrant of a BLOCKED buffer IS a dense (n/2)x(n/2) sub-array. */
static inline Opd blk(const Opd *o, int k, int h) {
  Opd d = *o; d.p = o->p + (size_t)k*h*h; d.ld = h; d.nt = 1;
  d.r[0] = 0; d.c[0] = 0; d.s[0] = 1; d.blk = 0; return d;
}
static inline Opd oq(const Opd *o, int qi, int qj, int h) {
  Opd d = *o;
  for (int t = 0; t < o->nt; t++) { d.r[t] += qi * h; d.c[t] += qj * h; }
  return d;
}
static inline Opd oadd(const Opd *a, const Opd *b, int sgn) {
  /* m155zb: the two operands may live in DIFFERENT blocks (different base pointers), so the
     offset of b's terms is folded into a's addressing frame: off = (b.p - a.p) + r*ld_b + c. */
  Opd d = *a;
  for (int t = 0; t < b->nt; t++) {
    const long long off = (long long)(b->p - d.p) + (long long)b->r[t]*b->ld + b->c[t];
    d.r[d.nt] = (int)(off / d.ld); d.c[d.nt] = (int)(off % d.ld);
    d.s[d.nt] = b->s[t] < 0 ? -sgn : sgn; d.nt++;
  }
  return d;
}
static inline Opd o1(const long long *p, int ld) {
  Opd d; d.p = p; d.ld = ld; d.nt = 1; d.r[0] = 0; d.c[0] = 0; d.s[0] = 1; d.blk = 0; return d;
}

/* A panel: rows [ic, ic+mc), k in [pc, pc+kc) of the operand */
static void packA(const Opd *o, int ic, int mc, int pc, int kc) {
  const int nt = o->nt;
  const long long *P = o->p; const int ld = o->ld;
  for (int r = 0; r < MC; r++) {
    long long *p = Apanel + (size_t)r * KCP;
    int k = 0;
    if (r < mc) {
      const long long *s0 = P + (size_t)(o->r[0] + ic + r) * ld + o->c[0] + pc;
      if (nt == 1 || nt == 2) {
        const long long *s1 = (nt == 2) ? P + (size_t)(o->r[1] + ic + r) * ld + o->c[1] + pc : 0;
        int sgn1 = (nt == 2) ? o->s[1] : 1;
        for (; k + 4 <= kc; k += 4) {
          __m256i v = LDQ((const __m256i *)(s0 + k));
          if (nt == 2) {
            __m256i w = LDQ((const __m256i *)(s1 + k));
            v = (sgn1 > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
          }
        STQ((__m256i *)p, v);
        p += 4;
        }
        for (; k < kc; k++) { long long v = s0[k];
          if (nt == 2) v = (sgn1 > 0) ? v + s1[k] : v - s1[k];
          *p++ = (v); }
      } else {
        for (; k + 4 <= kc; k += 4) {
          __m256i v = _mm256_setzero_si256();
          for (int t = 0; t < nt; t++) {
            __m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k));
            v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
          }
          STQ((__m256i *)p, v);
          p += 4;
        }
        for (; k < kc; k++) {
          long long v = 0;
          for (int t = 0; t < nt; t++)
            v += o->s[t] < 0 ? -P[(size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k]
                             :  P[(size_t)(o->r[t] + ic + r) * ld + o->c[t] + pc + k];
          *p++ = (v);
        }
      }
    }
    p += KC - k;   /* k >= kc is never read: mikro runs kc steps */
  }
}

static void packB(const Opd *o, int jc, int nc, int pc, int kc) {
  const int nt = o->nt;
  const long long *P = o->p; const int ld = o->ld;
  const int jbn = (nc + 3) & ~3;   /* mikro reads only jg < nc, step 4 */
  int jb = 0;
  /* PACK8: 8 columns per k.  The 4-column form loads 32 B at a stride of
     ld*8 (4 KB or more), so every load is a distinct, half-used line and the
     hardware prefetcher cannot see the stream.  Reading jb and jb+4 together
     makes each k one FULL 64 B line, halving the line count and the DRAM
     traffic.  All 8 columns come from one line because jc+jb is a multiple of
     8 (jb steps by 8) and P is 4096-byte aligned.  Byte-identical: the two
     stores write exactly the values the 4-column loop wrote to the same
     addresses. */
  if (nt <= 2) {
    for (; jb + 8 <= nc; jb += 8) {
      long long *q0 = Bpanel + (size_t)(jb / 4) * KC * 4;
      long long *q1 = q0 + (size_t)KC * 4;
      const long long *sq  = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
      const long long *sq1 = (nt >= 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
      int k = 0;
      /* 4 rows in flight: the panel write is k-major while the source is row-major, so
         each row is a separate 64 B line and the load-side MLP is what binds. */
      for (; k + 4 <= kc; k += 4) {
        const long long *a0 = sq + (size_t)(pc + k) * ld;
        const long long *a1 = a0 + ld, *a2 = a1 + ld, *a3 = a2 + ld;
        __m256i v0 = LDQ((const __m256i *)a0), u0 = LDQ((const __m256i *)(a0 + 4));
        __m256i v1 = LDQ((const __m256i *)a1), u1 = LDQ((const __m256i *)(a1 + 4));
        __m256i v2 = LDQ((const __m256i *)a2), u2 = LDQ((const __m256i *)(a2 + 4));
        __m256i v3 = LDQ((const __m256i *)a3), u3 = LDQ((const __m256i *)(a3 + 4));
        if (nt == 2) {
          const long long *b0 = sq1 + (size_t)(pc + k) * ld;
          const long long *b1 = b0 + ld, *b2 = b1 + ld, *b3 = b2 + ld;
          __m256i w0 = LDQ((const __m256i *)b0), x0 = LDQ((const __m256i *)(b0 + 4));
          __m256i w1 = LDQ((const __m256i *)b1), x1 = LDQ((const __m256i *)(b1 + 4));
          __m256i w2 = LDQ((const __m256i *)b2), x2 = LDQ((const __m256i *)(b2 + 4));
          __m256i w3 = LDQ((const __m256i *)b3), x3 = LDQ((const __m256i *)(b3 + 4));
          if (o->s[1] > 0) {
            v0 = _mm256_add_epi64(v0, w0); u0 = _mm256_add_epi64(u0, x0);
            v1 = _mm256_add_epi64(v1, w1); u1 = _mm256_add_epi64(u1, x1);
            v2 = _mm256_add_epi64(v2, w2); u2 = _mm256_add_epi64(u2, x2);
            v3 = _mm256_add_epi64(v3, w3); u3 = _mm256_add_epi64(u3, x3);
          } else {
            v0 = _mm256_sub_epi64(v0, w0); u0 = _mm256_sub_epi64(u0, x0);
            v1 = _mm256_sub_epi64(v1, w1); u1 = _mm256_sub_epi64(u1, x1);
            v2 = _mm256_sub_epi64(v2, w2); u2 = _mm256_sub_epi64(u2, x2);
            v3 = _mm256_sub_epi64(v3, w3); u3 = _mm256_sub_epi64(u3, x3);
          }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), u0);
        _mm256_store_si256((__m256i *)(q0 + (size_t)(k + 1) * 4), v1);
        _mm256_store_si256((__m256i *)(q1 + (size_t)(k + 1) * 4), u1);
        _mm256_store_si256((__m256i *)(q0 + (size_t)(k + 2) * 4), v2);
        _mm256_store_si256((__m256i *)(q1 + (size_t)(k + 2) * 4), u2);
        _mm256_store_si256((__m256i *)(q0 + (size_t)(k + 3) * 4), v3);
        _mm256_store_si256((__m256i *)(q1 + (size_t)(k + 3) * 4), u3);
      }
      for (; k + 2 <= kc; k += 2) {
        const long long *a0 = sq + (size_t)(pc + k) * ld;
        const long long *a1 = a0 + ld;
        __m256i v0 = LDQ((const __m256i *)a0), u0 = LDQ((const __m256i *)(a0 + 4));
        __m256i v1 = LDQ((const __m256i *)a1), u1 = LDQ((const __m256i *)(a1 + 4));
        if (nt == 2) {
          const long long *b0 = sq1 + (size_t)(pc + k) * ld;
          const long long *b1 = b0 + ld;
          __m256i w0 = LDQ((const __m256i *)b0), x0 = LDQ((const __m256i *)(b0 + 4));
          __m256i w1 = LDQ((const __m256i *)b1), x1 = LDQ((const __m256i *)(b1 + 4));
          if (o->s[1] > 0) { v0 = _mm256_add_epi64(v0, w0); u0 = _mm256_add_epi64(u0, x0);
                             v1 = _mm256_add_epi64(v1, w1); u1 = _mm256_add_epi64(u1, x1); }
          else             { v0 = _mm256_sub_epi64(v0, w0); u0 = _mm256_sub_epi64(u0, x0);
                             v1 = _mm256_sub_epi64(v1, w1); u1 = _mm256_sub_epi64(u1, x1); }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v0);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), u0);
        _mm256_store_si256((__m256i *)(q0 + (size_t)(k + 1) * 4), v1);
        _mm256_store_si256((__m256i *)(q1 + (size_t)(k + 1) * 4), u1);
      }
      for (; k < kc; k++) {
        const long long *a = sq + (size_t)(pc + k) * ld;
        __m256i v = LDQ((const __m256i *)a);
        __m256i u = LDQ((const __m256i *)(a + 4));
        if (nt == 2) {
          const long long *b = sq1 + (size_t)(pc + k) * ld;
          __m256i w = LDQ((const __m256i *)b);
          __m256i x = LDQ((const __m256i *)(b + 4));
          if (o->s[1] > 0) { v = _mm256_add_epi64(v, w); u = _mm256_add_epi64(u, x); }
          else             { v = _mm256_sub_epi64(v, w); u = _mm256_sub_epi64(u, x); }
        }
        _mm256_store_si256((__m256i *)(q0 + (size_t)k * 4), v);
        _mm256_store_si256((__m256i *)(q1 + (size_t)k * 4), u);
      }
    }
  }
  for (; jb < jbn; jb += 4) {
    long long *q = Bpanel + (size_t)(jb / 4) * KC * 4;
    int full = (jb + 4 <= nc);
    int lim = nc - jb; if (lim > 4) lim = 4; if (lim < 0) lim = 0;
    const long long *sq = P + (size_t)o->r[0] * ld + o->c[0] + jc + jb;
    const long long *sq1 = (nt >= 2) ? P + (size_t)o->r[1] * ld + o->c[1] + jc + jb : 0;
    if (full && nt <= 2) {
      for (int k = 0; k < kc; k++) {
        __m256i v = LDQ((const __m256i *)(sq + (size_t)(pc + k) * ld));
        if (nt == 2) {
          __m256i w = LDQ((const __m256i *)(sq1 + (size_t)(pc + k) * ld));
          v = (o->s[1] > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
        }
        _mm256_store_si256((__m256i *)(q + (size_t)k * 4), v);
      }
    } else if (full) {
      for (int k = 0; k < kc; k++) {
        __m256i v = _mm256_setzero_si256();
        for (int tt = 0; tt < nt; tt++) {
          __m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb));
          v = (o->s[tt] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
        }
        _mm256_store_si256((__m256i *)(q + (size_t)k * 4), v);
      }
    } else {
      for (int k = 0; k < KC; k++) {
        long long t[4] = {0, 0, 0, 0};
        if (k < kc) {
          int ncol = full ? 4 : lim;
          for (int c = 0; c < ncol; c++) {
            long long v = 0;
            for (int tt = 0; tt < nt; tt++)
              v += o->s[tt] < 0 ? -P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c]
                                :  P[(size_t)(o->r[tt] + pc + k) * ld + o->c[tt] + jc + jb + c];
            t[c] = v;
          }
        }
        __m256i v = LDQ((const __m256i *)t);
        _mm256_store_si256((__m256i *)(q + (size_t)k * 4), v);
      }
    }
  }
  /* lane_ch5k: duplicate the panel with every 64-bit lane's halves swapped, so
     the mikro needs no vpshufd in its k-loop.  Same volume as the pack itself. */
  {
    long long *d = BpanelS;
    const long long *sp = Bpanel;
    const size_t nv = (size_t)(jbn / 4) * KC * 4;
    for (size_t i = 0; i < nv; i += 4)
      _mm256_store_si256((__m256i *)(d + i),
                         _mm256_shuffle_epi32(LDQ((const __m256i *)(sp + i)), 0xB1));
  }
}


static inline int effop(int base, int term) {   /* term: 1 = +M, 2 = -M ; result: 0 store 1 add 2 sub */
  if (base == 0) return term == 1 ? 0 : 2;
  if (base == 1) return term == 2 ? 2 : 1;
  return term == 2 ? 1 : 2;
}

/* C1 op1 (+- second destination C2 op2, either may be absent) = Ao * Bo, block size n */
static void leaf(int n, const Opd *Ao, const Opd *Bo,
                 long long *C1, int ldc1, int op1, long long *C2, int ldc2, int op2) {
  int nb4 = n - (n & 3);
  for (int jc = 0; jc < nb4; jc += NC) {
    int nc = (nb4 - jc < NC) ? (nb4 - jc) : NC;
    for (int pc = 0; pc < n; pc += KC) {
      int kc = (n - pc < KC) ? (n - pc) : KC;
      packB(Bo, jc, nc, pc, kc);
      int e1 = (pc == 0) ? op1 : (op1 == 2 ? 2 : 1);
      int e2 = (pc == 0) ? op2 : (op2 == 2 ? 2 : 1);
      for (int ic = 0; ic < n; ic += MC) {
        int mc = (n - ic < MC) ? (n - ic) : MC;
        packA(Ao, ic, mc, pc, kc);
        for (int jgb = 0; jgb < nc; jgb += 64) {
          int jge = (nc - jgb < 64) ? nc : jgb + 64;
          for (int ir = 0; ir < mc; ir += MR) {
            int rows = (mc - ir < MR) ? (mc - ir) : MR;
            const long long *ap = Apanel + (size_t)ir * KCP * 1;
            for (int jg = jgb; jg < jge; jg += 4) {
              const long long *bp = Bpanel + (size_t)(jg / 4) * KC * 4;
              if (rows == 2) {                 /* TAIL: 2 valid rows, 2-row kernel */
                mikro2(kc / UU, ap, bp, C1 + (size_t)(ic + ir) * ldc1 + jc + jg, ldc1, e1,
                       C2 ? C2 + (size_t)(ic + ir) * ldc2 + jc + jg : 0, ldc2, e2, rows);
              } else {
                mikro(kc / UU, ap, bp, C1 + (size_t)(ic + ir) * ldc1 + jc + jg, ldc1, e1,
                      C2 ? C2 + (size_t)(ic + ir) * ldc2 + jc + jg : 0, ldc2, e2, rows);
              }
            }
          }
        }
      }
    }
  }
  for (int j = nb4; j < n; j++)
    for (int i = 0; i < n; i++) {
      unsigned long long acc = 0;
      for (int k = 0; k < n; k++) {
        long long av = 0, bv = 0;
        for (int t = 0; t < Ao->nt; t++)
          av += Ao->s[t] < 0 ? -Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k]
                             :  Ao->p[(size_t)(Ao->r[t] + i) * Ao->ld + Ao->c[t] + k];
        for (int t = 0; t < Bo->nt; t++)
          bv += Bo->s[t] < 0 ? -Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j]
                             :  Bo->p[(size_t)(Bo->r[t] + k) * Bo->ld + Bo->c[t] + j];
        acc += (unsigned long long)av * (unsigned long long)bv;
      }
      long long *cp = C1 + (size_t)i * ldc1 + j;
      if (op1 == 0) *cp = (long long)acc;
      else if (op1 == 1) *cp += (long long)acc;
      else *cp -= (long long)acc;
      if (C2) {
        long long *cq = C2 + (size_t)i * ldc2 + j;
        if (op2 == 0) *cq = (long long)acc;
        else if (op2 == 1) *cq += (long long)acc;
        else *cq -= (long long)acc;
      }
    }
}

/* recursive Strassen: C1 op1 op= Ao*Bo  and  C2 op2 op= Ao*Bo */
/* Scratch: one h*h product block per recursion level, carved from the tail of the
   buffer.  Only the levels that actually recurse consume space, so a walk of the
   recursion computes exactly the bytes matrix_multiply() hands us below. */
static long long g_scratch[6u * 1024u * 1024u] __attribute__((aligned(64)));

/* dst op= sign * T   (op: 0 store, 1 add, 2 sub) */
static void cmb(long long *dst, int ldd, const long long *T, int ldt, int h, int sign, int op) {
  __m256i sgn = _mm256_set1_epi64x(sign < 0 ? -1 : 0);
  for (int i = 0; i < h; i++) {
    const long long *t = T + (size_t)i * ldt;
    long long *c = dst + (size_t)i * ldd;
    int j = 0;
    if (op == 0) {
      for (; j + 4 <= h; j += 4) {
        __m256i v = LDQ((const __m256i *)(t + j));
        if (sign < 0) v = _mm256_sub_epi64(_mm256_setzero_si256(), v);
        STQ((__m256i *)(c + j), v);
      }
      for (; j < h; j++) c[j] = sign < 0 ? -t[j] : t[j];
    } else {
      for (; j + 4 <= h; j += 4) {
        __m256i v = LDQ((const __m256i *)(t + j));
        __m256i w = LDQ((const __m256i *)(c + j));
        if ((op == 1) != (sign < 0)) STQ((__m256i *)(c + j), _mm256_add_epi64(w, v));
        else                         STQ((__m256i *)(c + j), _mm256_sub_epi64(w, v));
      }
      for (; j < h; j++) { if ((op == 1) != (sign < 0)) c[j] += t[j]; else c[j] -= t[j]; }
    }
  }
}

/* fused two-destination combine: reads T ONCE for both destinations.
   Semantics identical to cmb(d1,s1,o1) then cmb(d2,s2,o2) on the same T. */
static void cmb2(long long *d1,int ld1,int op1,int s1,
                 long long *d2,int ld2,int op2,int s2,
                 const long long *T,int ldt,int h) {
  for (int i = 0; i < h; i++) {
    const long long *t = T + (size_t)i * ldt;
    long long *c1 = d1 + (size_t)i * ld1, *c2 = d2 + (size_t)i * ld2;
    int j = 0;
    for (; j + 16 <= h; j += 16) {
      __m256i v0=LDQ((const __m256i*)(t+j)),   v1=LDQ((const __m256i*)(t+j+4));
      __m256i v2=LDQ((const __m256i*)(t+j+8)), v3=LDQ((const __m256i*)(t+j+12));
      __m256i a0=LDQ((const __m256i*)(c1+j)),   a1=LDQ((const __m256i*)(c1+j+4));
      __m256i a2=LDQ((const __m256i*)(c1+j+8)), a3=LDQ((const __m256i*)(c1+j+12));
      __m256i b0=LDQ((const __m256i*)(c2+j)),   b1=LDQ((const __m256i*)(c2+j+4));
      __m256i b2=LDQ((const __m256i*)(c2+j+8)), b3=LDQ((const __m256i*)(c2+j+12));
      if ((op1 == 1) != (s1 < 0)) { STQ((__m256i*)(c1+j),_mm256_add_epi64(a0,v0));
        STQ((__m256i*)(c1+j+4),_mm256_add_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_add_epi64(a2,v2));
        STQ((__m256i*)(c1+j+12),_mm256_add_epi64(a3,v3)); }
      else { STQ((__m256i*)(c1+j),_mm256_sub_epi64(a0,v0));
        STQ((__m256i*)(c1+j+4),_mm256_sub_epi64(a1,v1)); STQ((__m256i*)(c1+j+8),_mm256_sub_epi64(a2,v2));
        STQ((__m256i*)(c1+j+12),_mm256_sub_epi64(a3,v3)); }
      if ((op2 == 1) != (s2 < 0)) { STQ((__m256i*)(c2+j),_mm256_add_epi64(b0,v0));
        STQ((__m256i*)(c2+j+4),_mm256_add_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_add_epi64(b2,v2));
        STQ((__m256i*)(c2+j+12),_mm256_add_epi64(b3,v3)); }
      else { STQ((__m256i*)(c2+j),_mm256_sub_epi64(b0,v0));
        STQ((__m256i*)(c2+j+4),_mm256_sub_epi64(b1,v1)); STQ((__m256i*)(c2+j+8),_mm256_sub_epi64(b2,v2));
        STQ((__m256i*)(c2+j+12),_mm256_sub_epi64(b3,v3)); }
    }
    for (; j + 4 <= h; j += 4) {
      __m256i v = LDQ((const __m256i *)(t + j));
      __m256i a = LDQ((const __m256i *)(c1 + j));
      STQ((__m256i *)(c1 + j), ((op1 == 1) != (s1 < 0)) ? _mm256_add_epi64(a,v)
                                                        : _mm256_sub_epi64(a,v));
      __m256i b = LDQ((const __m256i *)(c2 + j));
      STQ((__m256i *)(c2 + j), ((op2 == 1) != (s2 < 0)) ? _mm256_add_epi64(b,v)
                                                        : _mm256_sub_epi64(b,v));
    }
    for (; j < h; j++) {
      long long v = t[j];
      if ((op1 == 1) != (s1 < 0)) c1[j] += v; else c1[j] -= v;
      if ((op2 == 1) != (s2 < 0)) c2[j] += v; else c2[j] -= v;
    }
  }
}

/* single-destination recursive Strassen: d op= Ao*Bo, with T a scratch block of h*h */
/* materialise the operands of the level whose CHILDREN are leaves, so the leaf's packing
   reads <=2 contiguous sources per element instead of up to 8 strided ones.  Priced in
   work/lane_ax1_mmml123/LEDGER.md: this is what makes the third Strassen level pay. */
static long long matbufA[4][(size_t)256*CUTL*CUTL + 64], matbufB[4][(size_t)256*CUTL*CUTL + 64];

static void matop_blk(const Opd *o, int n, long long *dst) {
  const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
  const int h2 = n >> 1; const size_t hn = (size_t)h2*h2;
  for (int ii = 0; ii + 2 <= n; ii += 2) {
    const int bi = ii / h2, ir = ii % h2;
    /* m171: quadrant slot = index XOR 1 (sigma = [1,0,3,2], the lex-first optimum of the
       6-sum span problem: 12 -> 8, adjacent sums 2 -> 4).  Writer and readers share it. */
    long long *A0 = dst + (size_t)((bi*2)^1)*hn + (size_t)ir*h2;
    long long *B0 = A0 + h2;
    long long *A1 = dst + (size_t)((bi*2+1)^1)*hn + (size_t)ir*h2;
    long long *B1 = A1 + h2;
    const long long *b0 = P + (size_t)o->r[0]*ld + o->c[0] + (size_t)ii*ld;
    const long long *b1 = b0 + ld;
    for (int half = 0; half < 2; half++) {
      const int soff = half * h2;
      long long *w0 = half ? A1 : A0, *w1 = half ? B1 : B0;
      for (int j = 0; j + 16 <= h2; j += 16) {
        __m256i a0=_mm256_setzero_si256(),a1=a0,a2=a0,a3=a0;
        __m256i e0=a0,e1=a0,e2=a0,e3=a0;
        for (int t = 0; t < nt; t++) {
          const long long *p0 = b0 + (size_t)o->r[t]*ld + o->c[t] - (size_t)o->r[0]*ld - o->c[0] + j + soff;
          const long long *p1 = p0 + ld;
          __m256i x0=LDQ((const __m256i*)p0), x1=LDQ((const __m256i*)(p0+4));
          __m256i x2=LDQ((const __m256i*)(p0+8)), x3=LDQ((const __m256i*)(p0+12));
          __m256i y0=LDQ((const __m256i*)p1), y1=LDQ((const __m256i*)(p1+4));
          __m256i y2=LDQ((const __m256i*)(p1+8)), y3=LDQ((const __m256i*)(p1+12));
          if (o->s[t] < 0) { a0=_mm256_sub_epi64(a0,x0); a1=_mm256_sub_epi64(a1,x1);
                             a2=_mm256_sub_epi64(a2,x2); a3=_mm256_sub_epi64(a3,x3);
                             e0=_mm256_sub_epi64(e0,y0); e1=_mm256_sub_epi64(e1,y1);
                             e2=_mm256_sub_epi64(e2,y2); e3=_mm256_sub_epi64(e3,y3); }
          else             { a0=_mm256_add_epi64(a0,x0); a1=_mm256_add_epi64(a1,x1);
                             a2=_mm256_add_epi64(a2,x2); a3=_mm256_add_epi64(a3,x3);
                             e0=_mm256_add_epi64(e0,y0); e1=_mm256_add_epi64(e1,y1);
                             e2=_mm256_add_epi64(e2,y2); e3=_mm256_add_epi64(e3,y3); }
        }
        STQ((__m256i*)(w0+j),a0); STQ((__m256i*)(w0+j+4),a1);
        STQ((__m256i*)(w0+j+8),a2); STQ((__m256i*)(w0+j+12),a3);
        STQ((__m256i*)(w1+j),e0); STQ((__m256i*)(w1+j+4),e1);
        STQ((__m256i*)(w1+j+8),e2); STQ((__m256i*)(w1+j+12),e3);
      }
    }
  }
}

static void matop(const Opd *o, int n, long long *dst) {
  const int nt = o->nt; const long long *P = o->p; const int ld = o->ld;
  int i = 0;
  for (; i + 2 <= n; i += 2) {
    long long *d0 = dst + (size_t)i * n, *d1 = d0 + n;
    const long long *b0 = P + (size_t)o->r[0]*ld + o->c[0] + (size_t)i*ld;
    const long long *b1 = b0 + ld;
    int j = 0;
    for (; j + 16 <= n; j += 16) {
      __m256i a0=_mm256_setzero_si256(),a1=a0,a2=a0,a3=a0;
      __m256i e0=a0,e1=a0,e2=a0,e3=a0;
      for (int t = 0; t < nt; t++) {
        const long long *p0 = b0 + (size_t)o->r[t]*ld + o->c[t] - (size_t)o->r[0]*ld - o->c[0] + j;
        const long long *p1 = p0 + ld;
        __m256i x0=LDQ((const __m256i*)p0), x1=LDQ((const __m256i*)(p0+4));
        __m256i x2=LDQ((const __m256i*)(p0+8)), x3=LDQ((const __m256i*)(p0+12));
        __m256i y0=LDQ((const __m256i*)p1), y1=LDQ((const __m256i*)(p1+4));
        __m256i y2=LDQ((const __m256i*)(p1+8)), y3=LDQ((const __m256i*)(p1+12));
        if (o->s[t] < 0) { a0=_mm256_sub_epi64(a0,x0); a1=_mm256_sub_epi64(a1,x1);
                           a2=_mm256_sub_epi64(a2,x2); a3=_mm256_sub_epi64(a3,x3);
                           e0=_mm256_sub_epi64(e0,y0); e1=_mm256_sub_epi64(e1,y1);
                           e2=_mm256_sub_epi64(e2,y2); e3=_mm256_sub_epi64(e3,y3); }
        else             { a0=_mm256_add_epi64(a0,x0); a1=_mm256_add_epi64(a1,x1);
                           a2=_mm256_add_epi64(a2,x2); a3=_mm256_add_epi64(a3,x3);
                           e0=_mm256_add_epi64(e0,y0); e1=_mm256_add_epi64(e1,y1);
                           e2=_mm256_add_epi64(e2,y2); e3=_mm256_add_epi64(e3,y3); }
      }
      STQ((__m256i*)(d0+j),a0); STQ((__m256i*)(d0+j+4),a1);
      STQ((__m256i*)(d0+j+8),a2); STQ((__m256i*)(d0+j+12),a3);
      STQ((__m256i*)(d1+j),e0); STQ((__m256i*)(d1+j+4),e1);
      STQ((__m256i*)(d1+j+8),e2); STQ((__m256i*)(d1+j+12),e3);
    }
    for (; j + 4 <= n; j += 4) {          /* FIX: was `j < n`, which overran */
      for (int r = 0; r < 2; r++) {
        __m256i v = _mm256_setzero_si256();
        long long *dp = r ? d1 : d0;
        for (int t = 0; t < nt; t++) {
          __m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + i + r) * ld + o->c[t] + j));
          v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
        }
        STQ((__m256i *)(dp + j), v);
      }
    }
    for (; j < n; j++) {                  /* ragged columns handled scalar-ly */
      for (int r = 0; r < 2; r++) {
        unsigned long long v = 0;
        long long *dp = r ? d1 : d0;
        for (int t = 0; t < nt; t++) {
          unsigned long long w =
              (unsigned long long)P[(size_t)(o->r[t] + i + r) * ld + o->c[t] + j];
          v += (o->s[t] < 0) ? (0ULL - w) : w;
        }
        dp[j] = (long long)v;
      }
    }
  }
  for (; i < n; i++) {
    long long *dp = dst + (size_t)i * n;
    int j = 0;
    for (; j + 16 <= n; j += 16) {
      __m256i v0=_mm256_setzero_si256(),v1=v0,v2=v0,v3=v0;
      for (int t = 0; t < nt; t++) {
        const long long *sp = P + (size_t)(o->r[t] + i) * ld + o->c[t] + j;
        __m256i w0=LDQ((const __m256i *)sp), w1=LDQ((const __m256i *)(sp+4));
        __m256i w2=LDQ((const __m256i *)(sp+8)), w3=LDQ((const __m256i *)(sp+12));
        if (o->s[t] < 0) { v0=_mm256_sub_epi64(v0,w0); v1=_mm256_sub_epi64(v1,w1);
                           v2=_mm256_sub_epi64(v2,w2); v3=_mm256_sub_epi64(v3,w3); }
        else             { v0=_mm256_add_epi64(v0,w0); v1=_mm256_add_epi64(v1,w1);
                           v2=_mm256_add_epi64(v2,w2); v3=_mm256_add_epi64(v3,w3); }
      }
      STQ((__m256i *)(dp+j), v0); STQ((__m256i *)(dp+j+4), v1);
      STQ((__m256i *)(dp+j+8), v2); STQ((__m256i *)(dp+j+12), v3);
    }
    for (; j + 4 <= n; j += 4) {          /* FIX: was `j < n`, which overran */
      __m256i v = _mm256_setzero_si256();
      for (int t = 0; t < nt; t++) {
        __m256i w = LDQ((const __m256i *)(P + (size_t)(o->r[t] + i) * ld + o->c[t] + j));
        v = (o->s[t] < 0) ? _mm256_sub_epi64(v, w) : _mm256_add_epi64(v, w);
      }
      STQ((__m256i *)(dp + j), v);
    }
    for (; j < n; j++) {                  /* ragged columns handled scalar-ly */
      unsigned long long v = 0;
      for (int t = 0; t < nt; t++) {
        unsigned long long w = (unsigned long long)P[(size_t)(o->r[t] + i) * ld + o->c[t] + j];
        v += (o->s[t] < 0) ? (0ULL - w) : w;
      }
      dp[j] = (long long)v;
    }
  }
}
static void mmxR(int n, const Opd *Ao0, const Opd *Bo0, long long *d, int ld, int op,
                 long long *T, int lvl) {
  const Opd *Ao = Ao0, *Bo = Bo0; Opd A1, B1;
  if (n <= CUTL || (n & 1) || Ao->nt > 8 || Bo->nt > 8) {
    leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
  if (n <= 16 * CUTL && (Ao->nt > 1 || Bo->nt > 1)) {
    int L = (n > 8 * CUTL) ? 0 : (n > 4 * CUTL) ? 1 : (n > 2 * CUTL) ? 2 : 3;
    if (Ao->nt > 1) { if (BLKSEL && n >= 256 && ((n & 15) == 0)) matop_blk(Ao, n, matbufA[L]); else matop(Ao, n, matbufA[L]); A1 = o1(matbufA[L], n); A1.blk = (BLKSEL && n >= 256 && ((n & 15) == 0)); Ao = &A1; }
    if (Bo->nt > 1) { if (BLKSEL && n >= 256 && ((n & 15) == 0)) matop_blk(Bo, n, matbufB[L]); else matop(Bo, n, matbufB[L]); B1 = o1(matbufB[L], n); B1.blk = (BLKSEL && n >= 256 && ((n & 15) == 0)); Bo = &B1; }
  }
  if (n <= CUTL || (n & 1) || Ao->nt > 8 || Bo->nt > 8) {
    leaf(n, Ao, Bo, d, ld, op, 0, 0, 0); return; }
  const int h = n >> 1;
  Opd A11, A12, A21, A22, B11, B12, B21, B22;
  if (Ao->blk) { A11 = blk(Ao,1,h); A12 = blk(Ao,0,h); A21 = blk(Ao,3,h); A22 = blk(Ao,2,h); }
  else         { A11 = oq(Ao,0,0,h); A12 = oq(Ao,0,1,h); A21 = oq(Ao,1,0,h); A22 = oq(Ao,1,1,h); }
  if (Bo->blk) { B11 = blk(Bo,1,h); B12 = blk(Bo,0,h); B21 = blk(Bo,3,h); B22 = blk(Bo,2,h); }
  else         { B11 = oq(Bo,0,0,h); B12 = oq(Bo,0,1,h); B21 = oq(Bo,1,0,h); B22 = oq(Bo,1,1,h); }
  long long *Q[4];
  Q[0] = d; Q[1] = d + h; Q[2] = d + (size_t)h * ld; Q[3] = Q[2] + h;
  long long *T2 = T + (size_t)h * h;
  int wr[4] = {0, 0, 0, 0};
#define MK(qi, si, qj, sj, Ma, Mb) do { \
    int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
    int ob = wr[qj] ? effop(1, (sj)) : effop(op, (sj)); \
    wr[qi] = 1; wr[qj] = 1; \
    if (oa == 0) { \
      mmxR(h, &(Ma), &(Mb), Q[qi], ld, 0, T, lvl + 1); \
      cmb(Q[qj], ld, Q[qi], ld, h, (sj), ob); \
    } else { \
      mmxR(h, &(Ma), &(Mb), T, h, 0, T2, lvl + 1); \
      cmb2(Q[qi], ld, oa, (si), Q[qj], ld, ob, (sj), T, h, h); \
    } \
  } while (0) \

#define MK1(qi, si, Ma, Mb) do { \
    int oa = wr[qi] ? effop(1, (si)) : effop(op, (si)); \
    wr[qi] = 1; \
    mmxR(h, &(Ma), &(Mb), Q[qi], ld, oa, T, lvl + 1); \
  } while (0)
  { /* M1 = (A11+A22)(B11+B22) -> +C11, +C22 */
    Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
    MK(0, 1, 3, 1, Ma, Mb);
  }
  { /* M2 = (A21+A22)B11 -> +C21, -C22 */
    Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
    MK(2, 1, 3, 2, Ma, Mb);
  }
  { /* M3 = A11(B12-B22) -> +C12, +C22 */
    Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
    MK(1, 1, 3, 1, Ma, Mb);
  }
  { /* M4 = A22(B21-B11) -> +C11, +C21 */
    Opd Ma = A22, Mb = oadd(&B21, &B11, -1);
    MK(0, 1, 2, 1, Ma, Mb);
  }
  { /* M5 = (A11+A12)B22 -> -C11, +C12 */
    Opd Ma = oadd(&A11, &A12, 1), Mb = B22;
    MK(0, 2, 1, 1, Ma, Mb);
  }
  { /* M6 = (A21-A11)(B11+B12) -> +C22 */
    Opd Ma = oadd(&A21, &A11, -1), Mb = oadd(&B11, &B12, 1);
    MK1(3, 1, Ma, Mb);
  }
  { /* M7 = (A12-A22)(B21+B22) -> +C11 */
    Opd Ma = oadd(&A12, &A22, -1), Mb = oadd(&B21, &B22, 1);
    MK1(0, 1, Ma, Mb);
  }
#undef MK
#undef MK1
}

void matrix_multiply(int n, const long long *A, const long long *B, long long *C) {
  Opd Ao = o1(A, n), Bo = o1(B, n);
  if (n <= CUTL || (n & 1)) { leaf(n, &Ao, &Bo, C, n, 0, 0, 0, 0); return; }
  mmxR(n, &Ao, &Bo, C, n, 0, g_scratch, 0);
}
#pragma GCC pop_options

评测结果