为什么此代码 return 两个不同的值做同样的事情？

Question

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

int main()
{
    double a;
    double b;
    double q0 = 0.5 * M_PI + 0.5 * -2.1500000405000002;
    double q1 = 0.5 * M_PI + 0.5 * 0.0000000000000000;
    double w0 = 0.5 * M_PI + 0.5 * -43000.0008100000050000;
    double w1 = 0.5 * M_PI + 0.5 * -0.0000000000000000;
    double m = 1;
    double g = 43000000.81;
    double l1 = 0.1;
    double l2 = 0.1;
    double h = 0.0001;

    a = ((-g / l1) * sin(q0) + (sin(q1 - q0) * (cos(q1 - q0) * (w0 * w0 + (g / l1) * cos(q0)) + l2 * (w1 * w1 / l1))) / (m + pow(sin(q1 - q0), 2)));
    a = h * a;

    b = h * ((-g / l1) * sin(q0) + (sin(q1 - q0) * (cos(q1 - q0) * (w0 * w0 + (g / l1) * cos(q0)) + l2 * (w1 * w1 / l1))) / (m + pow(sin(q1 - q0), 2)));

    printf("%.20lf ", a);
    printf("%.20lf", b);
    return 0;
}

我对 a 和 b 进行了相同的计算，不同之处在于我分两步获得 a 的值，而一步获得 b。

我的代码returns： -629.47620126173774000000 -629.47620126173763000000

最后两位小数不同的原因是什么？

Answer 1

答案是因为在浮点数计算中方程 a=bcde 不必等于 x = bc y = de 和 a.= xy。

这是因为浮点运算的精度有限。

Answer 2

C 标准（99 和 11）说：

The values of operations with floating operands and values subject to the usual arithmetic conversions and of floating constants are evaluated to a format whose range and precision may be greater than required by the type.

因此，在 h*(X+Y) 之类的表达式中，就像您在 b 的赋值中一样，允许实现对 X+Y 的中间结果使用比可能更高的精度存储在 double 中，即使子表达式的类型仍被认为是 double。但是在 a=X+Y; a=h*a; 中，第一个赋值强制该值是实际可以存储在 double 中的值，导致结果略有不同。

另一种可能是编译器做了"floating point contraction"。再次引用 C 标准，

A floating expression may be contracted, that is, evaluated as though it were an atomic operation, thereby omitting rounding errors implied by the source code and the expression evaluation method.

如果处理器只有一条指令可以一步完成浮点加法和乘法运算，并且编译器决定使用它，则很可能会发生这种情况。

假设其中一个或两个是原因，您的值 b 可能更准确地表示您指定的计算（假设所有输入都限制为可以用 double).

关于宏的 cppreference 页面 FLT_EVAL_METHOD 更详细地讨论了这两个问题。找出 FLT_EVAL_METHOD 的值并使用 #pragma STDC FP_CONTRACT OFF.

可能会很有趣

为什么此代码 return 两个不同的值做同样的事情？

Why does this code return two different values doing the same?

c

precision

binary

numerical-computing

arbitrary-precision