为什么 R 在这里不正确地执行求和？

Question

（如果答案确实是 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？

是否有其他更适合此任务的免费计算引擎？

Answer 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

Answer 2

如果您正在寻找一种方法来解决您遇到的 overflow/underflow 问题，您可以使用日志（以保持中间计算的合理幅度），然后在最后取幂：

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 是免费和开放的来源。

为什么 R 在这里不正确地执行求和？

Why does R incorrectly perform sum here?

precision

r

sum