c - k x k bool 矩阵的快速乘法，其中 8 <= k <= 16

Question

我想找到一种尽可能快的方法来将两个小 bool 矩阵相乘，其中小意味着 8x8、9x9 ... 16x16。这个例程会被大量使用，所以需要非常高效，所以请不要建议直截了当的解决方案应该足够快。

对于 8x8 和 16x16 的特殊情况，基于 solution found here，我已经有了相当有效的实现。 ，我们将整个矩阵分别视为 uint64_t 或 uint64_t[4]。在我的机器上，这比直接实现快大约 70-80 倍。

但是，在 8 < k < 16 的情况下，我真的不知道如何利用任何合理的表示来启用上述这些巧妙的技巧。

所以基本上，我愿意接受任何使用任何类型的(矩阵的)表示和函数签名的建议。您可以假设这针对 32 位或 64 位架构(选择最适合您的建议)

Answer 1

给定两个 4x4 矩阵 a= 0010,0100,1111,0001，b=1100,0001,0100,0100，可以先计算转置 b' = 1000,1011,0000,0100。

然后得到矩阵 M(i,j)=a x b mod 2 == popcount(a[i]&b[j]) & 1;//或奇偶校验

从中可以注意到，只要位 vector 适合计算机字，复杂度只会在 n^2 中增长。

这至少可以加速 8x8 矩阵，前提是可以使用一些特殊的排列和位选择操作。可以用 vector 中的 NxN 位精确地迭代 N 次。 (所以 16x16 几乎是极限)。

每一步都由累加组成，即 Result(n+1) = Result(n) XOR A(n) .& B(n)，其中 Result(0) = 0，A(n) 是 A <<< n , 和 '<<<' == 元素的列旋转，其中 B(n) 从矩阵 B 复制对角线元素:

    a b c          a e i          d h c          g b f
B=  d e f  B(0) =  a e i  B(1) =  d h c   B(2) = g b f
    g h i          a e i          d h c          g b f

在进一步考虑之后，更好的选择是^^^(按行旋转)矩阵 B 并选择 A(n) == 从 A 复制的列对角线:

    a b c         a a a           b b b           c c c 
A=  d e f  A(0) = e e e , A(1) =  f f f,  A(2) =  d d d 
    g h i         i i i           g g g           h h h 

编辑 为了让后来的读者受益，我建议在可移植 C 中针对 W<=16 位矩阵乘法的完整解决方案。

#include <stdint.h>
void matrix_mul_gf2(uint16_t *a, uint16_t *b, uint16_t *c)
{
    // these arrays can be read in two successive xmm registers or in a single ymm
    uint16_t D[16];      // Temporary
    uint16_t C[16]={0};  // result
    uint16_t B[16];  
    uint16_t A[16];
    int i,j;
    uint16_t top_row;
    // Preprocess B (while reading from input) 
    // -- "un-tilt" the diagonal to bit position 0x8000
    for (i=0;i<W;i++) B[i]=(b[i]<<i) | (b[i]>>(W-i));
    for (i=0;i<W;i++) A[i]=a[i];  // Just read in matrix 'a'
    // Loop W times
    // Can be parallelized 4x with MMX, 8x with XMM and 16x with YMM instructions
    for (j=0;j<W;j++) {
        for (i=0;i<W;i++) D[i]=((int16_t)B[i])>>15;  // copy sign bit to rows
        for (i=0;i<W;i++) B[i]<<=1;                  // Prepare B for next round
        for (i=0;i<W;i++) C[i]^= A[i]&D[i];          // Add the partial product

        top_row=A[0];
        for (i=0;i<W-1;i++) A[i]=A[i+1];
        A[W-1]=top_row;
    }
    for (i=0;i<W;i++) c[i]=C[i];      // return result
}