c - 为什么这两个片段产生不同的值(value)？

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”为 void，数据类型仍然为 double)。要打印 long double，应该使用 %Lf。使用您的变体 %lf，可能会遇到格式、截断值等不正确的错误。(您似乎运行 32 位环境:在 64 位中，Unix 和 Windows 都通过 double在 XMM 寄存器中，但 long double 在其他地方 - Unix 的堆栈，Windows 的指针内存。在 Windows/x86_64 上，您的代码将出现段错误，因为被调用者需要指针。但是，使用 Visual Studio，据我所知，long double 别名为 double，因此您可以对这一变化一无所知。)

其次，您无法确定此代码没有被您的 C 编译器针对编译时计算进行优化(这可以比默认运行时计算更精确)。为了避免这种优化，请将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 切换为实长 double 运算:

4.47920533832942505776
4.47920533832942505820

5.17737751763962026144
5.17737751763962025971

与Linux情况相同；因此，在真实长 double 中，与 100 个被加数存在一些实际差异，并且根据常识，它比 50 个被加数大。但是您的平台默认四舍五入为双倍。

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

我现在还没有真正准备好调查 50 个被加数并四舍五入到双倍的结果来源，为什么它显示出如此巨大的差异以及为什么这种差异用 100 个被加数来补偿。这需要比我现在能承担的更深入的调查，但是，我希望这个答案清楚地向您展示了下一个需要挖掘的地方。

更新:如果是 Windows，您可以使用 _controlfp() 操作 FPU 模式。和 _controlfp_s().在 Linux 中，_FPU_SETCW 执行相同的操作。 This description详细阐述了一些细节并给出了示例代码。

更新2:使用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 转为长 double :

4.47920533832942505776 4.47920533832942524555
4.47920533832942505820 4.47920533832942524555

100 和 FPU 转为长 double :

5.17737751763962026144 5.17737751763961995266
5.17737751763962025971 5.17737751763961995266

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