c - 使用 AVX2 对带有十进制数字数组的 BigInt 进行矢量化随机初始化和打印？
问题描述
如何将我的代码传递给 AVX2 代码并获得与以前相同的结果？
是否可以__m256i在 LongNumInit、LongNumPrint 函数中使用，而不是uint8_t *L或一些类似类型的变量？
我对 AVX 的了解非常有限；我进行了很多调查，但是我不太了解如何转换我的代码，欢迎提出任何建议和解释。
我对 AVX2 中的这段代码非常感兴趣。
void LongNumInit(uint8_t *L, size_t N )
{
for(size_t i = 0; i < N; ++i){
L[i] = myRandom()%10;
}
}
void LongNumPrint( uint8_t *L, size_t N, uint8_t *Name )
{
printf("%s:", Name);
for ( size_t i=N; i>0;--i )
{
printf("%d", L[i-1]);
}
printf("\n");
}
int main (int argc, char **argv)
{
int i, sum1, sum2, sum3, N=10000, Rep=50;
seed = 12345;
// obtain parameters at run time
if (argc>1) { N = atoi(argv[1]); }
if (argc>2) { Rep = atoi(argv[2]); }
// Create Long Nums
unsigned char *V1= (unsigned char*) malloc( N);
unsigned char *V2= (unsigned char*) malloc( N);
unsigned char *V3= (unsigned char*) malloc( N);
unsigned char *V4= (unsigned char*) malloc( N);
LongNumInit ( V1, N ); LongNumInit ( V2, N ); LongNumInit ( V3, N );
//Print last 32 digits of Long Numbers
LongNumPrint( V1, 32, "V1" );
LongNumPrint( V2, 32, "V2" );
LongNumPrint( V3, 32, "V3" );
LongNumPrint( V4, 32, "V4" );
free(V1); free(V2); free(V3); free(V4);
return 0;
}
我在初始代码中获得的结果是：
V1:59348245908804493219098067811457
V2:24890422397351614779297691741341
V3:63392771324953818089038280656869
V4:00000000000000000000000000000000
解决方案
一般来说，这对于 BigInteger 来说是一种糟糕的格式，请参阅https://codereview.stackexchange.com/a/237764 ，了解对 BigInteger使用每字节一个十进制数字的设计缺陷的代码审查，以及您可以/应该做什么.
并查看长整数例程可以从 SSE 中受益吗？对于@Mysticial 关于存储数据的方法的注释，这些方法使 BigInteger 数学的 SIMD 变得实用，特别是部分字算术，其中您的临时变量可能没有“标准化”，让您进行惰性进位处理。
但显然你只是在问这个代码，random-init 和 print 函数，而不是如何在这种格式的两个数字之间进行数学运算。
我们可以很好地向量化这两者。我LongNumPrintName()的是你的替代品。
因为LongNumInit我只是展示了一个构建块，它存储两个 32 字节的块并返回一个递增的指针。循环调用它。（每次调用自然会产生 2 个向量，因此对于小 N，您可以制作一个替代版本。）
LongNumInit
生成包含随机数字的 1 GB 文本文件的最快方法是什么？在 4GHz Skylake 上以大约 33 GB/s 的速度生成以空格分隔的随机 ASCII 十进制数字，包括write()对/dev/null. （这高于 DRAM 带宽；128kiB 的缓存阻塞让存储命中 L2 缓存。内核驱动程序/dev/null甚至不读取用户空间缓冲区。）
它可以很容易地适应 AVX2 版本的void LongNumInit(uint8_t *L, size_t N ). 我的答案使用 AVX2 xorshift128+ PRNG（在 a 的 64 位元素中使用 4 个独立的 PRNG 进行矢量化__m256i），例如AVX/SSE 版本的 xorshift128+。这应该与您的rand() % 10.
它通过乘法逆元将其分解为十进制数字，以通过移位和除以 10 并取模vpmulhuw，使用为什么 GCC 在实现整数除法时使用乘以一个奇怪的数字？. （实际上使用 GNU C 原生向量语法让 GCC 确定魔法常数并发出乘法和移位以方便语法，例如v16u dig1 = v % ten;and v /= ten;）
您可以使用_mm256_packus_epi16将两个 16 位数字的向量打包成 8 位元素，而不是将奇数元素转换为 ASCII' '并将偶数元素转换为 ASCII '0'..'9'。（因此更改vec_store_digit_and_space为打包向量对而不是使用常量进行 ORing，见下文）
使用 gcc、clang 或 ICC（或者希望任何其他能够理解 C99 的 GNU C 方言和英特尔内部函数的编译器）编译它。
有关部分，请参阅https://gcc.gnu.org/onlinedocs/gcc/Vector-Extensions.html，有关内容，请参阅https://software.intel.com/sites/landingpage/IntrinsicsGuide/。还有https://stackoverflow.com/tags/sse/info。__attribute__((vector_size(32)))_mm256_*
#include <immintrin.h>
// GNU C native vectors let us get the compiler to do stuff like %10 each element
typedef unsigned short v16u __attribute__((vector_size(32)));
// returns p + size of stores. Caller should use outpos = f(vec, outpos)
// p must be aligned
__m256i* vec_store_digits(__m256i vec, __m256i *restrict p)
{
v16u v = (v16u)vec;
v16u ten = (v16u)_mm256_set1_epi16(10);
v16u divisor = (v16u)_mm256_set1_epi16(6554); // ceil((2^16-1) / 10.0)
v16u div6554 = v / divisor; // Basically the entropy from the upper two decimal digits: 0..65.
// Probably some correlation with the modulo-based values, especially dig3, but we do this instead of
// dig4 for more ILP and fewer instructions total.
v16u dig1 = v % ten;
v /= ten;
v16u dig2 = v % ten;
v /= ten;
v16u dig3 = v % ten;
// dig4 would overlap much of the randomness that div6554 gets
// __m256i or v16u assignment is an aligned store
v16u *vecbuf = (v16u*)p;
// pack 16->8 bits.
vecbuf[0] = _mm256_packus_epi16(div6554, dig1);
vecbuf[1] = _mm256_packus_epi16(dig2, dig3)
return p + 2; // always a constant number of full vectors
}
插入换行符的逻辑random_decimal_fill_buffer可以完全删除，因为您只需要一个十进制数字的平面数组。只需在循环中调用上述函数，直到填满缓冲区。
处理小尺寸（小于完整向量）：
将 malloc 填充到下一个 32 字节的倍数会很方便，因此在不检查是否可能进入未映射页面的情况下进行 32 字节加载总是安全的。
并使用 C11aligned_alloc获得 32 字节对齐的存储。例如，aligned_alloc(32, (size+31) & -32). 即使 N 是奇数，这也让我们可以进行完整的 32 字节存储。从逻辑上讲，只有缓冲区的前 N 个字节保存我们的真实数据，但我们可以随意添加填充以避免任何额外的条件检查 N 小于 32 或不是 32 的倍数。
不幸的是，缺少 ISO C 和 glibcaligned_realloc并且aligned_calloc. MSVC 确实提供了这些：为什么在大多数平台上没有“aligned_realloc”？允许您有时在对齐缓冲区的末尾分配更多空间而不复制它。如果非平凡可复制对象更改地址，则“try_realloc”对于可能需要运行复制构造函数的 C++ 来说是理想的。强制有时不必要的复制的非表达分配器 API 是我的一个小烦恼。
LongNumPrint
采用uint8_t *Namearg 是糟糕的设计。如果调用者想先打印一个"something:"字符串，他们可以这样做。您的函数应该只printf "%d"为int.
由于您以反向打印顺序存储数字，因此您需要将字节反转到 tmp 缓冲区并'0'..'9'通过 ORing 将 0..9 字节值转换为 ASCII 字符值'0'。然后将该缓冲区传递给fwrite.
具体来说，alignas(32) char tmpbuf[8192];用作局部变量。
您可以在固定大小的块（如 1kiB 或 8kiB）中工作，而不是分配一个潜在的巨大缓冲区。您可能仍希望通过 stdio（而不是write()直接管理自己的 I/O 缓冲）。使用 8kiB 缓冲区，高效fwrite可能只是将其write()直接传递给 stdio 缓冲区，而不是 memcpy。您可能想尝试调整它，但保持 tmp 缓冲区小于 L1d 缓存的一半将意味着它在您编写后重新读取时仍处于高速缓存中。
缓存阻塞使循环边界更加复杂，但对于非常大的 N 来说这是值得的。
一次反转 32 个字节：
您可以通过决定您的数字以 MSD 优先顺序存储来避免这项工作，但是如果您确实想要实现加法，则必须从末端向后循环。
您的函数可以用 SIMD 实现_mm_shuffle_epi8以反转 16 字节块，从数字数组的末尾开始并写入 tmp 缓冲区的开头。
或者更好的是，加载vmovdqu/ vinserti12816 字节加载以馈送_mm256_shuffle_epi8到通道内的字节反转，设置 32 字节存储。
在 Intel CPU 上，vinserti128解码为 load+ALU uop，但它可以在任何向量 ALU 端口上运行，而不仅仅是 shuffle 端口。因此，两个 128 位加载比 256 位加载更有效 -> vpshufb->vpermq如果数据在缓存中很热，这可能会成为 shuffle-port 吞吐量的瓶颈。Intel CPU 每个时钟周期最多可以执行 2 次加载 + 1 次存储（或在 IceLake 中，2 次加载 + 2 次存储）。如果没有内存瓶颈，我们可能会在前端出现瓶颈，因此在实践中不会使加载+存储和洗牌端口饱和。（https://agner.org/optimize/和https://uops.info/）
假设我们总是可以读取 32 个字节L而不会进入未映射的页面，这个函数也被简化了。但是在小 N 的 32 字节反转之后，输入的前 N 个字节成为 32 字节块中的最后 N 个字节。如果我们总是可以安全地在缓冲区末尾完成 32 字节的加载，那将是最方便的，但是在对象之前期望填充是不合理的。
#include <immintrin.h>
#include <stdalign.h>
#include <stddef.h>
#include <stdio.h>
#include <stdint.h>
// one vector of 32 bytes of digits, reversed and converted to ASCII
static inline
void ASCIIrev32B(void *dst, const void *src)
{
__m128i hi = _mm_loadu_si128(1 + (const __m128i*)src); // unaligned loads
__m128i lo = _mm_loadu_si128(src);
__m256i v = _mm256_set_m128i(lo, hi); // reverse 128-bit hi/lo halves
// compilers will hoist constants out of inline functions
__m128i byterev_lane = _mm_set_epi8(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15);
__m256i byterev = _mm256_broadcastsi128_si256(byterev_lane); // same in each lane
v = _mm256_shuffle_epi8(v, byterev); // in-lane reverse
v = _mm256_or_si256(v, _mm256_set1_epi8('0')); // digits to ASCII
_mm256_storeu_si256(dst, v); // Will usually be aligned in practice.
}
// Tested for N=32; could be bugs in the loop bounds for other N
// returns bytes written, like fwrite: N means no error, 0 means error in all fwrites
size_t LongNumPrint( uint8_t *num, size_t N)
{
// caller can print a name if it wants
const int revbufsize = 8192; // 8kiB on the stack should be fine
alignas(32) char revbuf[revbufsize];
if (N<32) {
// TODO: maybe use a smaller revbuf for this case to avoid touching new stack pages
ASCIIrev32B(revbuf, num); // the data we want is at the *end* of a 32-byte reverse
return fwrite(revbuf+32-N, 1, N, stdout);
}
size_t bytes_written = 0;
const uint8_t *inp = num+N; // start with last 32 bytes of num[]
do {
size_t chunksize = (inp - num >= revbufsize) ? revbufsize : inp - num;
const uint8_t *inp_stop = inp - chunksize + 32; // leave one full vector for the end
uint8_t *outp = revbuf;
while (inp > inp_stop) { // may run 0 times
inp -= 32;
ASCIIrev32B(outp, inp);
outp += 32;
}
// reverse first (lowest address) 32 bytes of this chunk of num
// into last 32 bytes of this chunk of revbuf
// if chunksize%32 != 0 this will overlap, which is fine.
ASCIIrev32B(revbuf + chunksize - 32, inp_stop - 32);
bytes_written += fwrite(revbuf, 1, chunksize, stdout);
inp = inp_stop - 32;
} while ( inp > num );
return bytes_written;
// caller can putchar('\n') if it wants
}
// wrapper that prints name and newline
void LongNumPrintName(uint8_t *num, size_t N, const char *name)
{
printf("%s:", name);
//LongNumPrint_scalar(num, N);
LongNumPrint(num, N);
putchar('\n');
}
// main() included on Godbolt link that runs successfully
这将编译并运行（在 Godbolt 上），gcc -O3 -march=haswell并为通过的 N=32 的标量循环产生相同的输出main。（我使用rand()了代替MyRandom()，所以我们可以使用您的 init 函数使用相同的种子进行测试并获得相同的数字。）
未经测试较大的 N，但 chunksize = min(ptrdiff, 8k) 的一般概念并使用它从末尾向下循环num[]应该是可靠的。
N%32如果我们转换第一个字节并将其传递给fwrite在开始主循环之前，我们可以加载（不仅仅是存储）对齐的向量。但这可能要么导致额外的write()系统调用，要么导致在 stdio 中进行笨重的复制。（除非已经有尚未打印的缓冲文本，例如Name:，在这种情况下，我们已经受到了惩罚。）
请注意，从技术上讲，inp从num. 因此inp -= 32，而不是inp = inp_stop-32将那个 UB 用于离开外循环的迭代。我实际上在这个版本中避免了这种情况，但它通常仍然有效，因为我认为 GCC 假设一个平面内存模型，并且 de-factor 定义了指针比较的行为。并且正常的操作系统保留零页，所以num绝对不能在物理内存开始的 32 个字节内（所以inp不能换行到高地址。）这一段大部分是我认为的第一次完全未经测试的尝试遗留下来的在内部循环中将指针递减到比实际更远的位置。
推荐阅读
- c++ - 使用 omp 进行并行化和同步
- android - 为什么在 Compose 中使用 jetpack 导航时视图一直闪烁？
- windows - 如何在这个混合的 cmd / powershell 代码中转义毒字符？
- javascript - jquery ajax 发送字符串数据和文件数据 PHP
- javascript - 未定义不和谐音乐机器人“连接”（discord.js）
- html - 如何通过 html 或 JavaScript 或 php 将文件下载到特定文件夹？
- operating-system - 虚拟地址与物理地址
- javascript - 如何从嵌套数组中取值
- javascript - 为什么我不能在 forEach() 中更改数组的元素？
- css - 页面向左移动