快速 1/X 除法(倒数)

Fast 1/X division (reciprocal)

本文关键字:倒数 除法 快速      更新时间:2023-10-16

如果精度不是关键,有没有办法提高速度方面的倒数(X 除法 1(?

所以,我需要计算 1/X。是否有一些解决方法,让我失去精度但做得更快?

以下是更有效地近似

的方法

我相信他正在寻找一种更有效的近似 1.0/x 的方法,而不是一些近似的技术定义,即你可以使用 1 作为非常不精确的答案。我也认为这满足了这一点。

#ifdef __cplusplus
    #include <cstdint>
#else
    #include <stdint.h>
#endif
__inline__ double __attribute__((const)) reciprocal( double x ) {
    union {
        double dbl;
        #ifdef __cplusplus
            std::uint_least64_t ull;
        #else
            uint_least64_t ull;
        #endif
    } u;
    u.dbl = x;
    u.ull = ( 0xbfcdd6a18f6a6f52ULL - u.ull ) >> 1;
                                // pow( x, -0.5 )
    u.dbl *= u.dbl;             // pow( pow(x,-0.5), 2 ) = pow( x, -1 ) = 1.0 / x
    return u.dbl;
}
__inline__ float __attribute__((const)) reciprocal( float x ) {
    union {
        float single;
        #ifdef __cplusplus
            std::uint_least32_t uint;
        #else
            uint_least32_t uint;
        #endif
    } u;
    u.single = x;
    u.uint = ( 0xbe6eb3beU - u.uint ) >> 1;
                                // pow( x, -0.5 )
    u.single *= u.single;       // pow( pow(x,-0.5), 2 ) = pow( x, -1 ) = 1.0 / x
    return u.single;
}

嗯。。。。。。。如果 CPU 制造商知道在设计 CPU 时只需一个乘法、减法和位移就可以近似倒数,我就知道......嗯。。。。。。。。。

至于基准测试,硬件 x2 指令与硬件减法指令相结合的速度与现代计算机上的硬件 1.0/x 指令一样快(我的基准测试是在英特尔 i7 上,但我假设其他处理器也有类似的结果(。但是,如果将此算法作为新的汇编指令实现到硬件中,那么速度的增加可能足以使该指令非常实用。

有关此方法的更多信息,此实现基于精彩的"快速"平方根反比算法。

正如 Pharap 提请我注意的那样,从工会中读取非活动属性是未定义的行为,因此我从他的有用评论中设计了两种可能的解决方案来避免未定义的行为。第一个解决方案似乎更像是一个令人讨厌的技巧,可以绕过实际上并不比原始解决方案更好的语言语义。

#ifdef __cplusplus
    #include <cstdint>
#else
    #include <stdint.h>
#endif
__inline__ double __attribute__((const)) reciprocal( double x ) {
    union {
        double dbl[2];
        #ifdef __cplusplus
            std::uint_least64_t ull[2];
        #else
            uint_least64_t ull[2];
        #endif
    } u;
    u.dbl[0] = x; // dbl is now the active property, so only dbl can be read now
    u.ull[1] = 0;//trick to set ull to the active property so that ull can be read
    u.ull][0] = ( 0xbfcdd6a18f6a6f52ULL - u.ull[0] ) >> 1;
    u.dbl[1] = 0; // now set dbl to the active property so that it can be read
    u.dbl[0] *= u.dbl[0];
    return u.dbl[0];
}
__inline__ float __attribute__((const)) reciprocal( float x ) {
    union {
        float flt[2];
        #ifdef __cplusplus
            std::uint_least32_t ull[2];
        #else
            uint_least32_t ull[2];
        #endif
    } u;
    u.flt[0] = x; // now flt is active
    u.uint[1] = 0; // set uint to be active for reading and writing
    u.uint[0] = ( 0xbe6eb3beU - u.uint[0] ) >> 1;
    u.flt[1] = 0; // set flt to be active for reading and writing
    u.flt[0] *= u.flt[0];
    return u.flt[0];
}

第二种可能的解决方案更可口,因为它完全摆脱了工会。但是,如果编译器未正确优化此解决方案,则此解决方案将慢得多。但是,从好的方面来说,下面的解决方案将完全与提供的字节顺序无关:

  1. 字节宽度为 8 位
  2. 该字节是目标计算机上最小的原子单元。
  3. 双精度为 8 字节宽,浮点数宽为 4 字节。

#ifdef __cplusplus
    #include <cstdint>
    #include <cstring>
    #define stdIntWithEightBits std::uint8_t
    #define stdIntSizeOfFloat std::uint32_t
    #define stdIntSizeOfDouble std::uint64_t
#else
    #include <stdint.h>
    #include <string.h>
    #define stdIntWithEightBits uint8_t
    #define stdIntSizeOfFloat uint32_t
    #define stdIntSizeOfDouble uint64_t
#endif

__inline__ double __attribute__((const)) reciprocal( double x ) {
    double byteIndexFloat = 1.1212798184631136e-308;//00 08 10 18 20 28 30 38 bits
    stdIntWithEightBits* byteIndexs = reinterpret_cast<stdIntWithEightBits*>(&byteIndexFloat);
    
    stdIntWithEightBits* inputBytes = reinterpret_cast<stdIntWithEightBits*>(&x);
    
    stdIntSizeOfDouble inputAsUll = (
        (inputBytes[0] << byteIndexs[0]) |
        (inputBytes[1] << byteIndexs[1]) |
        (inputBytes[2] << byteIndexs[2]) |
        (inputBytes[3] << byteIndexs[3]) |
        (inputBytes[4] << byteIndexs[4]) |
        (inputBytes[5] << byteIndexs[5]) |
        (inputBytes[6] << byteIndexs[6]) |
        (inputBytes[7] << byteIndexs[7])
    );
    inputAsUll = ( 0xbfcdd6a18f6a6f52ULL - inputAsUll ) >> 1;
    
    double outputDouble;
    
    const stdIntWithEightBits outputBytes[] = {
        inputAsUll >> byteIndexs[0],
        inputAsUll >> byteIndexs[1],
        inputAsUll >> byteIndexs[2],
        inputAsUll >> byteIndexs[3],
        inputAsUll >> byteIndexs[4],
        inputAsUll >> byteIndexs[5],
        inputAsUll >> byteIndexs[6],
        inputAsUll >> byteIndexs[7]
    };
    memcpy(&outputDouble, &outputBytes, 8);
    
    return outputDouble * outputDouble;
}

__inline__ float __attribute__((const)) reciprocal( float x ) {
    float byteIndexFloat = 7.40457e-40; // 0x00 08 10 18 bits
    stdIntWithEightBits* byteIndexs = reinterpret_cast<stdIntWithEightBits*>(&byteIndexFloat);
    
    stdIntWithEightBits* inputBytes = reinterpret_cast<stdIntWithEightBits*>(&x);
    
    stdIntSizeOfFloat inputAsInt = (
        (inputBytes[0] << byteIndexs[0]) |
        (inputBytes[1] << byteIndexs[1]) |
        (inputBytes[2] << byteIndexs[2]) |
        (inputBytes[3] << byteIndexs[3])
    );
    inputAsInt = ( 0xbe6eb3beU - inputAsInt ) >> 1;
    
    float outputFloat;
    
    const stdIntWithEightBits outputBytes[] = {
        inputAsInt >> byteIndexs[0],
        inputAsInt >> byteIndexs[1],
        inputAsInt >> byteIndexs[2],
        inputAsInt >> byteIndexs[3]
    };
    memcpy(&outputFloat, &outputBytes, 4);
    
    return outputFloat * outputFloat;
}

免责声明:最后,请注意,我更像是C++新手。因此,我张开双臂欢迎任何最佳实践、适当的格式或含义清晰的编辑,以改善所有阅读它的人的答案的质量,并在未来的岁月里扩展我对C++的了解。

首先,确保这不是过早优化的情况。你知道这是你的瓶颈吗?

正如Mystical所说,1/x可以非常快速地计算出来。确保没有对 1 或除数使用double数据类型。浮子要快得多。

也就是说,基准

,基准,基准。不要浪费时间在数值理论上花费数小时,只是为了发现性能不佳的根源是 IO 访问。

首先,如果打开编译器优化,编译器可能会在可能的情况下优化计算(例如,将其拉出循环(。 若要查看此优化,需要在发布模式下生成并运行。

除法可能比乘法重(但一位评论者指出,在现代 CPU 上,倒数与乘法一样快,在这种情况下,这对您的情况是不正确的(,所以如果你确实有1/X出现在循环中的某个地方(并且不止一次(,你可以通过在循环内缓存结果来提供帮助( float Y = 1.0f/X; (,然后使用 Y . (编译器优化在任何情况下都可能执行此操作。

此外,可以重新设计某些公式以删除除法或其他低效计算。 为此,您可以发布正在执行的较大计算。 即使在那里,程序或算法本身有时也可以重组,以防止频繁地命中耗时的循环。

可以牺牲多少精度? 如果很有可能你只需要一个数量级,你可以使用模运算符或按位运算轻松获得。

但是,一般来说,没有办法加快分裂。 如果有的话,编译器早就在这样做了。

我已经在Arduino NANO上测试了这些方法的速度和"准确性".
基本计算是设置变量,Y = 15,000,000 和 Z = 65,535
(在我的真实情况下,Y 是一个常数,Z 可以在 65353 和 3000 之间变化,所以这是一个有用的测试(
Arduino上的计算时间是通过将引脚放低,然后在计算时高,然后再次降低并与逻辑分析仪进行比较来确定的。 100 个周期。将变量作为无符号整数:-

Y * Z takes 0.231 msec
Y / Z takes  3.867 msec.  
With variables as floats:-  
Y * Z takes  1.066 msec
Y / Z takes  4.113 msec.  
Basic Bench Mark  and ( 15,000,000/65535 = 228.885 via calculator.) 

使用 {Jack G's} 浮点倒数算法:

Y * reciprocal(Z)  takes  1.937msec  which is a good improvement, but accuracy less so 213.68.  

使用 {nimig18's} float inv_fast:

Y* inv_fast(Z)  takes  5.501 msec  accuracy 228.116  with single iteration  
Y* inv_fast(Z)  takes  7.895 msec  accuracy 228.883  with second iteration 

使用维基百科的Q_rsqrt(由{Jack G}指向(

Y * Q*rsqrt(Z) takes  6.104 msec  accuracy   228.116  with single iteration  
All entertaining but ultimately disappointing!
<</div> div class="answers">

这应该通过许多预展开的牛顿迭代来完成,这些迭代被评估为霍纳多项式,该多项式使用融合乘法累加运算,大多数现代 CPU 在单个 Clk 周期中执行(每次(:

float inv_fast(float x) {
    union { float f; int i; } v;
    float w, sx;
    int m;
    sx = (x < 0) ? -1:1;
    x = sx * x;
    v.i = (int)(0x7EF127EA - *(uint32_t *)&x);
    w = x * v.f;
    // Efficient Iterative Approximation Improvement in horner polynomial form.
    v.f = v.f * (2 - w);     // Single iteration, Err = -3.36e-3 * 2^(-flr(log2(x)))
    // v.f = v.f * ( 4 + w * (-6 + w * (4 - w)));  // Second iteration, Err = -1.13e-5 * 2^(-flr(log2(x)))
    // v.f = v.f * (8 + w * (-28 + w * (56 + w * (-70 + w *(56 + w * (-28 + w * (8 - w)))))));  // Third Iteration, Err = +-6.8e-8 *  2^(-flr(log2(x)))
    return v.f * sx;
}

细则:接近 0,近似值表现不佳,因此程序员需要测试性能或在诉诸硬件划分之前限制输入达到低电平。即负责任!

据我所知,最快的方法是使用 SIMD 操作。 http://msdn.microsoft.com/en-us/library/796k1tty(v=vs.90(.aspx

rcpss 汇编指令计算近似倒数 |相对误差|≤ 1.5 ∗ 2^−12。

您可以在带有 -mrecip 标志的编译器上启用它(您可能还需要 -ffast-math (。

内在是_mm_rcp_ss.