代码

/* Additional credit: Karstadt and Schwartz, alternative-basis multiplication;
Juhani Karppa, On Bilinear Techniques for Matrix Multiplication (Aalto thesis),
Section 3.8, equations (3.19)-(3.21), for the exact 12-addition decomposition.
https://aaltodoc.aalto.fi/bitstreams/c72cf2da-7b39-43a8-b55b-6146765aa54c/download
Implementation and two-scratch product schedule below by pdoom.
*/
/*
Credit: pdoom #86352, https://duck.ac/submission/86352 : original exact
mod-65536 Strassen-Winograd engine, recursive packing and AVX2 pair-dot kernel.
saffah_cc_v41_agg1 #112006, https://duck.ac/submission/112006 : baseline code,
peeled first k pair, cache offsets, and packing that eliminates output shuffles.
saffah_codex_6s_agg2 #110387, https://duck.ac/submission/110387 : fused root
reconstruction/unpacking and constant-size recursion; #102935 : instruction order.
saffah_cc_v41_260924 #96771, https://duck.ac/submission/96771 : cache phase idea.
New experiments by pdoom: exact ring arithmetic only; see variant below.
*/
#pragma GCC optimize("O3,unroll-loops")
#pragma GCC target("avx2")
#include <immintrin.h>
#include <stdint.h>
#include <string.h>

#ifndef BASE
#define BASE 32
#endif
using U = uint16_t;
static constexpr int MAXN = 4096;
alignas(4096) static U pa_[MAXN * MAXN + 8192], pb_[MAXN * MAXN + 8192];
alignas(4096) static U pc_[MAXN * MAXN + 8192], ws_[MAXN * MAXN + 8192];

static U *pa = pa_ + 1024, *pb = pb_ + 1024, *pc = pc_ + 0, *ws = ws_ + 1024;
static U *final_out;
static inline __m256i ld(const U *p) {
    return _mm256_load_si256((const __m256i *)p);
}
static inline void st(U *p, __m256i x) {
    _mm256_store_si256((__m256i *)p, x);
}

// The recursive quadrant layout makes every Strassen addition contiguous.
static void pack(U *d, const U *s, int n, int stride, bool right) {
    if (n > BASE) {
        int h = n / 2, q = h * h;
        pack(d, s, h, stride, right);
        pack(d + q, s + h, h, stride, right);
        pack(d + 2*q, s + h*stride, h, stride, right);
        pack(d + 3*q, s + h*stride + h, h, stride, right);
        for(int i=0;i<q;i+=16) {
            __m256i b=ld(d+q+i), c=ld(d+2*q+i), e=ld(d+3*q+i);
            st(d+3*q+i,_mm256_add_epi16(e,right?_mm256_sub_epi16(c,b):_mm256_add_epi16(c,b)));
        }

    } else if (!right) {
        for(int i=0;i<n;i+=4) for(int k=0;k<n;k+=16) {
            __m256i a=ld(s+i*stride+k), b=ld(s+(i+1)*stride+k);
            __m256i c=ld(s+(i+2)*stride+k), e=ld(s+(i+3)*stride+k);
            __m256i t0=_mm256_unpacklo_epi32(a,b),t1=_mm256_unpacklo_epi32(c,e);
            __m256i t2=_mm256_unpackhi_epi32(a,b),t3=_mm256_unpackhi_epi32(c,e);
            __m256i u0=_mm256_unpacklo_epi64(t0,t1),u1=_mm256_unpackhi_epi64(t0,t1);
            __m256i u2=_mm256_unpacklo_epi64(t2,t3),u3=_mm256_unpackhi_epi64(t2,t3);
            st(d,_mm256_permute2x128_si256(u0,u1,0x20));
            st(d+16,_mm256_permute2x128_si256(u2,u3,0x20));
            st(d+32,_mm256_permute2x128_si256(u0,u1,0x31));
            st(d+48,_mm256_permute2x128_si256(u2,u3,0x31));d+=64;
        }
    } else {
        for(int j=0;j<n;) {
            int w=n-j>=24?24:n-j;
            for(int k=0;k<n;k+=2) {
                for(int t=0;t<w;t+=16) {
                    if(t+16<=w) {
                        __m256i x=_mm256_loadu_si256((const __m256i*)(s+k*stride+j+t));
                        __m256i y=_mm256_loadu_si256((const __m256i*)(s+(k+1)*stride+j+t));
                        st(d,_mm256_unpacklo_epi16(x,y));
                        st(d+16,_mm256_unpackhi_epi16(x,y));d+=32;
                    } else {
                        __m128i x0=_mm_load_si128((const __m128i*)(s+k*stride+j+t));
                        __m128i y0=_mm_load_si128((const __m128i*)(s+(k+1)*stride+j+t));
                        _mm_store_si128((__m128i*)d,     _mm_unpacklo_epi16(x0,y0));
                        _mm_store_si128((__m128i*)(d+8), _mm_unpackhi_epi16(x0,y0));
                        d+=16;
                    }
                }
            }
            j+=w;
        }
    }
}

static void unpack(U *d, const U *s, int n, int stride) {
    if (n > BASE) {
        int h = n/2, q = h*h;
        unpack(d,s,h,stride);
        unpack(d+h,s+q,h,stride);
        unpack(d+h*stride,s+2*q,h,stride);
        unpack(d+h*stride+h,s+3*q,h,stride);
    } else {
        for (int i=0;i<n;++i) {
            U *dp = d + (size_t)i*stride; const U *sp = s + (size_t)i*n;
            int k=0;
            for (;k+32<=n;k+=32) {
                _mm256_stream_si256((__m256i*)(dp+k),    _mm256_load_si256((const __m256i*)(sp+k)));
                _mm256_stream_si256((__m256i*)(dp+k+16), _mm256_load_si256((const __m256i*)(sp+k+16)));
            }
            for (;k<n;k+=16) _mm256_storeu_si256((__m256i*)(dp+k), _mm256_loadu_si256((const __m256i*)(sp+k)));
        }
    }
}

static inline void micro24(const U *a,const U *b,U *c,int n) {
    const U *end=a+4*n; U *c3=c+3*n; intptr_t stride=2*n;
    asm volatile(
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm2\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm4\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm5\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm8\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm10\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm11\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vmovdqa 128(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 288(%[b]), %%ymm12\n\t"
        "vmovdqa 320(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "add $64, %[a]\n\t"
        "add $384, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jae 2f\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 64(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vmovdqa 128(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 160(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 256(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "vmovdqa 288(%[b]), %%ymm12\n\t"
        "vmovdqa 320(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpaddd %%ymm14, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm7, %%ymm7\n\t"
        "vpaddd %%ymm14, %%ymm8, %%ymm8\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm9, %%ymm9\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm15\n\t"
        "vpmaddwd 352(%[b]), %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm10, %%ymm10\n\t"
        "vpaddd %%ymm14, %%ymm11, %%ymm11\n\t"
        "add $64, %[a]\n\t"
        "add $384, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jb 1b\n\t"
        "2:\n\t"
        "vpcmpeqd %%ymm13, %%ymm13, %%ymm13\n\t"
        "vpsrld $16, %%ymm13, %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpand %%ymm13, %%ymm2, %%ymm2\n\t"
        "vpackusdw %%ymm1, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%ymm0, 0(%[c])\n\t"
        "vpackusdw %%ymm2, %%ymm2, %%ymm2\n\t"
        "vpermq $216, %%ymm2, %%ymm2\n\t"
        "vmovdqu %%xmm2, 32(%[c])\n\t"
        "vpand %%ymm13, %%ymm3, %%ymm3\n\t"
        "vpand %%ymm13, %%ymm4, %%ymm4\n\t"
        "vpand %%ymm13, %%ymm5, %%ymm5\n\t"
        "vpackusdw %%ymm4, %%ymm3, %%ymm3\n\t"
        "vmovdqu %%ymm3, 0(%[c],%[s],1)\n\t"
        "vpackusdw %%ymm5, %%ymm5, %%ymm5\n\t"
        "vpermq $216, %%ymm5, %%ymm5\n\t"
        "vmovdqu %%xmm5, 32(%[c],%[s],1)\n\t"
        "vpand %%ymm13, %%ymm6, %%ymm6\n\t"
        "vpand %%ymm13, %%ymm7, %%ymm7\n\t"
        "vpand %%ymm13, %%ymm8, %%ymm8\n\t"
        "vpackusdw %%ymm7, %%ymm6, %%ymm6\n\t"
        "vmovdqu %%ymm6, 0(%[c],%[s],2)\n\t"
        "vpackusdw %%ymm8, %%ymm8, %%ymm8\n\t"
        "vpermq $216, %%ymm8, %%ymm8\n\t"
        "vmovdqu %%xmm8, 32(%[c],%[s],2)\n\t"
        "vpand %%ymm13, %%ymm9, %%ymm9\n\t"
        "vpand %%ymm13, %%ymm10, %%ymm10\n\t"
        "vpand %%ymm13, %%ymm11, %%ymm11\n\t"
        "vpackusdw %%ymm10, %%ymm9, %%ymm9\n\t"
        "vmovdqu %%ymm9, 0(%[c3])\n\t"
        "vpackusdw %%ymm11, %%ymm11, %%ymm11\n\t"
        "vpermq $216, %%ymm11, %%ymm11\n\t"
        "vmovdqu %%xmm11, 32(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end)
        : "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15");
}

static inline void micro16(const U *a,const U *b,U *c,int n) {
    const U *end=a+4*n; U *c3=c+3*n; intptr_t stride=2*n;
    asm volatile(
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm1\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm3\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm4\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm5\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm6\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm7\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vmovdqa 96(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 128(%[b]), %%ymm12\n\t"
        "vmovdqa 160(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "add $64, %[a]\n\t"
        "add $256, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jae 2f\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vmovdqa 32(%[b]), %%ymm13\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vmovdqa 96(%[b]), %%ymm13\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 128(%[b]), %%ymm12\n\t"
        "vmovdqa 160(%[b]), %%ymm13\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "vmovdqa 192(%[b]), %%ymm12\n\t"
        "vmovdqa 224(%[b]), %%ymm13\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpaddd %%ymm14, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpaddd %%ymm14, %%ymm3, %%ymm3\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm4, %%ymm4\n\t"
        "vpaddd %%ymm14, %%ymm5, %%ymm5\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpmaddwd %%ymm13, %%ymm14, %%ymm14\n\t"
        "vpaddd %%ymm15, %%ymm6, %%ymm6\n\t"
        "vpaddd %%ymm14, %%ymm7, %%ymm7\n\t"
        "add $64, %[a]\n\t"
        "add $256, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jb 1b\n\t"
        "2:\n\t"
        "vpcmpeqd %%ymm13, %%ymm13, %%ymm13\n\t"
        "vpsrld $16, %%ymm13, %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpackusdw %%ymm1, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%ymm0, 0(%[c])\n\t"
        "vpand %%ymm13, %%ymm2, %%ymm2\n\t"
        "vpand %%ymm13, %%ymm3, %%ymm3\n\t"
        "vpackusdw %%ymm3, %%ymm2, %%ymm2\n\t"
        "vmovdqu %%ymm2, 0(%[c],%[s],1)\n\t"
        "vpand %%ymm13, %%ymm4, %%ymm4\n\t"
        "vpand %%ymm13, %%ymm5, %%ymm5\n\t"
        "vpackusdw %%ymm5, %%ymm4, %%ymm4\n\t"
        "vmovdqu %%ymm4, 0(%[c],%[s],2)\n\t"
        "vpand %%ymm13, %%ymm6, %%ymm6\n\t"
        "vpand %%ymm13, %%ymm7, %%ymm7\n\t"
        "vpackusdw %%ymm7, %%ymm6, %%ymm6\n\t"
        "vmovdqu %%ymm6, 0(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end)
        : "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15");
}

static inline void micro8(const U *a,const U *b,U *c,int n) {
    const U *end=a+4*n; U *c3=c+3*n; intptr_t stride=2*n;
    asm volatile(
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm0\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm1\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm2\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm3\n\t"
        "vmovdqa 32(%[b]), %%ymm12\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "add $64, %[a]\n\t"
        "add $128, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jae 2f\n\t"
        "1:\n\t"
        "vmovdqa 0(%[b]), %%ymm12\n\t"
        "vpbroadcastd 0(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 4(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 8(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 12(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 32(%[b]), %%ymm12\n\t"
        "vpbroadcastd 16(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 20(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 24(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 28(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 64(%[b]), %%ymm12\n\t"
        "vpbroadcastd 32(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 36(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 40(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 44(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "vmovdqa 96(%[b]), %%ymm12\n\t"
        "vpbroadcastd 48(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm0, %%ymm0\n\t"
        "vpbroadcastd 52(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm1, %%ymm1\n\t"
        "vpbroadcastd 56(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm2, %%ymm2\n\t"
        "vpbroadcastd 60(%[a]), %%ymm14\n\t"
        "vpmaddwd %%ymm12, %%ymm14, %%ymm15\n\t"
        "vpaddd %%ymm15, %%ymm3, %%ymm3\n\t"
        "add $64, %[a]\n\t"
        "add $128, %[b]\n\t"
        "cmp %[end], %[a]\n\t"
        "jb 1b\n\t"
        "2:\n\t"
        "vpcmpeqd %%ymm13, %%ymm13, %%ymm13\n\t"
        "vpsrld $16, %%ymm13, %%ymm13\n\t"
        "vpand %%ymm13, %%ymm0, %%ymm0\n\t"
        "vpackusdw %%ymm0, %%ymm0, %%ymm0\n\t"
        "vpermq $216, %%ymm0, %%ymm0\n\t"
        "vmovdqu %%xmm0, 0(%[c])\n\t"
        "vpand %%ymm13, %%ymm1, %%ymm1\n\t"
        "vpackusdw %%ymm1, %%ymm1, %%ymm1\n\t"
        "vpermq $216, %%ymm1, %%ymm1\n\t"
        "vmovdqu %%xmm1, 0(%[c],%[s],1)\n\t"
        "vpand %%ymm13, %%ymm2, %%ymm2\n\t"
        "vpackusdw %%ymm2, %%ymm2, %%ymm2\n\t"
        "vpermq $216, %%ymm2, %%ymm2\n\t"
        "vmovdqu %%xmm2, 0(%[c],%[s],2)\n\t"
        "vpand %%ymm13, %%ymm3, %%ymm3\n\t"
        "vpackusdw %%ymm3, %%ymm3, %%ymm3\n\t"
        "vpermq $216, %%ymm3, %%ymm3\n\t"
        "vmovdqu %%xmm3, 0(%[c3])\n\t"
        : [a] "+&r"(a), [b] "+&r"(b)
        : [c] "r"(c), [c3] "r"(c3), [s] "r"(stride), [end] "r"(end)
        : "cc", "memory", "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15");
}

static void base(const U *a,const U *b,U *c,int n) {
    int j=0;
    for(;j+24<=n;j+=24)
        for(int i=0;i<n;i+=4) micro24(a+i*n,b+j*n,c+i*n+j,n);
    if(n-j==16)
        for(int i=0;i<n;i+=4) micro16(a+i*n,b+j*n,c+i*n+j,n);
    if(n-j==8)
        for(int i=0;i<n;i+=4) micro8(a+i*n,b+j*n,c+i*n+j,n);
}

template<bool Sub>
static inline void combine(U *d,const U *a,const U *b,int len) {
    for(int i=0;i<len;i+=16) {
        __m256i x=ld(a+i),y=ld(b+i);
        st(d+i,Sub?_mm256_sub_epi16(x,y):_mm256_add_epi16(x,y));
    }
}


template<int N>
static __attribute__((noinline)) void mul(const U *a,const U *b,U *c,U *work) {
 if constexpr(N<=BASE) base(a,b,c,N);
 else {
  constexpr int h=N/2,q=h*h;
  const U *a11=a,*a12=a+q,*a21=a+2*q,*a22=a+3*q;
  const U *b11=b,*b12=b+q,*b21=b+2*q,*b22=b+3*q;
  U *c11=c,*c12=c+q,*c21=c+2*q,*c22=c+3*q;
  U *x=work,*y=work+q,*next=work+2*q;
  mul<h>(a12,b21,c11,next); // Q2
  combine<true>(x,a22,a12,q); combine<false>(y,b12,b22,q);
  mul<h>(x,y,c21,next); // Q6
  combine<true>(x,a22,a21,q); combine<true>(y,b21,b22,q);
  mul<h>(x,y,c12,next); // Q7
  mul<h>(a22,b22,c22,next); // Q4
  for(int i=0;i<q;i+=16)
   st(c22+i,_mm256_sub_epi16(_mm256_sub_epi16(_mm256_add_epi16(ld(c22+i),ld(c12+i)),ld(c11+i)),ld(c21+i)));
  mul<h>(a11,b11,x,next); // Q1
  combine<false>(c11,c11,x,q);
  combine<true>(y,b11,b22,q);
  mul<h>(a21,y,x,next); // Q3
  combine<false>(c21,c21,x,q);
  combine<false>(x,a11,a22,q);
  mul<h>(x,b12,y,next); // Q5
  combine<true>(c12,c12,y,q);
 }
}
template<int N>
static void inverse(U *c) {
 if constexpr(N>BASE) {
  constexpr int h=N/2,q=h*h;
  for(int i=0;i<q;i+=16){
   __m256i b=ld(c+q+i),e=ld(c+3*q+i),d=ld(c+2*q+i);
   st(c+q+i,_mm256_sub_epi16(e,b));
   st(c+2*q+i,_mm256_add_epi16(d,e));
   st(c+3*q+i,_mm256_sub_epi16(_mm256_setzero_si256(),e));
  }
  inverse<h>(c);inverse<h>(c+q);inverse<h>(c+2*q);inverse<h>(c+3*q);
 }
}
void matrix_multiply(int n,const short *A,const short *B,short *C) {
 pack(pa,(const U*)A,n,n,false);pack(pb,(const U*)B,n,n,true);
 mul<4096>(pa,pb,pc,ws);inverse<4096>(pc);unpack((U*)C,pc,n,n);_mm_sfence();
}

评测结果