C中的浮点算术：epsilon比较

Question

我正在尝试使用 epsilon 比较具有双精度的值。但是，我有一个问题 - 最初我认为差异应该等于 epsilon，但事实并非如此。此外，当我尝试使用连续乘法检查二进制表示时，发生了一些奇怪的事情，我感到很困惑，因此我希望您能对问题和我的思维方式作出解释

#include <stdio.h>

#define EPSILON 1e-10

void double_equal(double a, double b) {
    printf("a: %.12f, b: %.12f, a - b = %.12f\n", a, b, a - b);
    printf("a: %.12f, b: %.12f, b - a = %.12f\n", a, b, b - a);
    if (a - b < EPSILON) printf("a - b < EPSILON\n");
    if (a - b == EPSILON) printf("a - b == EPSILON\n");
    if (a - b <= EPSILON) printf("a - b <= EPSILON\n");
    if (b - a <= EPSILON) printf("b - a <= EPSILON\n");
}

int main(void) {
    double wit1 = 1.0000000001;
    double wit2 = 1.0;
    double_equal(wit1, wit2);
    return 0;
}

输出是：

a: 1.000000000100, b: 1.000000000000, a - b = 0.000000000100
a: 1.000000000100, b: 1.000000000000, b - a = -0.000000000100
b - a <= EPSILON

如果我们在数字 (#define EPSILON 1e-10F) 之后不提供“F”/“f”符号，则 C 中的数字常量被声明为双精度数，因此我在这里看不到 @987654321 中的转换问题@。因此，我为这些特定示例创建了非常简单的程序（我知道它应该包括处理将整数部分转换为二进制数）。

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

char* convert(double a) {
    char* res = malloc(200);
    int count = 0;
    double integral;
    a = modf(a, &integral);

    if (integral == 1) {
        res[count++] = integral + '0';
        res[count++] = '.';
    } else {
        res[count++] = '0';
        res[count++] = '.';
    }

    while(a != 0 && count < 200) {
        printf("%.100f\n", a);
        a *= 2;
        a = modf(a, &integral);
        if (integral == 1) res[count++] = integral + '0';
        else res[count++] = '0';
    }

    res[count] = '\0';
    return res;
}

int main(void) {
    double wit1 = 1.0000000001;
    double diff = 0.0000000001;
    char* res = convert(wit1);
    char* di = convert(diff);
    printf("this: %s\n", res);
    printf("diff: %s\n", di);
    return 0;
}

直接输出：

this: 1.0000000000000000000000000000000001101101111100111
diff: 0.00000000000000000000000000000000011011011111001101111111011001110101111011110110111011

第一个问题：为什么差中有这么多结尾的零？ 为什么二进制点后的结果不同？

但是，如果我们看一下计算过程和小数部分，打印出来（我只展示前几行）：

1.0000000001:
0.0000000001000000082740370999090373516082763671875000000000000000000000000000000000000000000000000000
0.0000000002000000165480741998180747032165527343750000000000000000000000000000000000000000000000000000
0.0000000004000000330961483996361494064331054687500000000000000000000000000000000000000000000000000000

0.0000000001:
0.0000000001000000000000000036432197315497741579165547065599639608990401029586791992187500000000000000
0.0000000002000000000000000072864394630995483158331094131199279217980802059173583984375000000000000000
0.0000000004000000000000000145728789261990966316662188262398558435961604118347167968750000000000000000

第二个问题：为什么会有这么多奇怪的结尾数字？这是浮点运算无法精确表示十进制值的结果吗？

分析减法，我可以看到，为什么结果比 epsilon 大。我遵循程序：

1. 为要减去的序列准备一个由 0 组成的补序列
2. “添加”序列
3. 把开头的那个减去，加到最右边

因此：

   1.0000000000000000000000000000000001101101111100111
 - 1.0000000000000000000000000000000000000000000000000
               |
              \/
   1.0000000000000000000000000000000001101101111100111
"+"0.1111111111111111111111111111111111111111111111111    
 --------------------------------------------------------
  10.0000000000000000000000000000000001101101111100110    
          |
          \/
   0.0000000000000000000000000000000001101101111100111

与epsilon的计算值比较：

0.000000000000000000000000000000000110110111110011 0 1111111011001110101111011110110111011
0.000000000000000000000000000000000110110111110011 1

空格表示区别。

第三个问题：如果我无法比较等于 epsilon 的值，我是否需要担心？我认为这种情况表明了与 epsilon 的公差间隔是为了什么而制定的。但是，有什么我应该改变的吗？

Answer 1

此答案假定您的 C 实现使用 IEEE-754 binary64，也称为其 double 类型的“双”格式。这很常见。

如果 C 实现正确舍入，则 double wit1 = 1.0000000001; 将 wit1 初始化为 1.0000000001000000082740370999090373516082763671875。这款宝选择后者是因为它更接近。

如果正确舍入，用于EPSILON 的1e-10 将产生 0.00000000010000000000000000364321973154977415791655470655996396089904010295867910092180700000000000000000

显然wit1 - 1 超过EPSILON，因此double_equal 中的测试a - b < EPSILON 评估为假。

第一个问题：为什么差中有这么多结尾的零？

计算从第一个 1 到最后一个 1 的位数。在每个数字中，有 53 位。那是因为 double 的有效位有 53 位。您的数字恰好以 1 位结尾，这有点巧合。大约一半的时间，尾随位为 0，四分之一的时间，最后两位为零，依此类推。但是，由于double 的有效位有 53 位，因此从第一个 1 位到作为表示值的一部分的最后一位，将正好有 53 位。

由于您的第一个数字在整数位置以 1 开头，因此它之后最多有 52 位。此时，数字必须四舍五入到最接近的可表示值。

由于您的第二个数字介于 2^-34 和 2^-33 之间，因此它的第一个 1 位位于 2^-34 位置，并且它可以在必须四舍五入之前到达 2⁻⁸⁶ 位置。

第三个问题：如果我不能比较等于 epsilon 的值，我是否需要担心？

为什么要与 epsilon 进行比较？ There is no general solution for comparing floating-point numbers that contain errors from previous operations. 是否可以或应该使用“epsilon 比较”取决于应用程序以及所涉及的操作和数字。

Answer 2

为什么二进制点后的结果不同？

因为这就是区别。

期待别的东西来自认为1.0000000001 和0.0000000001 因为double 具有这两个值。他们不。它们的差异不是 1.0。它们的值接近这两个，每个都有大约 53 个二进制数字。它们的差异接近于unit in the last place 的1.0000000001。

为什么会有这么多奇怪的结尾数字？这是浮点运算无法精确表示十进制值的结果吗？

有点。
double 可以编码大约 2⁶⁴ 个不同的数字。 1.0000000001 和 0.0000000001 不在该集合中。取而代之的是附近的那些看起来像奇怪的结尾数字。

如果我无法比较等于 epsilon 的值，我是否需要担心？我认为这种情况表明了与 epsilon 的公差间隔是为了什么而制定的。但是，有什么我应该改变的吗？

是的，更改 epsilon 的用法。 epsilon 对相对差异有用，而不是绝对差异。非常大的连续double 值相距远大于epsilon。大约 45% 的 double，（小的）在数量级上都小于 epsilon。 if (a - b <= EPSILON) printf("a - b <= EPSILON\n"); 或 if (b - a <= EPSILON) printf("b - a <= EPSILON\n"); 将适用于小型 a, b，即使它们的大小相差数万亿倍。

过于简单化：

if (fabs(a-b) < EPSILON*fabs(a + b)) {
  return values_a_b_are_near_each_other;
}