c - 为什么这段代码返回两个不同的值做同样的事情？

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。

我的代码返回: -629.47620126173774000000 -629.47620126173763000000

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

Answer 1

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 中的值，导致结果略有不同。

另一种可能是编译器做了“浮点收缩”。再次引用 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 可能更准确地表示您指定的计算(假设所有输入都限制为可以在 a 中表示的值双).

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