是否存在两个相互相乘（或相除）的数字会引入错误？

Question

这是我正在进行的测试的bank，了解FP 基本操作（+、-、*、/）如何引入错误：

#include <iostream>
#include <math.h>

int main() {
    std::cout.precision(100);

    double a = 0.499999999999999944488848768742172978818416595458984375;

    double original = 47.9;
    double target = original * a;    
    double back = target / a;

    std::cout <<  original << std::endl;
    std::cout <<  back << std::endl;
    std::cout <<  fabs(original - back) << std::endl; // its always 0.0 for the test I did
}

你能告诉我两个值（original 和 a），一次是 *（或 /），由于 FP 数学，引入错误？

如果它们存在，是否可以确定该错误是由* 还是/ 引入的？如何？ （因为你需要两者来恢复价值；80 bit？）

使用+ 很容易（只需将0.499999999999999944488848768742172978818416595458984375 添加到0.5，您就会得到1.0，至于0.5 + 0.5）。

但我无法对 * 或 / 做同样的事情。

Answer 1

输出：

#include <cstdio>

int main(void)
{
    double a = 1000000000000.;
    double b = 1000000000000.;
    std::printf("a = %.99g.\n", a);
    std::printf("a = %.99g.\n", b);
    std::printf("a*b = %.99g.\n", a*b);
}

是：

a = 1000000000000。 a = 1000000000000。 a*b = 999999999999999983222784。

假设 IEEE-754 基本 64 位二进制浮点，正确舍入到最接近，与偶数相关。

显然，999999999999999983222784 与 1000000000000•1000000000000、1000000000000000000000000 的精确数学结果不同。

Answer 2

将任意两个大^†数相乘，可能会出现错误，因为可表示的值在高值范围内的距离很大。

虽然这个误差在绝对值上可能很大，但相对于数字本身的大小来说仍然很小，所以如果你进行反向除法，第一次运算的误差会以相同的比例缩小，并且完全消失。因此，这一系列操作是稳定的。

如果乘法的结果大于可表示的最大值，那么它将溢出到无穷大（可能取决于配置），在这种情况下，反向除法不会产生原始值，而是保持无穷大。

同样，如果你用一个很大的数除，你可能会下溢最小的可表示值，导致零或低于正常值。

^† 数字不一定要很大。在考虑巨大的价值时，更容易理解这个问题。该问题也适用于非常小的值。例如：

2.100000000000000088817841970012523233890533447265625 ×
2.100000000000000088817841970012523233890533447265625

正确结果：

4.410000000000000373034936274052605470949292688633679117285...

浮点结果示例：

4.410000000000000142108547152020037174224853515625

错误：

2.30926389122032568296724439173008679117285652827862296732064351090230047702789306640625
× 10^-16

Answer 3

是否存在两个相乘（或相除）的数字会产生误差？

使用"%a" 更容易看到这一点。

当结果的精度不足时，会发生舍入。通常double 具有 53 位 二进制 精度。将下面的 2 个 27 位数字相乘得到精确的 53 位答案，但 2 个 28 位的数字不能形成 55 位有效答案。

除法很容易演示，试试1.0/n*n。

int main(void) {
  double a = 1 + 1.0/pow(2,26);
  printf("%.15a,  %.17e\n", a, a);
  printf("%.15a,  %.17e\n", a*a, a*a);
  double b = 1 + 1.0/pow(2,27);
  printf("%.15a,  %.17e\n", b, b);
  printf("%.15a,  %.17e\n", b*b, b*b);

  for (int n = 47; n < 52; n += 2) {
    volatile double frac = 1.0/n;
    printf("%.15a,  %.17e %d\n", frac, frac, n);
    printf("%.15a,  %.17e\n", frac*n, frac*n);
  }
  return 0;
}

输出

//v-------v         27 significant bits.
0x1.000000400000000p+0,  1.00000001490116119e+00
//v-------------v   53 significant bits.
0x1.000000800000100p+0,  1.00000002980232261e+00
//v-------v         28 significant bits.
0x1.000000200000000p+0,  1.00000000745058060e+00
//v--------------v  not 55 significant bits.
0x1.000000400000000p+0,  1.00000001490116119e+00
//              ^^^ all zeros here, not the expected mathematical answer.

0x1.5c9882b93105700p-6,  2.12765957446808505e-02 47
0x1.000000000000000p+0,  1.00000000000000000e+00
0x1.4e5e0a72f053900p-6,  2.04081632653061208e-02 49
0x1.fffffffffffff00p-1,  9.99999999999999889e-01  <==== Not 1.0
0x1.414141414141400p-6,  1.96078431372549017e-02 51
0x1.000000000000000p+0,  1.00000000000000000e+00