50 次迭代后，常数“pi”的近似值并没有变得更好

在 R 中我写了这个函数

ifun <- function(m)  {
  o = c() 
  for (k in 1:m) {
    o[k] = prod(1:k) / prod(2 * (1:k) + 1)
    }
  o_sum = 2 * (1 + sum(o))  # Final result

  print(o_sum)
}

该函数近似于常数pi然而，之后m > 50近似值被卡住，即近似值是相同的并且不会变得更好。我怎样才能解决这个问题？谢谢。

我们进去吧：

o <- numeric(100)
for (k in 1:length(o)) {
  o[k] = prod(1:k) / prod(2 * (1:k) + 1)
  }
o
#  [1] 3.333333e-01 1.333333e-01 5.714286e-02 2.539683e-02 1.154401e-02
#  [6] 5.328005e-03 2.486402e-03 1.170072e-03 5.542445e-04 2.639260e-04
# [11] 1.262255e-04 6.058822e-05 2.917211e-05 1.408309e-05 6.814396e-06
# [16] 3.303950e-06 1.604776e-06 7.807016e-07 3.803418e-07 1.855326e-07
# [21] 9.060894e-08 4.429771e-08 2.167760e-08 1.061760e-08 5.204706e-09
# [26] 2.553252e-09 1.253415e-09 6.157124e-10 3.026383e-10 1.488385e-10
# [31] 7.323800e-11 3.605563e-11 1.775874e-11 8.750685e-12 4.313718e-12
# [36] 2.127313e-12 1.049474e-12 5.179224e-13 2.556832e-13 1.262633e-13
# [41] 6.237104e-14 3.081863e-14 1.523220e-14 7.530524e-15 3.723886e-15
# [46] 1.841922e-15 9.112667e-16 4.509361e-16 2.231906e-16 1.104904e-16
# [51] 5.470883e-17 2.709390e-17 1.342034e-17 6.648610e-18 3.294356e-18
# [56] 1.632601e-18 8.092024e-19 4.011431e-19 1.988861e-19 9.862119e-20
# [61] 4.890969e-20 2.425921e-20 1.203410e-20 5.970404e-21 2.962414e-21
# [66] 1.470070e-21 7.295904e-22 3.621325e-22 1.797636e-22 8.924434e-23
# [71] 4.431013e-23 2.200227e-23 1.092630e-23 5.426483e-24 2.695273e-24
# [76] 1.338828e-24 6.650954e-25 3.304296e-25 1.641757e-25 8.157799e-26
# [81] 4.053875e-26 2.014653e-26 1.001295e-26 4.976849e-27 2.473873e-27
# [86] 1.229786e-27 6.113795e-28 3.039627e-28 1.511323e-28 7.514865e-29
# [91] 3.736900e-29 1.858350e-29 9.242063e-30 4.596582e-30 2.286258e-30
# [96] 1.137206e-30 5.656871e-31 2.814078e-31 1.399968e-31 6.965017e-32

print(sum(o[1:49]), digits = 22)
#[1] 0.5707963267948963359544

print(sum(o[1:50]), digits = 22)
#[1] 0.5707963267948964469767

print(sum(o[1:51]), digits = 22)
#[1] 0.570796326794896557999

print(sum(o[1:52]), digits = 22)
#[1] 0.570796326794896557999

之后没有进一步改善51， 因为：

o[51] / o[1]
#[1] 1.641265e-16

o[52] / o[1]
#[1] 8.128169e-17

与第一项相比，其他项太小，很容易超出机器精度可以测量的范围。

.Machine$double.eps
#[1] 2.220446e-16

所以最终你只是添加零。

在这种情况下，总和o已经在数值上收敛，你的近似值也收敛了pi.

更多想法

双精度浮点格式的 IEEE 754 标准规定在 64 位机器上：11 位用于指数，53 位用于有效数字（包括符号位）。这给出了机器精度：1 / (2 ^ 52) = 2.2204e-16。也就是说，一个双精度浮点数最多有16个valid有效数字。 R函数print最多可显示 22 位数字，同时sprintf可以显示更多，但请记住，超过 16 位的任何数字都是无效的，垃圾值.

看看常数pi in R:

sprintf("%.53f", pi)
#[1] "3.14159265358979311599796346854418516159057617187500000"

如果你把它与如何打印 pi 值的小数点后 1000 位？，您会看到只有前 16 位数字是真正正确的：

3.141592653589793

可以采取什么替代措施，以便我可以使用我的方法计算更多数字？

No. 周围有许多疯狂的算法，因此我们可以计算出数量惊人的数字pi，但您无法修改您的方法以获得更有效的有效数字。

一开始我想到的是计算机sum(o[1:51]) and sum(o[52:100])分别，因为它们都给出 16 个有效有效数字。但我们不能只是将它们连接起来就得到 32 位数字。因为对于sum(o[1:51])，第 16 位之后的真实数字不是零，因此 16 位数字sum(o[52:100])不是第17~32位数字sum(o[1:100]).