在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个有效有效数字。R函数print最多可以显示22位数字，而sprintf可以显示更多，但请记住，超过第16位的任何数字都是无效的，垃圾值。
看看R:中的常数pi
sprintf("%.53f", pi)
#[1] "3.14159265358979311599796346854418516159057617187500000"
如果你将其与"如何打印1000个小数位数的π值？"进行比较？，你会看到只有前16位数字是真正正确的：
3.141592653589793
可以做什么替代方案，这样我就可以使用我的方法计算更多的数字？
否。有很多疯狂的算法，所以我们可以计算出数量惊人的pi的数字，但你不能修改你的方法来获得更有效的有效数字。
起初，我考虑分别计算sum(o[1:51])和sum(o[52:100])，因为它们都给出了16个有效的有效数字。但我们不能把它们连接起来得到32位数字。因为对于sum(o[1:51])，第16位之后的真数位不是零，所以sum(o[52:100])的16位不是sum(o[1:100])的第17～32位。