c++ - 解释了一种将双精度数舍入为 32 位整数的快速方法

Question

在阅读Lua的源代码时，我注意到 Lua 使用宏将double值四舍五入为 32 位int值。该宏在Llimits.h头文件中定义，内容如下：

union i_cast {double d; int i[2]};
#define double2int(i, d, t) \
    {volatile union i_cast u; u.d = (d) + 6755399441055744.0; \
    (i) = (t)u.i[ENDIANLOC];}

这里ENDIANLOC是根据字节序定义的：0代表小端，1代表大端架构；Lua 小心地处理字节序。该t 参数被替换为整数类型，如intor unsigned int。

我做了一些研究，发现该宏有一种更简单的格式，它使用相同的技术：

#define double2int(i, d) \
    {double t = ((d) + 6755399441055744.0); i = *((int *)(&t));}

或者，在 C++ 风格中：

inline int double2int(double d)
{
    d += 6755399441055744.0;
    return reinterpret_cast<int&>(d);
}

这个技巧可以在任何使用IEEE 754的机器上运行（这意味着今天几乎每台机器）。它适用于正数和负数，并且四舍五入遵循银行家规则。（这并不奇怪，因为它遵循 IEEE 754。）

我写了一个小程序来测试它：

int main()
{
    double d = -12345678.9;
    int i;
    double2int(i, d)
    printf("%d\n", i);
    return 0;
}

它-12345679按预期输出 。

我想详细了解这个棘手的宏是如何工作的。幻数6755399441055744.0实际上是 2 ⁵¹  + 2 ⁵²，即 1.5 × 2 ⁵²，二进制的 1.5 可以表示为 1.1。当任何 32 位整数与这个幻数相加时——</p>

好吧，我从这里迷路了。这个技巧是如何工作的？

更新

1. 正如@Mysticial 指出的那样，这种方法并不局限于 32-bit ，只要数字在 2 ⁵²int范围内，它也可以扩展到 64-bit 。（虽然宏需要一些修改。）int

2. 有些资料说这种方法不能在Direct3D中使用。

3. 当使用 Microsoft assembler for x86 时，有一个用汇编代码编写的更快的宏（以下也是从 Lua 源代码中提取的）：

    #define double2int(i,n)  __asm {__asm fld n   __asm fistp i}
   

4. 单精度数有一个类似的幻数： 1.5 × 2 ²³。

Answer 1

double浮点类型的值表示如下：

它可以看作是两个 32 位整数；现在，int所有版本的代码（假设它是 32-bit int）都是图中右侧的那个，所以你最终所做的只是取最低的 32 位尾数。

---

现在，到神奇的数字；正如您所说， 6755399441055744 是 2 ⁵¹  + 2 ⁵²；添加这样一个数字会强制double进入 2 ⁵²和 2 ⁵³之间的“甜蜜范围” ，正如Wikipedia 所解释的那样，它有一个有趣的属性：

在 2 ⁵²  = 4,503,599,627,370,496 和 2 ⁵³  = 9,007,199,254,740,992 之间，可表示的数字恰好是整数。

这是因为尾数是 52 位宽。

^{关于添加 2 51}  + 2 ⁵²的另一个有趣的事实是，它只影响尾数的两个最高位——无论如何都会被丢弃，因为我们只取它的最低 32 位。

---

最后但并非最不重要的：标志。

IEEE 754 浮点使用幅度和符号表示，而“普通”机器上的整数使用 2 的补码算法；这是如何处理的？

我们只讨论了正整数；现在假设我们正在处理由 32 位 表示的范围内的负数int，因此（绝对值）小于（-2 ³¹  + 1）；称之为-a。通过添加幻数显然可以使这样的数字变为正数，结果值为 2 ⁵²  + 2 ⁵¹  + (-a)。

现在，如果我们用 2 的补码表示来解释尾数，我们会得到什么？它必须是 (2 ⁵²  + 2 ⁵¹ ) 和 (−a) 的 2 补码和的结果。同样，第一项仅影响高两位，位 0-50 中剩下的是 (-a) 的 2 的补码表示（同样，减去高两位）。

由于将 2 的补码减少到更小的宽度只是通过切除左侧的额外位来完成，因此采用低 32 位可以在 32 位 2 的补码算术中正确地为我们提供 (-a)。

Answer 2

这种“技巧”来自较旧的 x86 处理器，使用 8087 指令/接口进行浮点运算。在这些机器上，有一条将浮点数转换为整数“拳头”的指令，但它使用当前的 fp 舍入模式。不幸的是，C 规范要求 fp->int 转换向零截断，而所有其他 fp 操作舍入到最近，因此进行
fp->int 转换需要首先更改 fp 舍入模式，然后进行拳头，然后恢复 fp舍入模式。

现在在最初的 8086/8087 上，这并不算太糟糕，但在后来开始获得超标量和乱序执行的处理器上，改变 fp 舍入模式通常会序列化 CPU 内核并且非常昂贵。所以在像 Pentium-III 或 Pentium-IV 这样的 CPU 上，这个总成本相当高——正常的 fp->int 转换比这个 add+store+load 技巧要贵 10 倍或更多。

然而，在 x86-64 上，浮点是使用 xmm 指令完成的，转换
fp->int 的成本非常小，因此这种“优化”可能比正常转换慢。

Answer 3

如果这有助于可视化，那么 lua 魔法值

  (2^52+2^51, or base2 of 110 then [50 zeros]

十六进制

  0x  0018 0000 0000 0000 (18e12)

八进制

  0 300 00000 00000 00000 ( 3e17)

Answer 4

下面是上述 Lua 技巧的一个更简单的实现：

/**
 * Round to the nearest integer.
 * for tie-breaks: round half to even (bankers' rounding)
 * Only works for inputs in the range: [-2^51, 2^51]
 */
inline double rint(double d)
{
    double x = 6755399441055744.0;  // 2^51 + 2^52
    return d + x - x;
}

该技巧适用于绝对值 < 2 ^ 51 的数字。

这是一个测试它的小程序：ideone.com

#include <cstdio>

int main()
{
    // round to nearest integer
    printf("%.1f, %.1f\n", rint(-12345678.3), rint(-12345678.9));

    // test tie-breaking rule
    printf("%.1f, %.1f, %.1f, %.1f\n", rint(-24.5), rint(-23.5), rint(23.5), rint(24.5));      
    return 0;
}

// output:
// -12345678.0, -12345679.0
// -24.0, -24.0, 24.0, 24.0