为什么这两个片段会产生不同的价值？

Question

当我运行以下代码时：

#include <stdio.h>

int main()
{
    int i = 0;
    volatile long double sum = 0;
    for (i = 1; i < 50; ++i) /* first snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf\n", sum);
    sum = 0;
    for (i = 49; i > 0; --i) /* second snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf", sum);
    return 0;
}

输出为：

4.47920533832942346919
4.47920533832942524555

这两个数字不应该相同吗？更有趣的是，下面的代码：

#include <stdio.h>

int main()
{
    int i = 0;
    volatile long double sum = 0;
    for (i = 1; i < 100; ++i) /* first snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf\n", sum);
    sum = 0;
    for (i = 99; i > 0; --i) /* second snippet */
    {
        sum += (long double)1 / i;
    }
    printf("%.20Lf", sum);
    return 0;
}

产生：

5.17737751763962084084
5.17737751763962084084

那么为什么他们现在不同而现在相同呢？

Answer 1

首先，请更正您的代码。按照 C 标准，%lf 不是 *printf 的主体（'l' 是无效的，数据类型保持双精度）。要打印 long double，应该使用 %Lf。使用您的变体 %lf，可能会遇到格式不正确、缩减值等错误。（您似乎运行ning 32 位环境：在 64 位中，Unix 和 Windows 在 XMM 寄存器中传递 double，但在其他地方传递 long double - Unix 的堆栈，Windows 的指针内存。在 Windows/x86_64 上，你的代码将出现段错误，因为被调用者需要指针。但是，Visual Studio, long double 是 AFAIK 的 double 别名，所以你可以对这个变化一无所知。）

其次，您不能确定这段代码没有被您的 C 编译器针对编译时计算进行优化（这可以比默认的运行-time 计算更精确）。为避免此类优化，请将 sum 标记为 volatile。

进行这些更改后，您的代码显示：

在 Linux/amd64，gcc4.8:

50 人：

4.47920533832942505776
4.47920533832942505820

100:

5.17737751763962026144
5.17737751763962025971

在 FreeBSD/i386，gcc4.8，没有精度设置或显式 fpsetprec(FP_PD):

4.47920533832942346919
4.47920533832942524555

5.17737751763962084084
5.17737751763962084084

（与您的示例相同）；

但是，在 FreeBSD 上使用 fpsetprec(FP_PE) 进行相同的测试，它将 FPU 切换为真正的 long double 操作：

4.47920533832942505776
4.47920533832942505820

5.17737751763962026144
5.17737751763962025971

与Linux大小写相同；所以，在real long double中，和100的加数是有一定区别的，按照常理来说，比50大。但是你的平台默认四舍五入到double。

最后，一般来说，这是众所周知的有限精度和随之而来的舍入效应。例如，在 this classical book 中，这种递减数列和的误舍入在最开始的章节中得到了解释。

我现在还没有真正准备好调查 50 个加数和四舍五入加倍的结果来源，为什么它显示出如此巨大的差异以及为什么这个差异用 100 个加数来补偿。这需要比我现在负担得起的更深入的调查，但是，我希望这个答案清楚地告诉你下一个要挖掘的地方。

更新：如果它是 Windows，您可以使用 _controlfp() and _controlfp_s(). In Linux, _FPU_SETCW does the same. This description 操作 FPU 模式详述了一些细节并给出了示例代码。

UPDATE2：使用Kahan summation 在所有情况下都能得到稳定的结果。下面显示4个值：升序i，没有KS；升序 i, KS;下降 i，没有 KS；降序 i, KS:

50 和 FPU 加倍：

4.47920533832942346919 4.47920533832942524555
4.47920533832942524555 4.47920533832942524555

100 和 FPU 加倍：

5.17737751763962084084 5.17737751763961995266
5.17737751763962084084 5.17737751763961995266

50 和 FPU 到 long double:

4.47920533832942505776 4.47920533832942524555
4.47920533832942505820 4.47920533832942524555

100 和 FPU 到 long double:

5.17737751763962026144 5.17737751763961995266
5.17737751763962025971 5.17737751763961995266

可以看到差异消失了，结果稳定了。我认为这几乎是可以在此处添加的最后一点:)

为什么这两个片段会产生不同的价值？

Why these two snippets generate different value?

c

floating-point

long-double