快速 1/X 除法（倒数）

Question

如果精度不重要，是否有某种方法可以提高速度的倒数（X 上的除法 1）？

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

原文由 klm123 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

𝗛𝗲𝗿𝗲’𝘀𝗧𝗼𝗔𝗽𝗽𝗿𝗼𝘅𝗶𝗺𝗮𝘁𝗲𝗠𝗼𝗿𝗲𝗘𝗳𝗳𝗶𝗰𝗶𝗲𝗻𝘁𝗹𝘆

我相信他正在寻找的是一种更有效的近似 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 时知道您可以只用一次乘法、减法和位移来近似倒数，那我就想知道了…. 嗯…….. .

至于基准测试，硬件 x 2指令结合硬件减法指令与现代计算机上的硬件 1.0/x 指令一样快（我的基准测试是在 Intel 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++ 的了解来。

原文由 Jack G 发布，翻译遵循 CC BY-SA 4.0 许可协议