如果用32位编译器(Virtual Pascal 2.1、Free Pascal 3.2.2、Delphi 3、Delphi 7、GCC 3.4.2、GFortran 3.4.2(或64位编译器(Free Pascal 3.2.0、GCC 9.2.0、GFortran9.2.0(编译，测试速度的非常简单的程序会给出不同的结果。程序代码为：
#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.15fn",j,yy[j]);
j=j+20;
printf("%7d %25.15fn",j,yy[1000000]);
return 0;
}
如果使用32位编译器编译，则最后两行输出为
200 1.105809961129190
221 1.283238771529100
如果使用64位编译器编译，则最后两行输出为
200 2.529262035756635
221 9.276179925149263
编译期间未应用任何优化。哪个结果是正确的？
您的程序尝试在for循环(当j为1000000时(和printf中访问yy[1000000]。数组有1000000个元素，所以最后一个元素的索引只有999999，而不是1000000。超越yy[999999]的访问不是由C标准定义的。这可能会损坏内存或导致其他影响，从而导致您观察到的输出。
要解决此问题，请将yy的定义更改为具有1000001个元素，或者将代码更改为使用索引0到999999，而不是1到1000000。
当我查看您的代码时，我看到了两个问题；
-
对数组的越界访问，这是@Eugene Sh的评论和@Eric Postpischil的回答。我相信这不是你在这里观察到的东西的来源。
-
程序执行的操作序列的性质。那里的操作序列非常敏感，因为小数字的子运算(加法也一样(。但这只是问题的间接根源。真正的问题很可能是，在64位程序中，有更多的寄存器可用于局部变量和中间结果，这使得可以以更高的精度进行累积。当您使用32位程序时，可用于变量和中间结果的寄存器较少，从而迫使累积的中间结果以较低的精度存储回内存(或缓存(。为了更好地理解这一点，有一点很重要，那就是CPU几乎总是使用比变量精度更大的寄存器进行计算，然后在将结果存储回内存时，结果会被截断。
如果你可以从y的最后一个值开始累积(从你的模型来看这似乎不容易(，你会看到更少的差异。同样，如果你做更少的迭代，你会看到更少的差异。
查看反汇编。有了gcc -m32，我在一定程度上看到：
11ef: db 45 f4 fildl -0xc(%ebp)
11f2: dd 83 40 00 00 00 fldl 0x40(%ebx)
11f8: de c9 fmulp %st,%st(1)
11fa: de c1 faddp %st,%st(1)
11fc: d9 e8 fld1
11fe: de f1 fdivp %st,%st(1)
1200: de e9 fsubrp %st,%st(1)
1202: dd 9b 40 00 00 00 fstpl 0x40(%ebx)
这是使用x87 FPU来计算结果。这些寄存器默认为80位浮点数字，精度为64位。
将其与gcc(无-m32选项(进行比较：
118b: f2 0f 59 cb mulsd %xmm3,%xmm1
118f: f2 0f 58 d1 addsd %xmm1,%xmm2
1193: f2 0f 10 0d 7d 0e 00 movsd 0xe7d(%rip),%xmm1 # 2018 <_IO_stdin_used+0x18>
119a: 00
119b: f2 0f 5e ca divsd %xmm2,%xmm1
119f: f2 0f 5c c1 subsd %xmm1,%xmm0
11a3: f2 0f 11 05 b5 2e 00 movsd %xmm0,0x2eb5(%rip) # 4060 <y>
11aa: 00
11ab: 8b 05 d3 40 7a 00 mov 0x7a40d3(%rip),%eax # 7a5284 <j>
11b1: f2 0f 10 05 a7 2e 00 movsd 0x2ea7(%rip),%xmm0 # 4060 <y>
此处使用SSE2寄存器。这些是精度为53位的64位浮点数。
我没有对此进行充分分析，但这可以解释差异。
关于上次打印结果的其他评论仍然正确。
为了看看正确的答案是什么，我使用了以下最大值程序：
fpprec: 200; /* 200 digits of precision */
y: 1b0; /* bigfloat number */
for j:1 to 200 do (y: y - 1/(1+j*y), printf(true, "~d: ~m~%", j, y));
最后两个值是：
1: 5.0b-1
2: 0.0b0
3: - 1.0b0
4: - 6.6666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667b-1
5: - 2.3809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809523809524b-1
...
199: 2.9134730438104486835684450669343781994684160775783702527498016474078580803950756073065445760100932978294411495273196565514143473753599906758526639607462515211042188884927808719790375365190366781832406b0
200: 2.9117598191464059260125139410787531875196505670673442081756913620541464298655107172850569895426676723484404517299155462989843760200675720715542525681413183321106093832508069466366241521724422505183095b0
Maxima可以将精确的值计算为有理数，但分数很快就会有大量的数字，所以速度慢了很多，我不想等待。因此，具有200位数字的bigfoats。更大数量的数字可以很容易地完成；我只是希望200位数就足够了。