为什么R在这里错误地执行求和？

(如果答案确实是R如何使用二进制数来存储小小数,那么仍然会喜欢听到有关为什么会发生这种情况的人的详细信息)

我用R来计算总和,特别是这个:

这是我的R代码:

r = 21 #Number of times the loop executes 
n = 52 
sum = 0 #Running total 
for(k in 0:r){
  sum = sum + (1-(choose(2^(k), n)*(factorial(n)/(2^(n*k)))))
  print(sum)
}

如果你看一下输出,你会注意到:

[1] 1
[1] 2
...
[1] 11.71419
[1] 11.71923
[1] 11.72176
[1] 12.72176
[1] 13.72176

为什么在第19次迭代后它开始递增1？

还有其他免费提供的计算引擎更适合这项任务吗？

所有计算都在这里使用浮点数,这通常不适合阶乘和提高功率(因为这些值变得非常快,这使得计算不太准确).

请试试:

> factorial(52)/(2^(52*19))
[1] 3.083278e-230
> factorial(52)/(2^(52*20))
[1] 0

并与Pari-GP进行比较:

? \p 256
realprecision = 269 significant digits (256 digits displayed)
? precision( (52!) / 2^(52*19) + .  , 24)
%1 = 3.08327794368108826958435659488643289724 E-230
? precision( (52!) / 2^(52*20) + .  , 24)
%2 = 6.8462523287873017654100595727727116496 E-246

如果您正在寻找一种方法来解决您遇到的溢出/下溢问题,您可以使用日志(以保持中间计算的合理范围),然后在结束时取幂:

options(digits=20)

for(k in 0:r){
  sum = sum + (1 - exp(lchoose(2^k, n) + log(factorial(n)) - k*n*log(2)))
  print(paste0(k,": ",sum))
}

[1] "0: 1"
[1] "1: 2"
[1] "2: 3"
...
[1] "19: 11.7217600143979"
[1] "20: 11.7230238079842"
[1] "21: 11.7236558993777"

为了检查这是否正确,我在Mathematica中运行了原始求和(不带日志)并得到了相同的结果到12位小数.

虽然你可以解决这个问题R,如果你想使用计算机代数系统(允许你进行符号数学和精确计算),Sage是免费的开源软件.