2

当我运行以下代码时：

#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

那么为什么它们当时不同而现在相同呢？

4

1 回答 1

4

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

其次，您不能确定您的 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

（与您的示例相同）；

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

4.47920533832942505776
4.47920533832942505820

5.17737751763962026144
5.17737751763962025971

与 Linux 案例相同；因此，在real long double 中，与 100 个加法有一些实际差异，按照常识，它大于 50。但是您的平台默认舍入为 double。

最后，一般来说，这是有限精度和随之而来的舍入的众所周知的效果。例如，在这本经典著作中，这种对递减数列和的错误舍入在第一章中进行了解释。

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

更新：如果是 Windows，您可以使用_controlfp()和_controlfp_s() 操作 FPU 模式。在 Linux 中，_FPU_SETCW 也是如此。此描述详细说明了一些细节并提供了示例代码。

UPDATE2：使用Kahan 求和在所有情况下都能得到稳定的结果。下面显示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要长一倍：

4.47920533832942505776 4.47920533832942524555
4.47920533832942505820 4.47920533832942524555

100和FPU长一倍：

5.17737751763962026144 5.17737751763961995266
5.17737751763962025971 5.17737751763961995266

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

于 2015-12-02T16:27:18.050 回答