Approximation to constant “pi” does not get any be

In R I have written this function

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)
}

This function approximates constant pi, however, after m > 50 the approximation gets stuck, i.e. the approximation is the same value and don't get better. How can I fix this? Thanks.

Let's go inside:

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

There is no further improvement after 51, because:

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

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

Further terms are too small compared with the 1st term, readily beyond what machine precision could measure.

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

So eventually you are just adding zeros.

In this case, the summation over o has numerically converged, so does your approximation to pi.

More thoughts

IEEE 754 standard for double-precision floating-point format states that on a 64-bit machine: 11 bits are used for exponential, 53 bits are used for significant digits (including a sign bit). This gives the machine precision: 1 / (2 ^ 52) = 2.2204e-16. In other words, a double-precision floating point number at most has 16 valid significant digits. R function print can display up to 22 digits, while sprintf can display more, but remember, any digits beyond the 16th are invalid, garbage values.

Have a look at the constant pi in R:

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

If you compare it with How to print 1000 decimals places of pi value?, you will see that only the first 16 digits are truly correct:

3.141592653589793

What could alternative be done, so I can calculate more digits using my approach?

No. There have been many crazy algorithms around so that we can compute a shockingly great many of digits of pi, but you cannot modify your approach to get more valid significant digits.

At first I was thinking about computing sum(o[1:51]) and sum(o[52:100]) separately, as both of them give 16 valid significant digits. But we can't just concatenate them to get 32 digits. Because for sum(o[1:51]), the true digits beyond the 16th are not zeros, so the 16 digits for sum(o[52:100]) are not the 17 ~ 32th digits of sum(o[1:100]).