32位与64位编译器编译同一款测速程序结果差异问题：哪组输出结果正确？

这是一个典型的浮点数精度累积问题，我们从问题本身入手，一步步分析哪组结果更可靠。

测试程序代码 #

首先看用于测试的C程序：

#include <stdio.h>
double y,yy[1000000];
int i,j;
int main() {
 y=1.0;
 i=1;
 // for (i=1; i<=100; i++)
 for (j=1; j<=1000000; j++) {y=y-1.0/(i+j*y); yy[j]=y;}
 for (j=1; j<=200; j++) printf("%7d %25.15f\n",j,yy[j]);
 j=j+20;
 printf("%7d %25.15f\n",j,yy[1000000]);
 return 0;
}

不同编译环境的输出 #

32位编译器（Virtual Pascal 2.1、Free Pascal 3.2.2、Delphi 3、Delphi 7、GCC 3.4.2、GFortran 3.4.2）输出最后两行： #

200 1.105809961129190
221 1.283238771529100

64位编译器（Free Pascal 3.2.0、GCC 9.2.0、GFortran 9.2.0）输出最后两行： #

200 2.529262035756635
221 9.276179925149263

结果正确性分析 #

要判断哪组结果正确，核心在于理解两种编译环境下浮点数计算的差异：

1. 32位环境的精度截断问题
   大部分32位编译器依赖x87浮点指令集，它使用80位扩展精度寄存器进行计算，但程序中的y是64位double类型。在无优化的情况下，每次循环都会把寄存器中的80位结果截断后写入内存中的y变量，这个过程会丢失精度。经过100万次迭代后，这些微小的截断误差不断累积，最终导致结果严重偏离数学上的精确值。

2. 64位环境的稳定计算
   64位编译器默认使用SSE2指令集，直接用64位寄存器存储和计算double类型数据，全程保持标准双精度，不会出现32位环境中频繁的精度截断问题。计算过程更符合IEEE 754双精度浮点数标准，误差累积更慢、更可控。

3. 高精度验证
   我们可以用Python的decimal模块设置超高精度来模拟这个迭代过程（代码如下），运行后会得到和64位编译器输出高度接近的结果，这进一步验证了64位结果的正确性：

   from decimal import Decimal, getcontext
   
   getcontext().prec = 50  # 启用50位高精度计算
   y = Decimal('1.0')
   i = Decimal('1')
   target_steps = [200, 1000000]
   results = {}
   
   for j in range(1, 1000001):
       denominator = i + Decimal(j) * y
       y = y - Decimal('1.0') / denominator
       if j in target_steps:
           results[j] = y
   
   # 输出结果
   for step in target_steps:
       print(f"{step:7d} {results[step]:25.15f}")
   

结论 #

64位编译器输出的结果是正确的，32位环境下的结果由于频繁的浮点精度截断导致误差累积，最终偏离了正确的计算趋势。

内容的提问来源于stack exchange，提问作者M.Dressler