// 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 1-5 of http://mathworld.wolfram.com/PiIterations.html
// This has quadratic (i.e. n^2) convergence. Super fast!
//
// 24 iterations will give 34 million digits
// 25 iterations will give 67 million digits
// 26 iterations will give 134 million digits
// 30 iterations will give > 1 billion digits

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

// n bits = n·log10(2) digits of 'pi'
#define PI_DIGITS (1000)
#define PRECISION_BITS (int)(PI_DIGITS * M_LN10 / M_LN2)
#define NUM_ITERATIONS (26)
#define BASE (10)

int main()
{
  mpf_set_default_prec(PRECISION_BITS);

  mpf_t x, y, pi;
  mpf_t tmp[14];

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

  mpf_sqrt_ui(tmp[0], 2);                        // √2
  mpf_sqrt(y, tmp[0]);                           // y = √√2
  mpf_ui_div(tmp[1], 1, y);                      // 1 / √√2
  mpf_add(tmp[2], y, tmp[1]);                    // (√√2 + 1/ √√2)
  mpf_div_ui(x, tmp[2], 2);                      // x = 1/2 (√√2 + 1/√√2)

  // Initial approximation of 'pi'
  mpf_add_ui(pi, tmp[0], 2);                      // pi = 2 + √2

  for (int n = 1; n < NUM_ITERATIONS; n++) {
    // Iteration for 'y'
    mpf_sqrt(tmp[3], x);                         // √x
    mpf_ui_div(tmp[4], 1, tmp[3]);               // 1/√x
    mpf_mul(tmp[5], tmp[3], y);                  // y√x
    mpf_add(tmp[6], tmp[4], tmp[5]);             // y√x + (1/√x)
    mpf_add_ui(tmp[7], y, 1);                    // 1+y
    mpf_div(y, tmp[6], tmp[7]);                  // y -> (y√x + (1/√x)) / (1+y)

    // Iteration for 'pi'
    mpf_add_ui(tmp[8], x, 1);                    // 1+x
    mpf_div(tmp[9], tmp[8], tmp[7]);             // (1+x)/(1+y)
    mpf_mul(tmp[10], pi, tmp[9]);                // pi * (1+x)/(1+y)
    mpf_set(pi, tmp[10]);                        // overwrite pi with above value

    // Iteration for 'x'
    mpf_sqrt(tmp[11], x);                        // √x
    mpf_ui_div(tmp[12], 1, tmp[11]);             // 1/√x
    mpf_add(tmp[13], tmp[11], tmp[12]);          // √x + 1/√x
    mpf_div_ui(x, tmp[13], 2);                   // x -> 1/2 (√x + 1/√x)
  }

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

  mpf_out_str(NULL, BASE, 0, pi); // print 'pi'

  return 0;
}