why this implementation of a Gauss Legendre algorithm in c++ is not producing result?

Question

I am learning c++ and I have found Gauss-Legendre Algorithm on wikipedia to approximate pi(link: https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm). I tried to implement it with c++ but it is not producing any result. Here is the code:

#include <iostream.h>
#include <conio.h>
#include <math.h>
#include <iomanip.h>
int main()
    {
        clrscr();
        long double a0 = 1,
                b0 = 1 / sqrt(2),
                t0 = 1 / 4,
                p0 = 1,
                an, bn, pn, tn;
    int i = 0;
    while(i < 3)
    {
        an = (a0 + b0) / 2;
        bn = sqrt(a0 * b0);
        tn = t0 - (p0 * (pow((a0 - an), 2)));
        pn = 2 * p0;
        a0 = an;
        b0 = bn;
        t0 = tn;
        p0 = pn;
    }
    long double pi = (pow((an + bn), 2)) / (4 * tn);
    cout<<pi;
    getch();
    return 0;
  }

When I searched for help I found this but it seem to me a different algorithm - gauss-legendre in c++

UPDATE: After adding i increment program gives wrong result.

Answer 1

Along with the fact that your while loop variable does not update you also declared your doubles incorrectly. It is good practice when declaring doubles, such as "1/4", to write them as "1.0/4.0" or "1/4.0" if you are being lazy.

The reason is that C/C++ will execute the "/" operator on the integers and perform the typecast after the fact. Essentially your t0 = 0 (you can check yourself).

Here is your code with a few modifications to print the full double precision at each while loop iteration.

#include <iostream>
#include <cmath>
#include <limits>

int main()
{
  long double a0=1.0, b0 = 1/sqrt(2),
   t0 = 1.0/4.0, p0 = 1.0;
  long double an,bn,pn,tn;
  int i = 0;
  long double pi;
  typedef std::numeric_limits<double> dbl;
  std::cout.precision(dbl::max_digits10);
  while(i < 4)
    {
      an = (a0 + b0)/2.0;      
      bn = sqrt(a0 * b0);
      tn = t0 - (p0 * (a0-an)*(a0-an));
      pn = 2*p0;
      a0 = an,b0 = bn,p0 = pn,t0 = tn;

      pi = (an+bn)*(an+bn) / (4*tn);
      std::cout << pi << std::endl;      
      i++;
    }
  return 0;
}