How can I get more accuracy in calculating PI?

Question

I'm supposed to calculate pi using numbers between 1 through 30 with this formula below

Picture

Here is what I've done so far:

int k, outside = 0;
double inside = 0;

cout << "Please enter in a value for k..." << endl;
cin >> k;

while (k < 1 || k > 30)
{

    cout << "Invalid --- k must be >=1 and <=30" << endl;
    cout << "Please enter in a value for k..." << endl;
    cin >> k;

}

outside = pow(2, k);
    
for (int loopNum = 1; loopNum < k; loopNum++)
{
        
    inside = sqrt(2 + inside);

}

inside = sqrt(2 - inside);
double ApproxPI = outside * inside;

cout << fixed << setprecision(20);
cout << "Approximation of PI = " << ApproxPI << endl;

I feel like I did the formula correctly, but the results aren't accurate. Also, when 'k' is set to 29 or 30, it only gives out the value 0. Does anyone know what I might be doing incorrectly?

Answer 1

Each sqrt() computation incurs 0.5 ULP error. Each addition can incur much more. @Yksisarvinen
Yet it is the severe cancellation of precision with 2 - inside that is the killer.

Although the formula mathematical converges, the raw code implementation does not due to accumulated errors.

inside approaches 2.0 and 2 - inside eventually that difference becomes 0.0. As the difference becomes smaller and smaller, the precision in the difference becomes less. Thus the best pi approximation occurs when inside has about half its most significant binary digits all set to 1 so the value is just less than 2.0, and the remaining digits form the basis for pi. I'd expect the best approximation to be within about 1/2 the usual 16 decimal digits of double.

Does anyone know what I might be doing incorrectly?

To get a better answer, code needs extended math or a different formula.

---

C example demo code

long long ULP_diff(double x, double ref) {
  union {
      double d;
      long long ll;
  } ux, uref;
  ux.d = x;
  uref.d = ref;
  return ux.ll - uref.ll;
}

double pip(int k) {
  double outside = pow(2, k);

  double inside = 0;
  for (int loopNum = 1; loopNum < k; loopNum++) {
    inside = sqrt(2 + inside);
  }
  inside = sqrt(2 - inside);

  double pi = 3.1415926535897932384626433832795;
  double ApproxPI = outside * inside;
  printf("%2d %a %a %15lld %.16f\n", k, pi, ApproxPI, ULP_diff(ApproxPI, pi), ApproxPI);
  return ApproxPI;
}

int main(void) {
  for (int i = 5; i < 32; i++)
    pip(i);
}

Output: Smallest ULP distance at iteration 14.

 5 0x1.921fb54442d18p+1 0x1.91f65f10dd7f5p+1  -2840637887779 3.1403311569547392
 6 0x1.921fb54442d18p+1 0x1.92155f7a3665ap+1   -710223644350 3.1412772509327569
 7 0x1.921fb54442d18p+1 0x1.921d1fcdec626p+1   -177559922418 3.1415138011441455
 8 0x1.921fb54442d18p+1 0x1.921f0fe670767p+1    -44390229425 3.1415729403678827
 9 0x1.921fb54442d18p+1 0x1.921f8beccbd5ep+1    -11097567162 3.1415877252799609
10 0x1.921fb54442d18p+1 0x1.921faaee60d6fp+1     -2774409129 3.1415914215046352
11 0x1.921fb54442d18p+1 0x1.921fb2aee1ba8p+1      -693506416 3.1415923456110768
12 0x1.921fb54442d18p+1 0x1.921fb49eceff7p+1      -173489441 3.1415925765450043
13 0x1.921fb54442d18p+1 0x1.921fb5190a222p+1       -45320950 3.1415926334632482
14 0x1.921fb54442d18p+1 0x1.921fb546e04f1p+1         2742233 3.1415926548075892
15 0x1.921fb54442d18p+1 0x1.921fb532811eap+1       -18619182 3.1415926453212153
16 0x1.921fb54442d18p+1 0x1.921fb4e1045cep+1      -104064842 3.1415926073757197
17 0x1.921fb54442d18p+1 0x1.921fb76cea693p+1       579500411 3.1415929109396727
18 0x1.921fb54442d18p+1 0x1.921fc19c82711p+1      3313760761 3.1415941251951911
19 0x1.921fb54442d18p+1 0x1.921fd5fbb1baap+1      8782278290 3.1415965537048196
20 0x1.921fb54442d18p+1 0x1.921fd5fbb1baap+1      8782278290 3.1415965537048196
21 0x1.921fb54442d18p+1 0x1.922261df7a054p+1    183772607292 3.1416742650217575
22 0x1.921fb54442d18p+1 0x1.9227799aa8796p+1    533740280446 3.1418296818892015
23 0x1.921fb54442d18p+1 0x1.923bd7e25164dp+1   1933437888821 3.1424512724941338
24 0x1.921fb54442d18p+1 0x1.923bd7e25164dp+1   1933437888821 3.1424512724941338
25 0x1.921fb54442d18p+1 0x1.94c583ada5b53p+1  46578493959739 3.1622776601683795
26 0x1.921fb54442d18p+1 0x1.94c583ada5b53p+1  46578493959739 3.1622776601683795
27 0x1.921fb54442d18p+1 0x1.bb67ae8584caap+1 726225619525522 3.4641016151377544
28 0x1.921fb54442d18p+1 0x1p+2 1932961502712552 4.0000000000000000
29 0x1.921fb54442d18p+1 0x0p+0 -4614256656552045848 0.0000000000000000
30 0x1.921fb54442d18p+1 0x0p+0 -4614256656552045848 0.0000000000000000
31 0x1.921fb54442d18p+1 0x0p+0 -4614256656552045848 0.0000000000000000