/*
The Gauss AGM algorithm using Schonhage variation.
[1] Jorg Arndt, Christoph Haenel, "Pi Unleashed", pp. 93, Springer, 2001.
(C) 2003 Hironobu SUZUKI, licensed by GPL2 
$Id: pi.c,v 1.1 2003/09/16 01:20:41 hironobu Exp hironobu $
*/
#include <stdio.h>
#include <gmp.h>
#include <sys/timeb.h>
#include <stdlib.h>

#define LOG_TEN_TWO  3.32192809488736234789
#define bprec(n) (int)(((n+10)*LOG_TEN_TWO)+2)

int
main (int ac, char *av[])
{
  long int k, loopmax;
  mpf_t A, B, a, b, s, s_1, t, t1, t2, t3, c2;
  long int prec, dprec;
  struct timeb start_t, end_t;

#ifdef VERBOSE
  printf("RUN WITH GMP VERSION: %s\n",gmp_version);
  ftime(&start_t);
#endif


  dprec = 1000000L;		/* decimal precision */
  prec = bprec(dprec);		/* binary precision (plus alpha) */
  mpf_set_default_prec (prec);
  loopmax = 21;

  mpf_init (A);			/* big A */
  mpf_init (B);			/* big B */
  mpf_init (a);			/* a */
  mpf_init (b);			/* b */
  mpf_init (s);			/* s(n) */
  mpf_init (s_1);		/* s(n-1) */
  mpf_init (t);			/* temporary */
  mpf_init (t1);		/* temporary */
  mpf_init (t2);		/* temporary */
  mpf_init (t3);		/* temporary */
  mpf_init (c2);		/* 2 constant */
  mpf_set_ui (A, 1);
  mpf_set_ui (a, 1);
  mpf_set_ui (t1, 1);
  mpf_div_ui (B, t1, 2);
  mpf_div_ui (s, t1, 2);
  mpf_set_ui (t1, 10);
  mpf_set_ui (c2, 2);
  for (k = 1; k <= loopmax; k++)
    {
      mpf_add (t1, A, B);	/* (A+B) */
      mpf_div_ui (t, t1, 4);	/*  t = (A+B)/4 */
      mpf_sqrt (b, B);		/* b = sqrt(B) */
      mpf_add (t1, a, b);	/* (a+b) */
      mpf_div_ui (a, t1, 2);	/* a = (a+b)/2 */
      mpf_mul (A, a, a);	/* A = a * a */
      mpf_sub (t1, A, t);	/*  (A-t) */
      mpf_mul_ui (B, t1, 2);	/* B = (A - t) * 2 */
      mpf_sub (t1, B, A);	/* (B-A) */
      mpf_pow_ui (t2, c2, k);	/* 2^k */
      mpf_mul (t3, t1, t2);	/* (B-A) * 2^k  */
      mpf_set (s_1, s);
      mpf_add (s, s_1, t3);	/* s = s + (B-A) * 2^k */
    }
  mpf_add (t1, a, b);
  mpf_pow_ui (t2, t1, 2);
  mpf_mul_ui (t1, s, 2);
  mpf_div (A, t2, t1);
  mpf_out_str (stdout, 10, dprec + 10, A);
#ifdef VERBOSE
  ftime(&end_t);
  {
    int sec, milli;
    sec=end_t.time - start_t.time;
    milli=end_t.millitm - start_t.millitm;
    if ( milli < 0 ) {
      sec--;
      milli =  1000 - milli;
    }
    fprintf(stderr,"Time: %d.%d sec\n",sec ,milli);
  }
#endif

  exit(0);
}