代码

/* 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>
#include <stddef.h>
#pragma GCC push_options
#pragma GCC target("avx2")

typedef long long TE;

#define MR 6
#define KC 512
#define MC 48
#define NC 512
#define UU 1
#define CUTL 512

static TE Apanel[(size_t)MC * KC * 1 + 64] __attribute__((aligned(64)));
static TE Bpanel[(size_t)2 * (NC / 4 + 1) * KC * 8 + 64] __attribute__((aligned(64)));
static __m256i AC[12] __attribute__((aligned(64)));

/* 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) {
  __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 %[cnt], %[cnt]\n\t"
    "jle 2f\n\t"
    "1:\n\t"
    "vmovdqu 0(%[bp]), %%ymm12\n\t"
    "vpshufd $0xB1, %%ymm12, %%ymm13\n\t"
    "vbroadcastsd 0(%[ap]), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpaddq %%ymm15, %%ymm0, %%ymm0\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
    "vbroadcastsd 4096(%[ap]), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpaddq %%ymm15, %%ymm2, %%ymm2\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
    "vbroadcastsd 8192(%[ap]), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpaddq %%ymm15, %%ymm4, %%ymm4\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm5, %%ymm5\n\t"
    "vbroadcastsd 12288(%[ap]), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpaddq %%ymm15, %%ymm6, %%ymm6\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
    "vbroadcastsd 16384(%[ap]), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpaddq %%ymm15, %%ymm8, %%ymm8\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
    "vbroadcastsd 20480(%[ap]), %%ymm14\n\t"
    "vpmuludq %%ymm14, %%ymm12, %%ymm15\n\t"
    "vpaddq %%ymm15, %%ymm10, %%ymm10\n\t"
    "vpmulld %%ymm13, %%ymm14, %%ymm15\n\t"
    "vpaddd %%ymm15, %%ymm11, %%ymm11\n\t"
    "add $8, %[ap]\n\t"
    "add $32, %[bp]\n\t"
    "dec %[cnt]\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"(ap), [bp] "+r"(bp), [cnt] "+r"(cnt)
    : [ac] "r"(&AC[0])
    : "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) _mm256_storeu_si256((__m256i *)cp, res);
    else if (op1 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)0 * ldc2;
      if (op2 == 0) _mm256_storeu_si256((__m256i *)cp, res);
      else if (op2 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
      else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((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) _mm256_storeu_si256((__m256i *)cp, res);
    else if (op1 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)1 * ldc2;
      if (op2 == 0) _mm256_storeu_si256((__m256i *)cp, res);
      else if (op2 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
      else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((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) _mm256_storeu_si256((__m256i *)cp, res);
    else if (op1 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)2 * ldc2;
      if (op2 == 0) _mm256_storeu_si256((__m256i *)cp, res);
      else if (op2 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
      else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((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) _mm256_storeu_si256((__m256i *)cp, res);
    else if (op1 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)3 * ldc2;
      if (op2 == 0) _mm256_storeu_si256((__m256i *)cp, res);
      else if (op2 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
      else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((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) _mm256_storeu_si256((__m256i *)cp, res);
    else if (op1 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)4 * ldc2;
      if (op2 == 0) _mm256_storeu_si256((__m256i *)cp, res);
      else if (op2 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
      else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((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) _mm256_storeu_si256((__m256i *)cp, res);
    else if (op1 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
    if (c2) { cp = c2 + (size_t)5 * ldc2;
      if (op2 == 0) _mm256_storeu_si256((__m256i *)cp, res);
      else if (op2 == 1) _mm256_storeu_si256((__m256i *)cp, _mm256_add_epi64(_mm256_loadu_si256((const __m256i *)cp), res));
      else _mm256_storeu_si256((__m256i *)cp, _mm256_sub_epi64(_mm256_loadu_si256((const __m256i *)cp), res)); } }
}

struct Opd {                 /* value(i,j) = sum_t s[t]*p[(r[t]+i)*ld + c[t]+j]  (nt <= 4) */
  const long long *p; int ld; int nt; int r[4], c[4], s[4];
};
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) {
  Opd d = *a;
  for (int t = 0; t < b->nt; t++) { d.r[d.nt] = b->r[t]; d.c[d.nt] = b->c[t]; 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; 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) {
  long long *p = Apanel;
  const int nt = o->nt;
  const long long *P = o->p; const int ld = o->ld;
  for (int r = 0; r < MC; r++) {
    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 = _mm256_loadu_si256((const __m256i *)(s0 + k));
          if (nt == 2) {
            __m256i w = _mm256_loadu_si256((const __m256i *)(s1 + k));
            v = (sgn1 > 0) ? _mm256_add_epi64(v, w) : _mm256_sub_epi64(v, w);
          }
        _mm256_storeu_si256((__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 < 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);
        }
      }
    }
    for (; k < KC; k++) { *p++ = 0; }
  }
}

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;
  for (int jb = 0; jb < NC; 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 = _mm256_loadu_si256((const __m256i *)(sq + (size_t)(pc + k) * ld));
        if (nt == 2) {
          __m256i w = _mm256_loadu_si256((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);
      }
      for (int k = kc; k < KC; k++) _mm256_store_si256((__m256i *)(q + (size_t)k * 4), _mm256_setzero_si256());
    } 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 = _mm256_loadu_si256((const __m256i *)t);
        _mm256_store_si256((__m256i *)(q + (size_t)k * 4), v);
      }
    }
  }
}


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 ir = 0; ir < mc; ir += MR) {
          int rows = (mc - ir < MR) ? (mc - ir) : MR;
          const long long *ap = Apanel + (size_t)ir * KC * 1;
          for (int jg = 0; jg < nc; jg += 4) {
            const long long *bp = Bpanel + (size_t)(jg / 4) * KC * 4;
            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 */
static void mmx(int n, const Opd *Ao, const Opd *Bo,
                long long *d0, int ld0, int op0, long long *d1, int ld1, int op1) {
  /* exactly ONE Strassen level: stop as soon as an operand is itself a combination or
     the product must reach two destinations (the leaf handles both) */
  if (n <= CUTL || (n & 1) || d1 || Ao->nt > 1 || Bo->nt > 1) {
    leaf(n, Ao, Bo, d0, ld0, op0, d1, ld1, op1); return; }
  const int h = n >> 1;
  Opd A11 = oq(Ao, 0, 0, h), A12 = oq(Ao, 0, 1, h), A21 = oq(Ao, 1, 0, h), A22 = oq(Ao, 1, 1, h);
  Opd 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] = d0; Q[1] = d0 + h; Q[2] = d0 + (size_t)h * ld0; Q[3] = Q[2] + h;
  int wr[4] = {0, 0, 0, 0}, opa[4] = {0, 0, 0, 0};
  (void)opa;

  { /* M1 */
    Opd Ma = oadd(&A11, &A22, 1), Mb = oadd(&B11, &B22, 1);
    opa[0] = wr[0] ? effop(1, 1) : effop(op0, 1);
    opa[3] = wr[3] ? effop(1, 1) : effop(op0, 1);
    wr[0] = 1; wr[3] = 1;
    mmx(h, &Ma, &Mb, Q[0], ld0, opa[0], Q[3], ld0, opa[3]);
  }
  { /* M2 */
    Opd Ma = oadd(&A21, &A22, 1), Mb = B11;
    opa[2] = wr[2] ? effop(1, 1) : effop(op0, 1);
    opa[3] = wr[3] ? effop(1, 2) : effop(op0, 2);
    wr[2] = 1; wr[3] = 1;
    mmx(h, &Ma, &Mb, Q[2], ld0, opa[2], Q[3], ld0, opa[3]);
  }
  { /* M3 */
    Opd Ma = A11, Mb = oadd(&B12, &B22, -1);
    opa[1] = wr[1] ? effop(1, 1) : effop(op0, 1);
    opa[3] = wr[3] ? effop(1, 1) : effop(op0, 1);
    wr[1] = 1; wr[3] = 1;
    mmx(h, &Ma, &Mb, Q[1], ld0, opa[1], Q[3], ld0, opa[3]);
  }
  { /* M4 */
    Opd Ma = A22, Mb = oadd(&B21, &B11, -1);
    opa[0] = wr[0] ? effop(1, 1) : effop(op0, 1);
    opa[2] = wr[2] ? effop(1, 1) : effop(op0, 1);
    wr[0] = 1; wr[2] = 1;
    mmx(h, &Ma, &Mb, Q[0], ld0, opa[0], Q[2], ld0, opa[2]);
  }
  { /* M5 */
    Opd Ma = oadd(&A11, &A12, 1), Mb = B22;
    opa[0] = wr[0] ? effop(1, 2) : effop(op0, 2);
    opa[1] = wr[1] ? effop(1, 1) : effop(op0, 1);
    wr[0] = 1; wr[1] = 1;
    mmx(h, &Ma, &Mb, Q[0], ld0, opa[0], Q[1], ld0, opa[1]);
  }
  { /* M6 */
    Opd Ma = oadd(&A21, &A11, -1), Mb = oadd(&B11, &B12, 1);
    opa[3] = wr[3] ? effop(1, 1) : effop(op0, 1);
    wr[3] = 1;
    mmx(h, &Ma, &Mb, Q[3], ld0, opa[3], 0, 0, 0);
  }
  { /* M7 */
    Opd Ma = oadd(&A12, &A22, -1), Mb = oadd(&B21, &B22, 1);
    opa[0] = wr[0] ? effop(1, 1) : effop(op0, 1);
    wr[0] = 1;
    mmx(h, &Ma, &Mb, Q[0], ld0, opa[0], 0, 0, 0);
  }
}

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);
  else mmx(n, &Ao, &Bo, C, n, 0, 0, 0, 0);
}
#pragma GCC pop_options

评测结果