/*
  Function: void mpf_pi (mpf_t rop)
    Set the value of pi.

  Reference:
    https://www.kurims.kyoto-u.ac.jp/~ooura/pi_fft-j.html

  ---- a formula based on the AGM (Arithmetic-Geometric Mean) ----
    c = sqrt (0.125);
    a = 1 + 3 * c;
    b = sqrt (a);
    e = b - 0.625;
    b = 2 * b;
    c = e - c;
    a = a + e;
    npow = 4;
    do {
        npow = 2 * npow;
        e = (a + b) / 2;
        b = sqrt (a * b);
        e = e - b;
        b = 2 * b;
        c = c - e;
        a = e + b;
    } while (e > SQRT_SQRT_EPSILON);
    e = e * e / 4;
    a = a + b;
    pi = (a * a - e - e / 2) / (a * c - e) / npow;
  ---- modification ----
    This is a modified version of Gauss-Legendre formula
    (by T.Ooura). It is faster than original version.
*/

#include <gmp.h>

void mpf_pi (mpf_t rop);

void mpf_pi (mpf_t a) {
  unsigned long int prec;
  mpf_t SQRT_SQRT_EPSILON, c, b, e;
  unsigned long int log2_npow;

  prec = mpf_get_prec (a);
  mpf_init2 (SQRT_SQRT_EPSILON, prec);
  mpf_init2 (c, prec);
  mpf_init2 (b, prec);
  mpf_init2 (e, prec);

  mpf_set_ui (SQRT_SQRT_EPSILON, 1);
  mpf_div_2exp (SQRT_SQRT_EPSILON, SQRT_SQRT_EPSILON, prec >> 2);

  mpf_set_d (c, 0.125);     /* c = sqrt (0.125); */
  mpf_sqrt (c, c);
  mpf_mul_ui (a, c, 3);     /* a = 1 + 3 * c; */
  mpf_add_ui (a, a, 1);
  mpf_sqrt (b, a);          /* b = sqrt (a); */
  mpf_set_d (e, 0.625);     /* e = b - 0.625; */
  mpf_sub (e, b, e);
  mpf_add (b, b, b);        /* b = 2 * b; */
  mpf_sub (c, e, c);        /* c = e - c; */
  mpf_add (a, a, e);        /* a = a + e; */
  log2_npow = 2;            /* npow = 4; */
  do {
    log2_npow++;            /* npow = 2 * npow; */
    mpf_add (e, a, b);      /* e = (a + b) / 2; */
    mpf_div_2exp (e, e, 1);
    mpf_mul (b, a, b);      /* b = sqrt (a * b); */
    mpf_sqrt (b, b);
    mpf_sub (e, e, b);      /* e = e - b; */
    mpf_add (b, b, b);      /* b = 2 * b; */
    mpf_sub (c, c, e);      /* c = c - e; */
    mpf_add (a, e, b);      /* a = e + b; */
  } while (mpf_cmp (e, SQRT_SQRT_EPSILON) > 0);  /* e > SQRT_SQRT_EPSILON */
  mpf_mul (e, e, e);        /* e = e * e / 4; */
  mpf_div_2exp (e, e, 2);
  mpf_add (a, a, b);        /* a = a + b; */
  mpf_mul (c, c, a);        /* pi = (a * a - e - e / 2) / (a * c - e) / npow; */
  mpf_sub (c, c, e);
  mpf_mul (a, a, a);
  mpf_sub (a, a, e);
  mpf_div_2exp (e, e, 1);
  mpf_sub (a, a, e);
  mpf_div (a, a, c);
  mpf_div_2exp (a, a, log2_npow);

  mpf_clear (e);
  mpf_clear (b);
  mpf_clear (c);
  mpf_clear (SQRT_SQRT_EPSILON);
}