在典型的浮点中，没有倒数的最小正整数是什么

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

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

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

Answer 1

双精度：

1/41=0x1.8f9c18f9c118fap-6，41*0x1.8f9 c18f9 c18 fap-6=0x1.0000000000000028，四舍五入为1。1/45=0x1.6c16c16c16c17p-6，45*0x1.6c16C16c16c17 p-6=0x1.000000000000002c，四舍五入为1。

然而，

1/49=0x1.4e5e0a72f0539p-6，49*0x1.4ee0a72 f0539P-6=0x0.fffffffffffa4，四舍五入为0x0.fffff ffffff8=0x1.fffffff ffff ff0p-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^五十三+d，则e*f>1和e*f=d*f+2⁵³*f<1-2^-53+2^-52=1+2^{>-53，再次四舍五入为1。}

编辑：上面的推理是错误的，因为连续两次双打之间的步幅相差了两倍。一个没有倒数的二重的例子是0x1.ffffffbffffe。0x1.00000002000001p-1太小，但0x1.000000020000 02p-1太大。没有倒数的整数的最小例子是237。1/237大致为0x1.1485f0e0acd3B68c6Bp-8，四舍五入为0x1.1485 f0e0Acd58p-8。这个数字太小，而后面的下一个双位数太大。

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++的故障点，并证实了38780位浮点的使用。但是将(一个非常接近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;