如何将代码传递给AVX2代码并获得与以前相同的结果？
是否可以在LongNumInit、LongNumPrint函数中使用__m256i，而不是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实用，特别是部分单词算术，其中您的临时变量可能不是"；"归一化"；，让你懒散地搬运。
但很明显，您只是在问这个代码、随机初始化和打印函数，而不是如何在这种格式的两个数字之间进行数学运算
我们可以很好地将这两者矢量化。我的LongNumPrintName()是你的临时替代品。
对于LongNumInit，我只是展示了一个构建块，它存储两个32字节的块，并返回一个递增的指针。循环调用它。(每次调用自然会产生2个矢量，因此对于小N，您可以制作备用版本。)
LongNumInit
生成包含随机数字的1GB文本文件的最快方法是什么？在4GHz Skylake上以大约33GB/s的速度生成空间分隔的随机ASCII十进制数字，包括write()系统调用/dev/null的开销。(这高于DRAM带宽；128kiB的缓存阻塞可以让存储命中二级缓存。/dev/null的内核驱动程序甚至不读取用户空间缓冲区。)
它可以很容易地适应于void LongNumInit(uint8_t *L, size_t N )的AVX2版本。我的答案是使用AVX2 xorshift128+PRNG(在__m256i的64位元素中用4个独立的PRNG进行矢量化)，类似于xorshift228+的AVX/SSE版本。这应该与您的rand() % 10具有相似的随机性。
它通过乘法逆将其分解为十进制数字，用移位和vpmulhuw除以10并取模，使用"为什么GCC在实现整数除法时使用奇异数乘法？"？。(实际上，使用GNUC原生向量语法，让GCC确定魔术常数，并发出乘法和移位，以获得方便的语法，如v16u dig1 = v % ten;和v /= ten;)
您可以使用_mm256_packus_epi16将两个16位数字的矢量打包为8位元素，而不是将奇数元素转换为ASCII' '，将偶数元素转换为ASCII'0'..'9'。(因此，将vec_store_digit_and_space更改为打包向量对，而不是与常数进行"或"运算，见下文)
使用gcc、clang或ICC(或者任何其他理解C99的GNU C方言和英特尔内部语言的编译器)进行编译。
请参阅https://gcc.gnu.org/onlinedocs/gcc/Vector-Extensions.html对于__attribute__((vector_size(32)))部分，以及https://software.intel.com/sites/landingpage/IntrinsicsGuide/对于CCD_ 20的东西。而且https://stackoverflow.com/tags/sse/info.
#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和glibc缺少aligned_realloc和aligned_calloc。MSVC确实提供了这些：为什么没有'；aligned_realloc'；在大多数平台上？允许您有时在对齐的缓冲区的末尾分配更多的空间而不复制它；try_ alloc"；对于C++来说是理想的选择，因为如果不常见的可复制对象更改了地址，则可能需要运行复制构造函数。有时强制进行不必要的复制的非表达型分配器API是我的一大苦恼。
LongNumPrint
采用uint8_t *Name参数是糟糕的设计。如果调用者想先打印一个"something:"字符串，他们可以这样做。您的函数应该只做printf"%d"对int所做的操作。
由于您以反向打印顺序存储数字，因此需要将字节反向存储到tmp缓冲区中，并通过与'0'进行"或"运算将0..9字节值转换为'0'..'9'ASCII字符值。然后将该缓冲区传递给fwrite。
具体来说，使用alignas(32) char tmpbuf[8192];作为局部变量。
您可以在固定大小的块(如1kiB或8kiB)中工作，而不是分配潜在的巨大缓冲区。您可能仍然希望通过stdio(而不是直接使用write()并管理自己的I/O缓冲)。对于8kiB缓冲区，高效的fwrite可能只是将其直接传递给write()，而不是将memcpy传递到stdio缓冲区。您可能想对此进行调整，但将tmp缓冲区保持在L1d缓存的一半以下将意味着在写入后重新读取时，它在缓存中仍然很热。
缓存阻塞使循环边界更加复杂，但对于非常大的N.来说，这是值得的
字节一次反转32个字节：
您可以通过决定您的数字按第一顺序存储在MSD中来避免这项工作，但如果您确实想实现加法，则必须从末尾向后循环。
您的函数可以用SIMD_mm_shuffle_epi8来实现，以反转16字节的块，从数字数组的末尾开始，并写入tmp缓冲区的开头。
或者更好的方法是，加载vmovdqu/vinserti12816字节负载，以将_mm256_shuffle_epi8馈送到通道内的字节反向，设置32字节存储。
在英特尔CPU上，vinserti128解码为加载+ALU uop，但它可以在任何矢量ALU端口上运行，而不仅仅是混洗端口。因此，两个128位加载比256位加载更有效->CCD_ 44->vpermq，如果数据在缓存中很热，它可能会成为混洗端口吞吐量的瓶颈。英特尔CPU每个时钟周期最多可以进行2次加载+1次存储(或者在冰湖，2次加载+2次存储)。如果没有内存瓶颈，我们可能会在前端遇到瓶颈，所以在实践中不会饱和加载+存储和混洗端口。(https://agner.org/optimize/和https://uops.info/)
这个函数也被简化了，因为我们总是可以从L读取32个字节，而不会进入未映射的页面。但是，在对小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
它使用gcc -O3 -march=haswell编译和运行(在Godbolt上)，并为main传递的N=32生成与标量循环相同的输出。(我使用了rand()而不是MyRandom()，所以我们可以使用init函数使用相同的种子进行测试并获得相同的数字。)
对于较大的N未进行测试，但chunksize=min(ptrdiff，8k)并使用它从num[]的末尾向下循环的总体想法应该是可靠的。
如果我们转换了第一个N%32字节，并在开始主循环之前将其传递给fwrite，那么我们可以加载(而不仅仅是存储)对齐的向量。但这可能会导致额外的write()系统调用，或者导致stdio内部笨拙的复制。(除非已经有缓冲文本尚未打印，如Name:，在这种情况下，我们已经有了惩罚。)
请注意，从技术上讲，在num开始之后递减inp是C UB。因此inp -= 32而不是inp = inp_stop-32将具有用于离开外循环的迭代的UB。在这个版本中，我实际上避免了这种情况，但它通常还是有效的，因为我认为GCC假设了一个平面内存模型，并且去因子定义了指针比较的行为。普通操作系统保留零页，因此num肯定不能在物理内存开始的32字节内(因此inp不能换行到高地址。