快速 1/X 除法（倒数）
Posted
技术标签:
【中文标题】快速 1/X 除法（倒数）【英文标题】：Fast 1/X division (reciprocal) 【发布时间】：2012-03-30 08:19:52 【问题描述】：如果精度不重要，是否有某种方法可以提高速度的倒数（X 的除法 1）？
所以，我需要计算 1/X。是否有一些解决方法让我失去精度但做得更快？
【问题讨论】：
这在很大程度上取决于您正在使用的硬件平台。此外，这还取决于您准备损失多少精度。显然，float recip(float x) return 1; 速度很快，但不是很准确……
Single-precision reciprocals run in 5 cycles on the lastest processors. A floating-point multiplication is also 5 cycles. 所以我严重怀疑你会比(float)1/(float)x 更快。
对于初学者，您的平台和编译器是什么？您正在处理什么样的数据？
@Mysticial 小心 5 个周期绝对是最佳情况下的最低延迟，但另一个数字是 37 个周期左右的最坏情况数字？请记住，硬件实现了迭代求根逼近算法，如牛顿法，直到精度足够
【参考方案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 时知道您可以只用一次乘法、减法和位移来近似倒数，那我就想知道了......嗯.. .......
至于基准测试，硬件 x2 指令结合硬件减法指令与现代计算机上的硬件 1.0/x 指令一样快（我的基准测试是在 Intel i7 上，但我会假设其他处理器的结果类似）。但是，如果将此算法作为新的汇编指令实现到硬件中，那么速度的提高可能足以使该指令非常实用。
有关此方法的更多信息，此实现基于精彩的"fast" inverse square root algorithm。
正如 Pharap 引起我注意的那样，从联合中读取非活动属性是未定义的行为，因此我从他有用的 cmets 中设计了两种可能的解决方案来避免未定义的行为。第一个解决方案似乎更像是一个令人讨厌的技巧，以绕过实际上并不比原始解决方案更好的语言语义。
#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];
第二种可能的解决方案更受欢迎，因为它完全摆脱了工会。但是，如果编译器没有正确优化，这个解决方案会慢很多。但是，从好的方面来说，下面的解决方案将完全不知道所提供的字节顺序：
-
字节宽度为 8 位
字节是目标机器上最小的原子单位。
双精度数为 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++ 的了解来吧。
【讨论】：
你可以解释一下这个幻数，以及它假设的浮点表示。 这很有趣。谢谢！你们有精度和速度对比测试的结果吗？ 您是否针对 x86 的近似倒数指令RCPSS 在您的 i7 上对此进行了测试？它与整数乘法一样快，并且不需要将数据从 XMM 寄存器移动到整数。您可以通过 _mm_rcp_ss( _mm_set_ss(x) ) 在 C++ 中使用它。如果您使用 -ffast-math，gcc 和 clang 会将 1.0/x 转换为 RCPSS + Newton-Raphson 迭代，但我认为如果您想要没有近似步骤的值，则必须手动使用内在函数。
像这样使用union 是未定义的行为。
@Pharap 您能否推断一下这个 UB 案例。我对 C++ 相当缺乏经验，并且在摆脱 UB 方面付出了巨大的努力。我会喜欢解释和/或在线资源来了解更多信息。非常感谢您提请我注意这个问题。【参考方案2】：
首先，确保这不是过早优化的情况。你知道这是你的瓶颈吗？
正如 Mystical 所说，1/x 可以很快计算出来。确保您没有为 1 或除数使用 double 数据类型。浮点数要快得多。
也就是说，基准，基准，基准。不要浪费时间在数值理论上花费数小时来发现性能不佳的根源是 IO 访问。
【讨论】：
“浮点数更快” - 真的吗？做出如此笼统的陈述是危险的。您可以做很多事情来更改编译器生成的代码。它还取决于编译器所针对的硬件。例如，在 IA32 上，gcc 在不使用 SSE 时生成的代码（我认为是 -mfpmath=387 选项）对于 double 和 float 的速度将是相同的，因为 FPU 只处理 80bit 值，任何速度差异都会下降到内存带宽。 是的，显然这是一个笼统的声明。但这个问题同样笼统。让 OP 给出细节，我就能做出更“雄辩”的回应。 1/x 可以快速计算出来......但是如何让编译器真正发出 RCPSS？ Mystical 在该声明中是“神秘的”，他引用了 5-clk 的下限，这实际上是反向吞吐量，应该是 6-clk 周期最佳情况场景另一端是最坏情况的情况是大约 38 时钟。 6-21,7-21,10-24,10-15,14-16,11-25,12-26,10-24,9-38,9-37,19,22,19,38,6 -38 在这些情况下，必须使用平均每指令时钟数 (CPI) 来比较苹果和橙子。 投了反对票，因为这个答案的大部分实际上并没有回答任何问题，甚至没有尝试过的微小部分。【参考方案3】：首先，如果您打开编译器优化，编译器可能会尽可能优化计算（例如，将其拉出循环）。要查看此优化，您需要在发布模式下构建和运行。
除法可能比乘法重（但评论者指出，倒数与现代 CPU 上的乘法一样快，在这种情况下，这不适合您的情况），所以如果您确实有 1/X 出现在循环内的某处（并且不止一次），您可以通过在循环内缓存结果（float Y = 1.0f/X;）然后使用Y 来提供帮助。 （编译器优化在任何情况下都可能这样做。）
此外，可以重新设计某些公式以消除除法或其他低效计算。为此，您可以发布正在执行的更大的计算。即便如此，有时也可以对程序或算法本身进行重组，以防止频繁触发耗时的循环。
可以牺牲多少准确性？如果您只需要一个数量级的机会，您可以使用模数运算符或按位运算轻松获得。
但是，一般来说，没有办法加速除法。如果有，编译器早就在做。
【讨论】：
如果您只需要一个数量级的机会，您可以使用模运算符或按位运算轻松获得。 如何？ 我并不是要暗示“微不足道”。另外，我应该添加 X >> 1 的警告（见评论结尾）。这种情况下，可以利用X^-1 = 2^(-log_2(X))，利用en.wikipedia.org/wiki/Find_first_set#Algorithms得到log_2(X)的数量级，得到数量级形式为2^-n。如果 X 的上限和下限已知，这可以用来提高速度。如果计算中的其他量（未在问题中显示）具有已知界限并且在数量级上有些相称，则可以将它们缩放并转换为整数。 如果您使用-ffast-math，编译器只能将 Y=1.0f/X 提升出循环，因此如果您不打算启用 @，最好在源代码中执行此操作987654326@ 告诉编译器你不关心舍入发生的地点/时间/方式，或以什么顺序。【参考方案4】：
我已经在 Arduino NANO 上测试了这些方法的速度和“准确性”。 基本计算是设置变量，Y = 15,000,000 和 Z = 65,535 （在我的真实情况下，Y 是一个常数，Z 可以在 65353 和 3000 之间变化，所以这是一个有用的测试） Arduino 上的计算时间是通过将引脚置低来建立的，然后在 calc 制作时将其置高，然后再次置低，并与逻辑分析仪比较时间。 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 Giffin's 浮点倒数算法：
Y * reciprocal(Z) takes 1.937msec which is a good improvement, but accuracy less so 213.68.
使用 nimig18's 浮动 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 Giffin 指点）
Y * Q*rsqrt(Z) takes 6.104 msec accuracy 228.116 with single iteration
All entertaining but ultimately disappointing!
【讨论】：
【参考方案5】：这应该通过许多预先展开的牛顿迭代来评估为霍纳多项式，该多项式使用大多数现代 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;
Fine Print：接近 0 时，近似值不太好，因此程序员需要测试性能或限制输入在诉诸硬件除法之前变低。 即负责！
【讨论】：
你能解释一下位模式 0x7EF127EA 的含义吗？我想是将 x 标准化为 ~ 1.0，但如何实现呢？【参考方案6】：据我所知，最快的方法是使用 SIMD 操作。 http://msdn.microsoft.com/en-us/library/796k1tty(v=vs.90).aspx
【讨论】：
还是买更快的cpu？:) 问题是算法问题。 RCPSS/ RCPPS 是个好建议。快速近似逆（和逆 sqrt）在 x86 上的硬件中可用，用于标量或浮点的 SIMD 向量。您不必在循环的其余部分使用 SIMD 即可利用。如果这个答案已经解释了其中的任何一个，它就不会让 cmets 如此困惑。以上是关于快速 1/X 除法（倒数）的主要内容，如果未能解决你的问题，请参考以下文章