// Copyright (C) 2016, Satish B. Setty 
// This Source Code Form is subject to the terms of the Mozilla Public License,
// version 2.0, which you can obtain at https://mozilla.org/MPL/2.0/.

// Pi Computation
// Based on Formulas 8-13 of http://mathworld.wolfram.com/PiIterations.html
// This has quartic (i.e. n^4) convergence. Supersonic jet fast!
//
// 10 iterations will give 2 million digits
// 12 iterations will give 32 million digits
// 13 iterations will give 129 million digits
// 15 iterations will give > 2 billion digits

// NOTE: Mathworld's formula (11) has error. See original paper for correct formula (1989)
// "Ramanujan, Modular Equations, and Approximations to Pi, or How to Compute One Billion Digits of Pi"
// URL: http://www.cecm.sfu.ca/~pborwein/PAPERS/P40.pdf

#include <stdio.h>
#include <math.h>
#include <gmp.h>

// n bits = n·log10(2) digits of 'pi'
#define PI_DIGITS (10000000)
#define PRECISION_BITS (int)(PI_DIGITS * M_LN10 / M_LN2)
#define NUM_ITERATIONS (12)
#define BASE (62)

int main()
{
  mpf_set_default_prec(PRECISION_BITS);

  mpf_t x, y, sqrt2, pi;
  mpf_t tmp[13];

  mpf_inits(x, y, sqrt2, pi, NULL);
  for (int i = 0; i < 13; i++) { mpf_init(tmp[i]); }

  mpf_sqrt_ui(sqrt2, 2);                                // sqrt2 = √2
  mpf_sub_ui(y, sqrt2, 1);                              // y0 = √2 - 1
  mpf_mul_ui(tmp[0], sqrt2, 4);                         // 4√2
  mpf_ui_sub(x, 6, tmp[0]);                             // x0 = 6 -  4√2
  mpf_clear(sqrt2);

  // 'n' must start from 1, not 0
  for (int n = 1; n <= NUM_ITERATIONS; n++) {
    // Iteration for 'y'
    mpf_pow_ui(tmp[0], y, 4);                           // y⁴
    mpf_ui_sub(tmp[1], 1, tmp[0]);                      // 1 - y⁴
    mpf_sqrt(tmp[2], tmp[1]);                           // √(1-y⁴)
    mpf_sqrt(tmp[3], tmp[2]);                           // √√(1-y⁴)

    mpf_ui_sub(tmp[4], 1, tmp[3]);                      // numer = 1 - √√(1-y⁴)
    mpf_add_ui(tmp[5], tmp[3], 1);                      // denom = 1 + √√(1-y⁴)
    mpf_div(y, tmp[4], tmp[5]);                         // y -> numer/denom

    // Iteration for 'x'
    //.. first term
    mpf_add_ui(tmp[6], y, 1);                           // 1+y
    mpf_pow_ui(tmp[7], tmp[6], 4);                      // (1+y)⁴
    mpf_mul(tmp[8], x, tmp[7]);                         // x(1+y)⁴
    //.. second term
    mpf_pow_ui(tmp[9], y, 2);                           // y²
    mpf_add(tmp[10], tmp[6], tmp[9]);                   // (1+y)+y²
    mpf_mul(tmp[11], y, tmp[10]);                       // y(1+y+y²)
    mpf_mul_ui(tmp[12], tmp[11], (1 << (2*n+1)));       // (2^(2n+1))·y(1+y+y²)
    //.. (first term) - (second term)
    mpf_sub(x, tmp[8], tmp[12]);
  }

  // Release memory before dumping 'pi' to file
  for (int i = 0; i < 13; i++) { mpf_clear(tmp[i]); }

  mpf_ui_div(pi, 1, x);                                  // pi = 1/x
  mpf_out_str(NULL, BASE, 0, pi);                        // print 'pi' to stdout

  return 0;
}