典型浮点数中没有倒数的最小正整数是多少？

Question

一个常见的假设是1 / x * x == 1。在常见的符合 IEEE 754 的硬件上打破这一点的最小正整数是多少？

当乘法逆的假设失败时，写得不好的有理算术就不再起作用了。因为包括 C 和 C++ 在内的许多语言默认使用舍入到零将浮点数转换为整数，所以即使是很小的错误也可能导致整数结果偏位。

快速测试程序会产生各种结果。

#include <iostream>

int main () {
    {
        double n;
        for ( n = 2; 1 / n * n == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
        for ( ; (int) ( 1 / n * n ) == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
    }
    {
        float n;
        for ( n = 2; 1 / n * n == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
        for ( ; (int) ( 1 / n * n ) == 1; ++ n ) ;
        std::cout << n << " (" << 1 - 1/n*n << ")\n";
    }
}

在 ideone.com 上使用 GCC 4.3.4 结果是

41 (5.42101e-20)
45 (5.42101e-20)
41 (5.42101e-20)
45 (5.42101e-20)

使用 GCC 4.5.1 会产生相同的结果，但报告的误差范围恰好为零。

在我的机器上（GCC 4.7.2 或 Clang 4.1），结果是

49 (1.11022e-16)
49 (1.11022e-16)
41 (5.96046e-08)
41 (5.96046e-08)

这与--fast-math 选项无关。使用-mfpmath=387 出乎意料地产生了

41 (5.42101e-20)
41 (5.42101e-20)
41 (5.42101e-20)
41 (5.42101e-20)

值 5×10^-20 似乎暗示 epsilon 对应于 64 位尾数，即使用 Intel 80 位扩展精度的内部计算。

这似乎高度依赖于 FPU 硬件。是否有适合测试的可靠值？

注意：我不在乎语言标准或编译器对浮点数系统的保证，尽管我认为在任何常见的编程系统中都没有很多有意义的保证。我想知道数字和现实世界计算机之间的交互。

Answer 1

双精度：

1/41 = 0x1.8f9c18f9c18fap-6 和 41*0x1.8f9c18f9c18fap-6 = 0x1.000000000000028，四舍五入为 1。 1/45 = 0x1.6c16c16c16c17p-6，45*0x1.6c16c16c16c17p-6 = 0x1.00000000000002c，取整为 1。

然而，

1/49 = 0x1.4e5e0a72f0539p-6，49*0x1.4e5e0a72f0539p-6 = 0x0.fffffffffffffa4，四舍五入为0x0.ffffffffffff8 = 0x1.ffffffffffff0p-1

49 确实有倒数！它是 0x1.4e5e0a72f053ap-6。

更一般地，如果 f 是 [1, 2) 中的浮点数，则 f 有倒数。在通常的四舍五入算术中，如果一个数字位于 [1 - 2^-54, 1 + 2^-53] 中，则该数字将四舍五入为 1。 请注意，最接近 1/f 的双精度数，例如 d，距离 1/f 小于 2^-54。如果 d > 1/f，那么我们就是黄金； 1 -54) -54 * f -53，因此 f*d 舍入为 1。如果 d -53。如果是，则 f*d 位于 [1 - 2^-53, 1 - 2^-54)。如果你取 e = 2^-53 + d，那么 e*f > 1 并且 e*f = d*f + 2^-53*f -53 + 2^-52 = 1 + 2^-53，再次四舍五入为 1。

编辑：上面的推理是错误的，因为两个连续双打之间的步幅相差了两倍。没有倒数的双精度示例是 0x1.ffffffbffffffe。 0x1.0000002000001p-1 太小，但 0x1.0000002000002p-1 太大。没有倒数的整数的最小示例是 237。1/237 大约是 0x1.1485f0e0acd3B68c6Bp-8，四舍五入为 0x1.1485f0e0acd58p-8。这个数字太小了，而它之后的下一个 double 太大了。

Answer 2

这个问题似乎与 C++ 选择转换为整数的方法有关。

这是一个用于比较的 Ada 版本，测试 32 位、64 位和 80 位浮点数（只要求 7、15 和 18 位数字，或者使用前两个的内置类型）。

结果和注释第一，代码如下。

$ gnatmake fp_torture.adb
gcc -c fp_torture.adb
gnatbind -x fp_torture.ali
gnatlink fp_torture.ali
$ ./fp_torture
 41 ( 5.96046E-08)
Error representing float  2.14748E+09 as integer
 49 ( 1.11022302462516E-16)
 2147483647 ( 0.00000000000000E+00)
 41 ( 5.42101086242752217E-20)
 2147483647 ( 0.00000000000000000E+00)
$

正如我们所见，浮点计算重现了 C++ 故障点，并确认使用了 387 个 80 位浮点数。但是将（非常接近 1 的数字）转换回整数，比较有效。

看到这一点，在 C++ 示例中添加适当的舍入确实可以进行比较。在 MAX_INT 处添加终止条件，“double n”然后起作用。

当++n 未能递增n 时，“float n”中出现了一点，因此迭代器停止迭代，但这是另一回事！

下面的 Ada 版本创建了一个泛型，因此我可以使用任何浮点类型对其进行实例化。 （异常处理程序是必要的，因为 2^31 - 1 转换为 32 位浮点数并返回溢出......）

with Ada.Text_IO;   
use Ada.Text_IO;

procedure FP_Torture is

    generic
       type Float_Type is digits <>;
    procedure Test_FP;

    procedure Test_FP is
       F : Float_Type;
    begin
       -- for ( n = 2; 1 / n * n == 1; ++ n ) ;
       for i in 2 .. Natural'Last loop
          F := Float_Type(i);
          exit when 1.0 / F * F /= 1.0;
       end loop;
       Put_Line(natural'image(natural(F)) & " (" 
               & Float_Type'image(1.0 - (1.0 / F * F)) & ")");

       -- for ( ; (int) ( 1 / n * n ) == 1; ++ n ) ;
       for i in 1 .. Natural'Last  loop
          F := Float_Type(i);
          exit when natural(1.0 / F * F) /= 1;
       end loop;
       Put_Line(Natural'image(Natural(F)) & " (" 
               & Float_Type'image(1.0 - (1.0 / F * F)) & ")");
    exception
       when Constraint_Error => 
           Put_Line("Error representing float " & Float_Type'image(F) 
                    & " as integer");
    end;

    type Big_Float is digits 18;

    procedure Test7 is new Test_FP(Float);
    procedure Test15 is new Test_FP(Long_Float);
    procedure Test18 is new Test_FP(Big_Float);

begin
    Test7;
    Test15;
    Test18;
end FP_Torture;