在C语言中进行二进制向量和矩阵操作

Question

在C语言中进行二进制向量和矩阵操作

cmatrixbinarytime-complexityspace-complexity

6

我正在尝试用C实现一种数据结构，可以有效地操作**二进制**矩阵（仅包含1或0）。我将解释我需要对此矩阵应用的操作，并想知道使用什么最佳数据结构？

这些操作在F_2域中进行（这意味着1 + 1 = 0，其他操作保持不变）。我有一个称为H的k*n矩阵（其中k < n）。最多，k = 2325，n = 3009。

我将使用仅使用行交换和行加法部分对角化它。完成后，我将不再使用行操作，并在此矩阵上进行大量（！）的列加法操作（我的意思是约((n-k)/2)³个列加法）。 我考虑用于矩阵的数据结构： 对于矩阵系数，我考虑一次在单个无符号整数中存储多个位的序列。例如，我可以将序列(11001011)存储到uint8_t 203中（从二进制转换为十进制）

这是一个好主意吗？

如果我这样做，我有两个选择：

我可以使用uint16_t或uint64_t系数将我的矩阵H分成许多4*4或8*8子矩阵。

这是一个好的选择（在时间效率方面），如果是，更好使用uint16_t还是uint64_t？

否则，我考虑将每行存储在多个uint32_t或uint64_t中，然后进行部分对角化。接下来切换到一种结构，将矩阵编码为n列向量以处理剩余操作。

您认为这更有效吗？

无论我使用哪种方法，我都需要有效地访问无符号整数（uint16、32或64）的第n位。我该怎么做？

- Celerio

请问您能否解释一下二项式系数的作用是什么？最终结果是一个F_2矩阵吗？ - pablo1977

我在编辑时加入了一些明确的内容。我所说的二项式系数是时间复杂度的指示，与数据结构无关。在这种情况下，没有“最终结果”，我只是在寻找存储仅包含1和0的矩阵的最佳方法，同时给出后续需要应用的操作。 - Celerio

3个回答

3

当你提到类型uint16_t，uint64_t等时，我猜这是为矩阵系数而设的。因此，你应该比我们更清楚你正在操作哪个值：如果你可能会生成巨大的数字，那么你需要一个大型的类型来避免溢出。关于效率，我怀疑你在速度方面不会感觉到差异，但你可以通过选择较小的类型来节省一些空间。

总之，这都是优化问题：你不应该太纠结于强类型。首先，使用char（或uint8_t）应该就可以了，因为你只处理1和0。

我不认为从char matrix[][]切换到typedef struct matrix_columns有任何意义，我认为你可以通过明智地使用行和列索引来执行你的操作。

最后，要获取unsigned int coef中位置的位：

unsigned int bit = (coef>>i) & 1;

[编辑] 此外，这里有一个关于二进制矩阵运算（乘法、加法和异或）的相关问题。

- n0p

是的，我在谈论矩阵系数。这些系数为1或0，因此它们不是很大的数字。但我的问题不太清楚。我编辑了一下，以使其更加明确。而uint8_t、uint32_t、uint64_t是否会改变处理器应用异或运算的速度呢？ - Celerio

我认为这取决于您的处理器架构：它可以在一个时钟周期内处理多少位。如果速度对您来说是至关重要的问题，您可以执行基准测试：尝试使用不同类型的系数进行多次运行。 - n0p

2

+1 使你的问题更清晰：你想做什么（将几个位存储在 unsigned int 中）实际上是一个 bitfield。在你的情况下，这可能是一个不错的优化，但它很棘手且难以阅读/维护。如果你想管理一个矩阵-> 实现一个矩阵：char matrix[][]。优化可以在完成并分析你的程序后进行 ;) - n0p

我刚刚取消了你的答案，因为我想听取其他人的意见。我对位域进行了一些研究，并尝试了一些代码片段来比较暴露的思路，我仍然觉得我在帖子中提出的一些想法可能比char matrix[][]更好。顺便说一下，我已经编写了需要这种数据结构的程序。 :) - Celerio

2

我认为将每行存储在多个uintX_t中是一个好主意，但我会选择X以匹配处理器的字长。这样，您可以一次对多个矩阵条目求和：

uint8_t a = 3; // 00000011
uint8_t b = 2; // 00000010
b ^= a; //results 00000001 
        //for XOR is just addition in F_2

我认为你的问题瓶颈在于列运算，对吗？所以，在进行部分对角化后，为什么不像对部分对角化那样对其转置进行列加法操作呢？虽然这需要逐位完成转置，但之后的列加法会更加简单。假设你有一个64x64的矩阵：

uint64_t matrix[64]; // This means that matrix[i] is the i-th row

你如何对第j列到第l列求和（模2），其中j和l的取值范围在1到64之间？您需要迭代所有行（假设矩阵行中更高的位置位于较不显著的位数）：

int i;
for (i=0; i<64; i++){
    matrix[i] ^= (matrix[i] & (1<<(64-j))) << (64-l);
}

但是，如果你已经把matrix转置为matrixT，那么上述for循环将等同于：

matrixT[l] ^= matrixT[j];

处理器会在一步完成这个操作（我猜是这样，但不确定）。因此，如果你在之后花时间进行转置，你将从处理器的算术运算能力中受益，因为它可以在一步中对其字长数据执行算术运算。

- Tarc

谢谢您确认我的直觉。我相当确定处理器使用单个步骤对两个uint64_t执行XOR操作（我的处理器是64位的）。我将尝试这样做，并与coconop建议的在char矩阵上执行所需的执行时间进行比较。与其转置，我正在考虑反转系数描述矩阵的方式（按列或按行），但我想这两种方法都是等效的。 - Celerio

也许最初的部分对角化可以像我在答案中描述的那样进行，只需使用C矩阵来存储列而不是行。我想这种准对角化并不值得为了获取矩阵转置而麻烦吧？ - Tarc

我的意思是，如果你反转矩阵系数存储的方式并选择以列为单位进行存储，那么最初的对角化就会受到惩罚。 - Tarc

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Nominal Animal · Accepted Answer

为了获得最佳性能，在行交换和行加法中使用一组行指针。使用<stdint.h>，以及支持的字长最小的快速无符号整数类型 -- 我建议使用uint_fast32_t，除非您打算在16位或8位处理器上运行。

当所有行交换和行加法完成后，转置数组。虽然这个操作比较“慢”，但以下列操作将非常快，可以抵消转置成本。

请考虑以下内容：

#include <stdint.h>
#include <limits.h>

typedef uint_fast32_t  word;
#define WORD_BITS      (CHAR_BIT * sizeof (word))

typedef struct {
    int    rows;  /* Number of row vectors */
    int    cols;  /* Number of defined bits in each row */
    int    words; /* Number of words per row vector */
    word **row;   /* Array of pointers */
} row_matrix;

typedef struct {
    int    rows;  /* Number of defined bits in each column */
    int    cols;  /* Number of col vectors */
    int    words; /* Number of words per col vector */
    word **col;
} col_matrix;

虽然您可以使用单个类型来描述这两种矩阵形式，但使用不同的类型使代码和函数更易于维护。最终您可能会有一些重复的代码，但与具有清晰、直观类型相比，这是一个小问题。

在32位系统上，uint_fast32_t通常是32位类型。在64位系统上，它通常是64位。 WORD_BITS宏扩展为word中的位数 - 它并不总是32！

编号位的最简单方法是将矩阵中最左边的位指定为位0，并将位存储在每个字中的最低有效位。如果您有row_matrix * rm，则行row，列col处的位为：

!!(rm->row[row][col / WORD_BITS] & ((word)1U << (col % WORD_BITS)))

!!是非非运算符：如果参数非零，则返回1，否则返回0。因为我们从单词中屏蔽了一个位，所以“位被设置”的值否则将是2的某个幂次方（1、2、4、8、16、32、64等）。

要设置位，请使用：

rm->row[row][col / WORD_BITS] |= (word)1U << (col % WORD_BITS);

为了清除一个比特位，你需要用一个掩码进行按位与运算，该掩码除目标比特位外都为1。可以使用not操作符~轻松实现这一点：

rm->row[row][col / WORD_BITS] &= ~((word)1U << (col % WORD_BITS));

col_matrix *cm的相应操作如下：

!!(cm->col[col][row / WORD_BITS] & ((word)1U << (row % WORD_BITS)))
cm->col[col][row / WORD_BITS] |= (word)1U << (row % WORD_BITS);
cm->col[col][row / WORD_BITS] &= ~((word)1U << (row % WORD_BITS));

尽管除法/和模数（或余数）%通常比加法、减法甚至乘法慢，但在所有广泛使用的架构上，WORD_BITS将是一个二次方编译时常量。我知道的所有编译器都会将上述内容转换为快速位移和二进制AND操作。

要将行srcrow添加到行dstrow，您只需对所有单词进行二进制异或：

{
    const word *const src = rm->row[srcrow];
    word *const       dst = rm->row[dstrow];
    int               w = rm->words;
    while (w-->0)
        dst[w] ^= src[w];
}

同样地，对于列矩阵，

{
    const word *const src = cm->col[srccol];
    word *const       dst = cm->col[dstcol];
    int               w = cm->words;
    while (w-->0)
        dst[w] ^= src[w];
}

请注意，如果你合并超过两行，你可以非常高效地完成；这比连续执行加法要快得多。英特尔和AMD的CPU非常擅长预测上述模式，所以你可以使用多个源行/列。此外，目标不一定必须参与结果，尽管如果我猜对了你正在实现的算法，我想你希望它参与其中。

如果你知道目标架构具有SSE2或更好的性能，甚至是AVX，你可以分别使用代码emmintrin.h或immintrin.h头文件来为编译器内置类型和运算符提供支持，允许你一次性异或128位和256位；这有时会给你带来相当大的提升。

由于向量类型需要C标准所谓的“超出对齐”，因此您还需要包括mm_malloc.h，并使用_mm_malloc()和_mm_free()来为单词数据分配行/列向量-显然将单词四舍五入以便您可以将行/列作为适当的整数字类型(适用于SSE*的__m128i，适用于AVX的__m256i)访问。

个人而言，我总是先实现未向量化的版本，然后进行一些“恶意”测试用例进行单元测试，然后再考虑向量化它。这样做的好处是，您可以将未向量化的版本作为初步版本提供给那些将使用它的人，并且您可以比较向量化和未向量化用例之间的测试结果，以查看哪一个有错误。

转置操作非常简单，尽管我建议使用三重循环：最内部循环在单个字中循环位。此外，你可能想检查哪个顺序——行或列主要——对于外循环效果最好；取决于矩阵大小，你可能会看到很大的差异。(这是由于缓存行为：你希望CPU能够预测访问模式，而不必重新加载相同的缓存行。在最好的情况下，在至多几年前的AMD和Intel x86-64处理器上，如果两个矩阵都适合缓存，你可以获得接近缓存速度的速度。)

以上所有内容都可以实现在一个头文件中——甚至包括向量化版本，如果目标架构支持SSE2/AVX——因此应该不难实现。

有问题吗？