什么是最大的"非浮动"整数,可以存储在IEEE 754 double类型而不会丢失精度？
我试图找到潜在因子素数的除数(形式n!+ - 1的数量),因为我最近买了Skylake-X工作站,我认为我可以使用AVX512指令加快速度.
算法简单,主要步骤是对同一个除数重复取模.主要是循环大范围的n值.这是用c写的天真的方法(P是素数表):
uint64_t factorial_naive(uint64_t const nmin, uint64_t const nmax, const uint64_t *restrict P)
{
uint64_t n, i, residue;
for (i = 0; i < APP_BUFLEN; i++){
residue = 2;
for (n=3; n <= nmax; n++){
residue *= n;
residue %= P[i];
// Lets check if we found factor
if (nmin <= n){
if( residue == 1){
report_factor(n, -1, P[i]);
}
if(residue == P[i]- 1){
report_factor(n, 1, P[i]);
}
}
}
}
return EXIT_SUCCESS;
}
这里的想法是针对同一组除数检查大范围的n,例如1,000,000 - > 10,000,000.因此,我们将以数百万次对相同的除数进行模数尊重.使用DIV非常慢,因此根据计算范围有几种可能的方法.在我的情况下,n很可能小于10 …
SSE2具有在单精度浮点数和32位整数之间转换向量的指令.
_mm_cvtps_epi32()_mm_cvtepi32_ps()但是没有双精度和64位整数的等价物.换句话说,他们失踪了:
_mm_cvtpd_epi64()_mm_cvtepi64_pd()似乎AVX也没有它们.
模拟这些内在函数的最有效方法是什么？
我想矢量化两个内存对齐数组的乘法.我没有找到任何方法在AVX/AVX2中乘以64*64位,所以我只是循环展开和AVX2加载/存储.有更快的方法吗？
注意:我不想保存每次乘法的高半结果.
void multiply_vex(long *Gi_vec, long q, long *Gj_vec){
int i;
__m256i data_j, data_i;
__uint64_t *ptr_J = (__uint64_t*)&data_j;
__uint64_t *ptr_I = (__uint64_t*)&data_i;
for (i=0; i<BASE_VEX_STOP; i+=4) {
data_i = _mm256_load_si256((__m256i*)&Gi_vec[i]);
data_j = _mm256_load_si256((__m256i*)&Gj_vec[i]);
ptr_I[0] -= ptr_J[0] * q;
ptr_I[1] -= ptr_J[1] * q;
ptr_I[2] -= ptr_J[2] * q;
ptr_I[3] -= ptr_J[3] * q;
_mm256_store_si256((__m256i*)&Gi_vec[i], data_i);
}
for (; i<BASE_DIMENSION; i++)
Gi_vec[i] -= Gj_vec[i] * q;
}
更新:
我正在使用Haswell微体系结构和ICC/GCC编译器.所以AVX和AVX2都很好.我在乘法循环展开后-=用C inrisic 替换_mm256_sub_epi64它,在那里得到一些加速.目前,它是ptr_J[0] *= q; ...
我用, …
我们是否仍然需要在软件中模拟128位整数,或者现在平均桌面处理器中是否有硬件支持？
SSE/AVX寄存器可以被视为整数或浮点BigNums.也就是说,人们可以忽视存在通道.是否有一种简单的方法可以利用这种观点并将这些寄存器单独或组合用作BigNum？我问,因为我从BigNum库中看到的很少,它们几乎普遍存储并对数组进行算术运算,而不是SSE/AVX寄存器.可移植性？
例:
假设您将SSE寄存器的内容存储为a中的键std::set,您可以将这些内容作为BigNum进行比较.
我创建了一个使用SIMD进行64位*64位到128位的功能.目前我已经使用SSE2(acutally SSE4.1)实现了它.这意味着它可以同时运行两个64b*64b到128b的产品.同样的想法可以扩展到AVX2或AVX512,同时提供四个或八个64b*64到128b的产品.我的算法基于http://www.hackersdelight.org/hdcodetxt/muldws.c.txt
该算法进行一次无符号乘法,一次有符号乘法和两次有符号*无符号乘法.签名的*signed和unsigned*unsigned操作很容易使用_mm_mul_epi32和_mm_mul_epu32.但混合签名和未签名的产品给我带来了麻烦.例如,考虑一下.
int32_t x = 0x80000000;
uint32_t y = 0x7fffffff;
int64_t z = (int64_t)x*y;
双字产品应该是0xc000000080000000.但是如果你假设你的编译器知道如何处理混合类型,你怎么能得到这个呢？这就是我想出的:
int64_t sign = x<0; sign*=-1; //get the sign and make it all ones
uint32_t t = abs(x); //if x<0 take two's complement again
uint64_t prod = (uint64_t)t*y; //unsigned product
int64_t z = (prod ^ sign) - sign; //take two's complement based on the sign
使用SSE可以这样做
__m128i xh; //(xl2, xh2, xl1, xh1) high is signed, low unsigned
__m128i …当我第一次使用Haswell处理器时,我尝试使用FMA来确定Mandelbrot集.主要算法是这样的:
intn = 0;
for(int32_t i=0; i<maxiter; i++) {
floatn x2 = square(x), y2 = square(y); //square(x) = x*x
floatn r2 = x2 + y2;
booln mask = r2<cut; //booln is in the float domain non integer domain
if(!horizontal_or(mask)) break; //_mm256_testz_pd(mask)
n -= mask
floatn t = x*y; mul2(t); //mul2(t): t*=2
x = x2 - y2 + cx;
y = t + cy;
}
这确定n像素是否在Mandelbrot集中.因此对于双浮点,它运行超过4个像素(floatn = __m256d,intn = __m256i).这需要4个SIMD浮点乘法和4个SIMD浮点加法.
然后我修改了这个就像这样使用FMA
intn n = 0; …在Hacker的喜悦中,有一种算法来计算两个(带符号)单词的双字产品.
该函数muldws1使用四次乘法和五次加法来计算两个单词的双字.
在该代码的末尾有一行注释掉
/* w[1] = u*v; // Alternative. */
该替代方案使用五次乘法和四次加法,即它为乘法交换加法.
但我认为这种替代方法可以改进.我还没有说过硬件.让我们假设一个假设的CPU,它可以计算两个字但不是高位字的乘积的低位字(例如,对于32位字32x32到低32).在这种情况下,在我看来,这个算法可以改进.这是我假设32位字的想法(相同的概念适用于64位字).
void muldws1_improved(int w[], int32_t x, int32_t y) {
uint16_t xl = x; int16_t xh = x >> 16;
uint16_t yl = y; int16_t yh = y >> 16;
uint32 lo = x*y;
int32_t t = xl*yh + xh*yl;
uint16_t tl = t; int16_t th = t >>16;
uint16_t loh = lo >> 16;
int32_t cy = loh<tl; //carry
int32_t hi = xh*yh + th …从该值我们可以推断它使用与双精度浮点硬件相同的组件。但是 double 有 53 位有效数，那么为什么 AVX512-IFMA 限制为 52 位呢？当然尾数只有 52 位并且隐藏了一位，但它仍然对值有贡献，需要输入加法器/乘法器/除法器......
我在AVX2上工作,需要计算64位x64位 - > 128位加宽乘法,并以最快的方式获得64位高位.由于AVX2没有这样的指令,使用Karatsuba算法提高效率和提高速度是否合理？