/* Copyright (C) 2014 by Alexandru Cojocaru */

/* This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>. */


/* Compile: gcc -std=c11 -lgmp */

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

#include <error.h>

void
salamin_brent (int k, mpf_t pi)
{
  mpf_t a, a_k_1, b, c, s, t;
  
  mpf_inits (a, a_k_1, b, c, s, t, NULL);

  mpf_set_ui (a_k_1, 1);

  mpf_sqrt_ui (b, 2);
  mpf_ui_div (b, 1, b);

  mpf_set_ui (s, 2);
  mpf_ui_div (s, 1, s);

  for (int i = 1; i <= k; ++i) {
    mpf_add (a, a_k_1, b);
    mpf_div_ui (a, a, 2);

    mpf_mul (b, a_k_1, b);
    mpf_sqrt (b, b);

    mpf_mul (c, a, a);
    mpf_mul (t, b, b);
    mpf_sub (c, c, t);

    mpf_set_ui (t, 2);
    mpf_pow_ui (t, t, i);
    mpf_mul (t, t, c);
    mpf_sub (s, s, t);

    mpf_set (a_k_1, a);
  }

  mpf_set (pi, a);
  mpf_pow_ui (pi, a, 2);
  mpf_mul_ui (pi, pi, 2);
  mpf_div (pi, pi, s);
  
}

int
main (void)
{
  mpf_set_default_prec (140);

  mpf_t pi;
  mpf_init (pi);
  
  salamin_brent (5, pi);
  mpf_sub_ui (pi, pi, 3);

  if (mpf_out_str (stdout, 10, 30, pi) == 0)
    error (1, 0, "mpf_out_str failed");
}