你会如何转置二进制矩阵？
How would you transpose a binary matrix?
我在 C++ 中有二进制矩阵，我用 8 位值的向量表示。
例如下面的矩阵：
1 0 1 0 1 0 1
0 1 1 0 0 1 1
0 0 0 1 1 1 1
表示为：
const uint8_t matrix[] = {
0b01010101,
0b00110011,
0b00001111,
};
我之所以这样做是因为计算这样一个矩阵和一个 8 位向量的乘积变得非常简单和高效（每行只需一个按位 AND 和一个奇偶校验计算），这比单独计算每个位要好得多。
我现在正在寻找一种转置此类矩阵的有效方法，但我一直无法弄清楚如何在无需手动计算每一位的情况下进行转置。
澄清一下，对于上面的例子，我想从转置中得到以下结果：
const uint8_t transposed[] = {
0b00000000,
0b00000100,
0b00000010,
0b00000110,
0b00000001,
0b00000101,
0b00000011,
0b00000111,
};
注意：我更喜欢可以用任意大小的矩阵计算这个的算法，但我也对只能处理特定大小的算法感兴趣。
我的建议是，你不做转置，而是在矩阵数据中添加一位信息，指示矩阵是否转置。
现在，如果要将转置矩阵与向量相乘，则与将左侧矩阵与向量相乘（然后转置）相同。这很简单：只需对 8 位数字进行一些 xor 运算即可。
然而，这会使一些其他操作变得复杂（例如，添加两个矩阵）。但是在评论中你说乘法正是你想要优化的。
我的建议是使用查找 table 来加快处理速度。
另一件需要注意的事情是矩阵的当前定义最大尺寸为 8x8 位。这适合 uint64_t，因此我们可以利用它来发挥我们的优势，尤其是在使用 64 位平台时。
我已经使用查找 table 找到了一个简单的示例，您可以在下面找到它，运行 使用：http://www.tutorialspoint.com/compile_cpp11_online.php 在线编译器。
示例代码
#include <iostream>
#include <bitset>
#include <stdint.h>
#include <assert.h>
using std::cout;
using std::endl;
using std::bitset;
/* Static lookup table */
static uint64_t lut[256];
/* Helper function to print array */
template<int N>
void print_arr(const uint8_t (&arr)[N]){
for(int i=0; i < N; ++i){
cout << bitset<8>(arr[i]) << endl;
}
}
/* Transpose function */
template<int N>
void transpose_bitmatrix(const uint8_t (&matrix)[N], uint8_t (&transposed)[8]){
assert(N <= 8);
uint64_t value = 0;
for(int i=0; i < N; ++i){
value = (value << 1) + lut[matrix[i]];
}
/* Ensure safe copy to prevent misalignment issues */
/* Can be removed if input array can be treated as uint64_t directly */
for(int i=0; i < 8; ++i){
transposed[i] = (value >> (i * 8)) & 0xFF;
}
}
/* Calculate lookup table */
void calculate_lut(void){
/* For all byte values */
for(uint64_t i = 0; i < 256; ++i){
auto b = std::bitset<8>(i);
auto v = std::bitset<64>(0);
/* For all bits in current byte */
for(int bit=0; bit < 8; ++bit){
if(b.test(bit)){
v.set((7 - bit) * 8);
}
}
lut[i] = v.to_ullong();
}
}
int main()
{
calculate_lut();
const uint8_t matrix[] = {
0b01010101,
0b00110011,
0b00001111,
};
uint8_t transposed[8];
transpose_bitmatrix(matrix, transposed);
print_arr(transposed);
return 0;
}
工作原理
您的 3x8 矩阵将转置为 8x3 矩阵，以 8x8 数组表示。
问题是你想转换位，你的 "horizontal" 表示为垂直的，分为几个字节。
正如我上面提到的，我们可以利用输出 (8x8) 始终适合 uint64_t 这一事实。我们将利用它来发挥我们的优势，因为现在我们可以使用 uint64_t 来编写 8 字节数组，但我们也可以将它用于加法、异或等，因为我们可以在 64 位上执行基本算术运算整数。
3x8 矩阵（输入）中的每个条目都是 8 位宽，为了优化处理，我们首先生成 256 个条目查找 table（对于每个字节值）。该条目本身是一个 uint64_t 并将包含位的旋转版本。
示例：
byte = 0b01001111 = 0x4F
lut[0x4F] = 0x0001000001010101 = (uint8_t[]){ 0, 1, 0, 0, 1, 1, 1, 1 }
现在进行计算：
我们使用 uint64_t 进行计算，但请记住，在水下它将代表一个 uint8_t[8] 数组。我们简单地移动当前值（从 0 开始），查找我们的第一个字节并将其添加到当前值。
这里的'magic'是查找table中uint64_t的每个字节要么是1要么是0所以它只会设置最低有效位（每个字节的).移动 uint64_t 将移动每个字节，只要我们确保 我们不这样做超过 8 次！ 我们可以单独对每个字节进行操作。
问题
正如有人在评论中指出的那样：Translate(Translate(M)) != M 所以如果你需要这个，你需要做一些额外的工作。
可以通过直接映射 uint64_t 而不是 uint8_t[8] 数组来提高性能，因为它省略了 "safe-copy" 以防止对齐问题。
我花了更多时间寻找解决方案，并且找到了一些不错的解决方案。
SSE2 方式
在现代 x86 CPU 上，使用 SSE2 指令可以非常高效地完成二进制矩阵的转置。使用此类指令可以处理 16×8 矩阵。
此解决方案的灵感来自 this blog post by mischasan，并且远远优于我迄今为止对这个问题的所有建议。
思路很简单：
#include <emmintrin.h>- 将16个
uint8_t个变量打包成一个__m128i
- 使用
_mm_movemask_epi8获取每个字节的MSB，产生uint16_t
- 使用
_mm_slli_epi64将128位寄存器移位一位
- 重复直到你得到所有 8
uint16_ts
通用 32 位解决方案
不幸的是，我还需要在 ARM 上进行这项工作。实施 SSE2 版本后，很容易找到 NEON 等效项，但是 Cortex-M CPU，（与 Cortex-M 相反A) 没有 SIMD 功能，所以 NEON 目前对我来说用处不大。
注意：因为Cortex-M没有原生的64位算法，我无法在任何建议将 8x8 块视为 uint64_t 的答案中使用这些想法。大多数具有 Cortex-M CPU 的微控制器也没有太多内存，所以我更喜欢在没有查找的情况下完成所有这些操作 table.
经过一番思考，可以使用普通的 32 位算法和一些巧妙的编码来实现相同的算法。这样，我一次可以处理 4×8 个块。一位同事建议，神奇之处在于 32 位乘法的工作方式：您可以找到一个 32 位数字，您可以将其乘以，然后每个字节的 MSB 在高 32 位中彼此相邻结果。
- 在 32 位变量中打包 4
uint8_ts
- 屏蔽每个字节的第一位（使用
0x80808080）
- 乘以
0x02040810
- 乘法的高32位取低4位
- 一般情况下，可以屏蔽每个字节中的第N位（屏蔽右移N位），然后乘以幻数，左移N位。这里的优点是，如果您的编译器足够聪明以展开循环，则掩码和 'magic number' 都会成为编译时常量，因此移动它们不会导致任何性能损失。最后一系列的 4 位有一些问题，因为一个 LSB 丢失了，所以在那种情况下我需要将输入左移 8 位并使用与第一个 4 位系列相同的方法。
如果你用两个 4×8 的块来做这个，那么你就可以完成一个 8x8 的块并排列结果位，以便一切都进入正确的位置。
我添加了一个新的 awnser 而不是编辑我原来的 awnser 以使其更加可见（不幸的是没有评论权）。
在您自己的 awnser 中，您添加了第一个中没有的附加要求：它必须在 ARM Cortex-M 上工作
我确实在我原来的 awnser 中为 ARM 提出了一个替代解决方案，但省略了它，因为它不是问题的一部分并且似乎偏离了主题（主要是因为 C++ 标签）。
ARM具体解决Cortex-M：
部分或大部分 Cortex-M 3/4 有一个位带区域，可用于您所需要的，它将位扩展为 32 位字段，该区域可用于执行原子位操作。
如果您将数组放在位带区域中，它将在位带区域中有一个 'exploded' 镜像，您可以在其中对位本身使用移动操作。如果你做一个循环，编译器肯定能够展开和优化到只移动操作。
如果你真的想要，你甚至可以设置一个 DMA 控制器来处理整批转置操作，并完全从 cpu 卸载它:)
也许这对您仍然有帮助。
这有点晚了，但我今天偶然发现了这个交流。
如果您查看 Hacker's Delight，第 2 版，从第 141 页开始，有几种有效转置布尔数组的算法。
它们非常有效：我的一位同事获得了大约 10 倍的因子
在 X86 上与简单编码相比加速。
这是 Jay Foad 给我的关于快速布尔矩阵的电子邮件正文
转置：
布尔转置算法的核心是一个我称之为 transpose8x8 的函数，它转置一个 64 位字中的 8x8 布尔矩阵（按从 MSB 到 LSB 的行主要顺序）。要转置任何宽度和高度为 8 的倍数的矩形矩阵，请将其分解为 8x8 块，单独转置每个块并将它们存储在输出中的适当位置。要加载 8x8 块，您必须加载 8 个单独的字节并将它们移位和或运算为 64 位字。同样的存储方式。
transpose8x8 的普通 C 实现依赖于这样一个事实，即与前导对角线平行的任何对角线上的所有位都移动相同的距离 up/down 和 left/right。例如，前导对角线上方的所有位都必须向左移动一位并向下移动一位，即在打包的 64 位字中向右移动 7 位。这导致了这样的算法：
transpose8x8(word) {
return
(word & 0x0100000000000000) >> 49 // top right corner
| (word & 0x0201000000000000) >> 42
| ...
| (word & 0x4020100804020100) >> 7 // just above diagonal
| (word & 0x8040201008040201) // leading diagonal
| (word & 0x0080402010080402) << 7 // just below diagonal
| ...
| (word & 0x0000000000008040) << 42
| (word & 0x0000000000000080) << 49; // bottom left corner
}
这 运行 比以前的实现快大约 10 倍，后者从内存中的源字节中单独复制每个位并将其合并到内存中的目标字节中。
或者，如果您有 PDEP 和 PEXT 指令，您可以实现完美的随机播放，并使用它来执行 Hacker's Delight 中提到的转置。这要快得多（但我手边没有时间）：
shuffle(word) {
return pdep(word >> 32, 0xaaaaaaaaaaaaaaaa) | pdep(word, 0x5555555555555555);
} // outer perfect shuffle
transpose8x8(word) { return shuffle(shuffle(shuffle(word))); }
POWER 的 vgbbd 指令在单个指令中有效地实现了整个 transpose8x8（并且由于它是一个 128 位向量指令，它独立地在低 64 位和高 64 位）。这比普通的 C 实现提速了大约 15%。 （只有 15%，因为虽然位旋转速度快得多，但总体 运行 时间现在主要由加载 8 个字节和 assemble 将它们加载到 transpose8x8 的参数所花费的时间支配, 并将结果存储为 8 个单独的字节。)
这是我在 gitub 上发布的内容 (mischasan/sse2/ssebmx.src)
将 INP() 和 OUT() 更改为使用归纳变量可以分别保存一个 IMUL。
AVX256 的速度是它的两倍。
AVX512 不是一个选项，因为没有 _mm512_movemask_epi8().
#include <stdint.h>
#include <emmintrin.h>
#define INP(x,y) inp[(x)*ncols/8 + (y)/8]
#define OUT(x,y) out[(y)*nrows/8 + (x)/8]
void ssebmx(char const *inp, char *out, int nrows, int ncols)
{
int rr, cc, i, h;
union { __m128i x; uint8_t b[16]; } tmp;
// Do the main body in [16 x 8] blocks:
for (rr = 0; rr <= nrows - 16; rr += 16)
for (cc = 0; cc < ncols; cc += 8) {
for (i = 0; i < 16; ++i)
tmp.b[i] = INP(rr + i, cc);
for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
*(uint16_t*)&OUT(rr, cc + i) = _mm_movemask_epi8(tmp.x);
}
if (rr == nrows) return;
// The remainder is a row of [8 x 16]* [8 x 8]?
// Do the [8 x 16] blocks:
for (cc = 0; cc <= ncols - 16; cc += 16) {
for (i = 8; i--;)
tmp.b[i] = h = *(uint16_t const*)&INP(rr + i, cc),
tmp.b[i + 8] = h >> 8;
for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
OUT(rr, cc + i) = h = _mm_movemask_epi8(tmp.x),
OUT(rr, cc + i + 8) = h >> 8;
}
if (cc == ncols) return;
// Do the remaining [8 x 8] block:
for (i = 8; i--;)
tmp.b[i] = INP(rr + i, cc);
for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
OUT(rr, cc + i) = _mm_movemask_epi8(tmp.x);
}
HTH.
受罗伯茨回答的启发，可以利用Arm Neon中的多项式乘法来分散比特--
inline poly8x16_t mull_lo(poly8x16_t a) {
auto b = vget_low_p8(a);
return vreinterpretq_p8_p16(vmull_p8(b,b));
}
inline poly8x16_t mull_hi(poly8x16_t a) {
auto b = vget_high_p8(a);
return vreinterpretq_p8_p16(vmull_p8(b,b));
}
auto a = mull_lo(word);
auto b = mull_lo(a), c = mull_hi(a);
auto d = mull_lo(b), e = mull_hi(b);
auto f = mull_lo(c), g = mull_hi(c);
然后 vsli 可用于成对组合位。
auto ab = vsli_p8(vget_high_p8(d), vget_low_p8(d), 1);
auto cd = vsli_p8(vget_high_p8(e), vget_low_p8(e), 1);
auto ef = vsli_p8(vget_high_p8(f), vget_low_p8(f), 1);
auto gh = vsli_p8(vget_high_p8(g), vget_low_p8(g), 1);
auto abcd = vsli_p8(ab, cd, 2);
auto efgh = vsli_p8(ef, gh, 2);
return vsli_p8(abcd, efgh, 4);
Clang 优化此代码以避免 vmull2 指令，大量使用 ext q0,q0,8 到 vget_high_p8。
迭代方法可能不仅速度更快，而且使用的寄存器更少，并且还可以简化 2 倍或更多的吞吐量。
// transpose bits in 2x2 blocks, first 4 rows
// x = a b|c d|e f|g h a i|c k|e m|g o | byte 0
// i j|k l|m n|o p b j|d l|f n|h p | byte 1
// q r|s t|u v|w x q A|s C|u E|w G | byte 2
// A B|C D|E F|G H r B|t D|v F|h H | byte 3 ...
// ----------------------
auto a = (x & 0x00aa00aa00aa00aaull);
auto b = (x & 0x5500550055005500ull);
auto c = (x & 0xaa55aa55aa55aa55ull) | (a << 7) | (b >> 7);
// transpose 2x2 blocks (first 4 rows shown)
// aa bb cc dd aa ii cc kk
// ee ff gg hh -> ee mm gg oo
// ii jj kk ll bb jj dd ll
// mm nn oo pp ff nn hh pp
auto d = (c & 0x0000cccc0000ccccull);
auto e = (c & 0x3333000033330000ull);
auto f = (c & 0xcccc3333cccc3333ull) | (d << 14) | (e >> 14);
// Final transpose of 4x4 bit blocks
auto g = (f & 0x00000000f0f0f0f0ull);
auto h = (f & 0x0f0f0f0f00000000ull);
x = (f & 0xf0f0f0f00f0f0f0full) | (g << 28) | (h >> 28);
在 ARM 中，每个步骤现在可以由 3 条指令组成：
auto tmp = vrev16_u8(x);
tmp = vshl_u8(tmp, plus_minus_1); // 0xff01ff01ff01ff01ull
x = vbsl_u8(mask_1, x, tmp); // 0xaa55aa55aa55aa55ull
tmp = vrev32_u16(x);
tmp = vshl_u16(tmp, plus_minus_2); // 0xfefe0202fefe0202ull
x = vbsl_u8(mask_2, x, tmp); // 0xcccc3333cccc3333ull
tmp = vrev64_u32(x);
tmp = vshl_u32(tmp, plus_minus_4); // 0xfcfcfcfc04040404ull
x = vbsl_u8(mask_4, x, tmp); // 0xf0f0f0f00f0f0f0full
我在 C++ 中有二进制矩阵，我用 8 位值的向量表示。
例如下面的矩阵：
1 0 1 0 1 0 1
0 1 1 0 0 1 1
0 0 0 1 1 1 1
表示为：
const uint8_t matrix[] = {
0b01010101,
0b00110011,
0b00001111,
};
我之所以这样做是因为计算这样一个矩阵和一个 8 位向量的乘积变得非常简单和高效（每行只需一个按位 AND 和一个奇偶校验计算），这比单独计算每个位要好得多。
我现在正在寻找一种转置此类矩阵的有效方法，但我一直无法弄清楚如何在无需手动计算每一位的情况下进行转置。
澄清一下，对于上面的例子，我想从转置中得到以下结果：
const uint8_t transposed[] = {
0b00000000,
0b00000100,
0b00000010,
0b00000110,
0b00000001,
0b00000101,
0b00000011,
0b00000111,
};
注意：我更喜欢可以用任意大小的矩阵计算这个的算法，但我也对只能处理特定大小的算法感兴趣。
我的建议是，你不做转置，而是在矩阵数据中添加一位信息，指示矩阵是否转置。
现在，如果要将转置矩阵与向量相乘，则与将左侧矩阵与向量相乘（然后转置）相同。这很简单：只需对 8 位数字进行一些 xor 运算即可。
然而，这会使一些其他操作变得复杂（例如，添加两个矩阵）。但是在评论中你说乘法正是你想要优化的。
我的建议是使用查找 table 来加快处理速度。
另一件需要注意的事情是矩阵的当前定义最大尺寸为 8x8 位。这适合 uint64_t，因此我们可以利用它来发挥我们的优势，尤其是在使用 64 位平台时。
我已经使用查找 table 找到了一个简单的示例，您可以在下面找到它，运行 使用：http://www.tutorialspoint.com/compile_cpp11_online.php 在线编译器。
示例代码
#include <iostream>
#include <bitset>
#include <stdint.h>
#include <assert.h>
using std::cout;
using std::endl;
using std::bitset;
/* Static lookup table */
static uint64_t lut[256];
/* Helper function to print array */
template<int N>
void print_arr(const uint8_t (&arr)[N]){
for(int i=0; i < N; ++i){
cout << bitset<8>(arr[i]) << endl;
}
}
/* Transpose function */
template<int N>
void transpose_bitmatrix(const uint8_t (&matrix)[N], uint8_t (&transposed)[8]){
assert(N <= 8);
uint64_t value = 0;
for(int i=0; i < N; ++i){
value = (value << 1) + lut[matrix[i]];
}
/* Ensure safe copy to prevent misalignment issues */
/* Can be removed if input array can be treated as uint64_t directly */
for(int i=0; i < 8; ++i){
transposed[i] = (value >> (i * 8)) & 0xFF;
}
}
/* Calculate lookup table */
void calculate_lut(void){
/* For all byte values */
for(uint64_t i = 0; i < 256; ++i){
auto b = std::bitset<8>(i);
auto v = std::bitset<64>(0);
/* For all bits in current byte */
for(int bit=0; bit < 8; ++bit){
if(b.test(bit)){
v.set((7 - bit) * 8);
}
}
lut[i] = v.to_ullong();
}
}
int main()
{
calculate_lut();
const uint8_t matrix[] = {
0b01010101,
0b00110011,
0b00001111,
};
uint8_t transposed[8];
transpose_bitmatrix(matrix, transposed);
print_arr(transposed);
return 0;
}
工作原理
您的 3x8 矩阵将转置为 8x3 矩阵，以 8x8 数组表示。 问题是你想转换位，你的 "horizontal" 表示为垂直的，分为几个字节。
正如我上面提到的，我们可以利用输出 (8x8) 始终适合 uint64_t 这一事实。我们将利用它来发挥我们的优势，因为现在我们可以使用 uint64_t 来编写 8 字节数组，但我们也可以将它用于加法、异或等，因为我们可以在 64 位上执行基本算术运算整数。
3x8 矩阵（输入）中的每个条目都是 8 位宽，为了优化处理，我们首先生成 256 个条目查找 table（对于每个字节值）。该条目本身是一个 uint64_t 并将包含位的旋转版本。
示例：
byte = 0b01001111 = 0x4F
lut[0x4F] = 0x0001000001010101 = (uint8_t[]){ 0, 1, 0, 0, 1, 1, 1, 1 }
现在进行计算：
我们使用 uint64_t 进行计算，但请记住，在水下它将代表一个 uint8_t[8] 数组。我们简单地移动当前值（从 0 开始），查找我们的第一个字节并将其添加到当前值。
这里的'magic'是查找table中uint64_t的每个字节要么是1要么是0所以它只会设置最低有效位（每个字节的).移动 uint64_t 将移动每个字节，只要我们确保 我们不这样做超过 8 次！ 我们可以单独对每个字节进行操作。
问题
正如有人在评论中指出的那样：Translate(Translate(M)) != M 所以如果你需要这个，你需要做一些额外的工作。
可以通过直接映射 uint64_t 而不是 uint8_t[8] 数组来提高性能，因为它省略了 "safe-copy" 以防止对齐问题。
我花了更多时间寻找解决方案，并且找到了一些不错的解决方案。
SSE2 方式
在现代 x86 CPU 上，使用 SSE2 指令可以非常高效地完成二进制矩阵的转置。使用此类指令可以处理 16×8 矩阵。
此解决方案的灵感来自 this blog post by mischasan，并且远远优于我迄今为止对这个问题的所有建议。
思路很简单：
#include <emmintrin.h>- 将16个
uint8_t个变量打包成一个__m128i - 使用
_mm_movemask_epi8获取每个字节的MSB，产生uint16_t - 使用
_mm_slli_epi64将128位寄存器移位一位 - 重复直到你得到所有 8
uint16_ts
通用 32 位解决方案
不幸的是，我还需要在 ARM 上进行这项工作。实施 SSE2 版本后，很容易找到 NEON 等效项，但是 Cortex-M CPU，（与 Cortex-M 相反A) 没有 SIMD 功能，所以 NEON 目前对我来说用处不大。
注意：因为Cortex-M没有原生的64位算法，我无法在任何建议将 8x8 块视为 uint64_t 的答案中使用这些想法。大多数具有 Cortex-M CPU 的微控制器也没有太多内存，所以我更喜欢在没有查找的情况下完成所有这些操作 table.
经过一番思考，可以使用普通的 32 位算法和一些巧妙的编码来实现相同的算法。这样，我一次可以处理 4×8 个块。一位同事建议，神奇之处在于 32 位乘法的工作方式：您可以找到一个 32 位数字，您可以将其乘以，然后每个字节的 MSB 在高 32 位中彼此相邻结果。
- 在 32 位变量中打包 4
uint8_ts - 屏蔽每个字节的第一位（使用
0x80808080） - 乘以
0x02040810 - 乘法的高32位取低4位
- 一般情况下，可以屏蔽每个字节中的第N位（屏蔽右移N位），然后乘以幻数，左移N位。这里的优点是，如果您的编译器足够聪明以展开循环，则掩码和 'magic number' 都会成为编译时常量，因此移动它们不会导致任何性能损失。最后一系列的 4 位有一些问题，因为一个 LSB 丢失了，所以在那种情况下我需要将输入左移 8 位并使用与第一个 4 位系列相同的方法。
如果你用两个 4×8 的块来做这个，那么你就可以完成一个 8x8 的块并排列结果位，以便一切都进入正确的位置。
我添加了一个新的 awnser 而不是编辑我原来的 awnser 以使其更加可见（不幸的是没有评论权）。
在您自己的 awnser 中，您添加了第一个中没有的附加要求：它必须在 ARM Cortex-M 上工作
我确实在我原来的 awnser 中为 ARM 提出了一个替代解决方案，但省略了它，因为它不是问题的一部分并且似乎偏离了主题（主要是因为 C++ 标签）。
ARM具体解决Cortex-M：
部分或大部分 Cortex-M 3/4 有一个位带区域，可用于您所需要的，它将位扩展为 32 位字段，该区域可用于执行原子位操作。
如果您将数组放在位带区域中，它将在位带区域中有一个 'exploded' 镜像，您可以在其中对位本身使用移动操作。如果你做一个循环，编译器肯定能够展开和优化到只移动操作。
如果你真的想要，你甚至可以设置一个 DMA 控制器来处理整批转置操作，并完全从 cpu 卸载它:)
也许这对您仍然有帮助。
这有点晚了，但我今天偶然发现了这个交流。 如果您查看 Hacker's Delight，第 2 版，从第 141 页开始，有几种有效转置布尔数组的算法。
它们非常有效：我的一位同事获得了大约 10 倍的因子 在 X86 上与简单编码相比加速。
这是 Jay Foad 给我的关于快速布尔矩阵的电子邮件正文 转置：
布尔转置算法的核心是一个我称之为 transpose8x8 的函数，它转置一个 64 位字中的 8x8 布尔矩阵（按从 MSB 到 LSB 的行主要顺序）。要转置任何宽度和高度为 8 的倍数的矩形矩阵，请将其分解为 8x8 块，单独转置每个块并将它们存储在输出中的适当位置。要加载 8x8 块，您必须加载 8 个单独的字节并将它们移位和或运算为 64 位字。同样的存储方式。
transpose8x8 的普通 C 实现依赖于这样一个事实，即与前导对角线平行的任何对角线上的所有位都移动相同的距离 up/down 和 left/right。例如，前导对角线上方的所有位都必须向左移动一位并向下移动一位，即在打包的 64 位字中向右移动 7 位。这导致了这样的算法：
transpose8x8(word) {
return
(word & 0x0100000000000000) >> 49 // top right corner
| (word & 0x0201000000000000) >> 42
| ...
| (word & 0x4020100804020100) >> 7 // just above diagonal
| (word & 0x8040201008040201) // leading diagonal
| (word & 0x0080402010080402) << 7 // just below diagonal
| ...
| (word & 0x0000000000008040) << 42
| (word & 0x0000000000000080) << 49; // bottom left corner
}
这 运行 比以前的实现快大约 10 倍，后者从内存中的源字节中单独复制每个位并将其合并到内存中的目标字节中。
或者，如果您有 PDEP 和 PEXT 指令，您可以实现完美的随机播放，并使用它来执行 Hacker's Delight 中提到的转置。这要快得多（但我手边没有时间）：
shuffle(word) {
return pdep(word >> 32, 0xaaaaaaaaaaaaaaaa) | pdep(word, 0x5555555555555555);
} // outer perfect shuffle
transpose8x8(word) { return shuffle(shuffle(shuffle(word))); }
POWER 的 vgbbd 指令在单个指令中有效地实现了整个 transpose8x8（并且由于它是一个 128 位向量指令，它独立地在低 64 位和高 64 位）。这比普通的 C 实现提速了大约 15%。 （只有 15%，因为虽然位旋转速度快得多，但总体 运行 时间现在主要由加载 8 个字节和 assemble 将它们加载到 transpose8x8 的参数所花费的时间支配, 并将结果存储为 8 个单独的字节。)
这是我在 gitub 上发布的内容 (mischasan/sse2/ssebmx.src) 将 INP() 和 OUT() 更改为使用归纳变量可以分别保存一个 IMUL。 AVX256 的速度是它的两倍。 AVX512 不是一个选项，因为没有 _mm512_movemask_epi8().
#include <stdint.h>
#include <emmintrin.h>
#define INP(x,y) inp[(x)*ncols/8 + (y)/8]
#define OUT(x,y) out[(y)*nrows/8 + (x)/8]
void ssebmx(char const *inp, char *out, int nrows, int ncols)
{
int rr, cc, i, h;
union { __m128i x; uint8_t b[16]; } tmp;
// Do the main body in [16 x 8] blocks:
for (rr = 0; rr <= nrows - 16; rr += 16)
for (cc = 0; cc < ncols; cc += 8) {
for (i = 0; i < 16; ++i)
tmp.b[i] = INP(rr + i, cc);
for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
*(uint16_t*)&OUT(rr, cc + i) = _mm_movemask_epi8(tmp.x);
}
if (rr == nrows) return;
// The remainder is a row of [8 x 16]* [8 x 8]?
// Do the [8 x 16] blocks:
for (cc = 0; cc <= ncols - 16; cc += 16) {
for (i = 8; i--;)
tmp.b[i] = h = *(uint16_t const*)&INP(rr + i, cc),
tmp.b[i + 8] = h >> 8;
for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
OUT(rr, cc + i) = h = _mm_movemask_epi8(tmp.x),
OUT(rr, cc + i + 8) = h >> 8;
}
if (cc == ncols) return;
// Do the remaining [8 x 8] block:
for (i = 8; i--;)
tmp.b[i] = INP(rr + i, cc);
for (i = 8; i--; tmp.x = _mm_slli_epi64(tmp.x, 1))
OUT(rr, cc + i) = _mm_movemask_epi8(tmp.x);
}
HTH.
受罗伯茨回答的启发，可以利用Arm Neon中的多项式乘法来分散比特--
inline poly8x16_t mull_lo(poly8x16_t a) {
auto b = vget_low_p8(a);
return vreinterpretq_p8_p16(vmull_p8(b,b));
}
inline poly8x16_t mull_hi(poly8x16_t a) {
auto b = vget_high_p8(a);
return vreinterpretq_p8_p16(vmull_p8(b,b));
}
auto a = mull_lo(word);
auto b = mull_lo(a), c = mull_hi(a);
auto d = mull_lo(b), e = mull_hi(b);
auto f = mull_lo(c), g = mull_hi(c);
然后 vsli 可用于成对组合位。
auto ab = vsli_p8(vget_high_p8(d), vget_low_p8(d), 1);
auto cd = vsli_p8(vget_high_p8(e), vget_low_p8(e), 1);
auto ef = vsli_p8(vget_high_p8(f), vget_low_p8(f), 1);
auto gh = vsli_p8(vget_high_p8(g), vget_low_p8(g), 1);
auto abcd = vsli_p8(ab, cd, 2);
auto efgh = vsli_p8(ef, gh, 2);
return vsli_p8(abcd, efgh, 4);
Clang 优化此代码以避免 vmull2 指令，大量使用 ext q0,q0,8 到 vget_high_p8。
迭代方法可能不仅速度更快，而且使用的寄存器更少，并且还可以简化 2 倍或更多的吞吐量。
// transpose bits in 2x2 blocks, first 4 rows
// x = a b|c d|e f|g h a i|c k|e m|g o | byte 0
// i j|k l|m n|o p b j|d l|f n|h p | byte 1
// q r|s t|u v|w x q A|s C|u E|w G | byte 2
// A B|C D|E F|G H r B|t D|v F|h H | byte 3 ...
// ----------------------
auto a = (x & 0x00aa00aa00aa00aaull);
auto b = (x & 0x5500550055005500ull);
auto c = (x & 0xaa55aa55aa55aa55ull) | (a << 7) | (b >> 7);
// transpose 2x2 blocks (first 4 rows shown)
// aa bb cc dd aa ii cc kk
// ee ff gg hh -> ee mm gg oo
// ii jj kk ll bb jj dd ll
// mm nn oo pp ff nn hh pp
auto d = (c & 0x0000cccc0000ccccull);
auto e = (c & 0x3333000033330000ull);
auto f = (c & 0xcccc3333cccc3333ull) | (d << 14) | (e >> 14);
// Final transpose of 4x4 bit blocks
auto g = (f & 0x00000000f0f0f0f0ull);
auto h = (f & 0x0f0f0f0f00000000ull);
x = (f & 0xf0f0f0f00f0f0f0full) | (g << 28) | (h >> 28);
在 ARM 中，每个步骤现在可以由 3 条指令组成：
auto tmp = vrev16_u8(x);
tmp = vshl_u8(tmp, plus_minus_1); // 0xff01ff01ff01ff01ull
x = vbsl_u8(mask_1, x, tmp); // 0xaa55aa55aa55aa55ull
tmp = vrev32_u16(x);
tmp = vshl_u16(tmp, plus_minus_2); // 0xfefe0202fefe0202ull
x = vbsl_u8(mask_2, x, tmp); // 0xcccc3333cccc3333ull
tmp = vrev64_u32(x);
tmp = vshl_u32(tmp, plus_minus_4); // 0xfcfcfcfc04040404ull
x = vbsl_u8(mask_4, x, tmp); // 0xf0f0f0f00f0f0f0full