3

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

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

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

#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";
    }
}

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

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似乎暗示对应于 64 位尾数的 epsilon，即使用 Intel 80 位扩展精度的内部计算。

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

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

4

2 回答 2

4

双精度：

1/41 = 0x1.8F9C18F9C18FAP-6和41 * 0x1.8F9C18F9C18FAP-6 = 0x1.000000000000028，其舍入到1.1 / 45 = 0x1.6c16c16c16c16p-6，45 * 0x1.6c16c16c16c17p-6 = 0x1.0000000000000002c , 四舍五入为 1。

然而，

1/49 = 0x1.4e5e0a72f0539p-6，49*0x1.4e5e0a72f0539p-6 = 0x0.fffffffffffffa4，四舍五入为 0x0.ffffffffffff8 = 0x1.ffffffffffffff0p-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 < f*d < f*(1/f+2 -54 ) <= 1 + 2 -54 * f < 1 + 2 -53，所以 f*d 舍入为 1。如果 d < 1/f，则 f *d 可能会舍入到 1 - 2 -53。如果是，则 f*d 位于 [1 - 2 -53 , 1 - 2 -54 )。如果你取 e = 2 -53 + d，那么 e*f > 1 并且 e*f = d*f + 2 -53 *f < 1 - 2 -53+ 2 -52 = 1 + 2 -53，再次四舍五入为 1。

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

于 2012-11-26T07:41:53.497 回答
0

这个问题似乎与 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;
于 2012-11-26T16:24:36.847 回答