/* 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 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); }