/* hex_pi.c
   pi = sum k=0..inf  1/16^k (4/(8k+1)-2/(8k+4)-1/(8k+5)-1/(8k+6)
*/
#include <stdio.h>
#include <gmp.h>

int main(int argc, char *argv[])
{
  int k;
  mpf_t r16k, t1, t2, t3, t4, sum, term, tmp;
  mpf_t mp1, mp2, mp4, mp5, mp6, mp8, mp16;
  mpf_t mp8k;
  int digits=100;

  printf("hex_pi.c digits=%d \n", digits);
  mpf_set_default_prec(digits*3.32);

  mpf_init(r16k);
  mpf_init(t1);
  mpf_init(t2);
  mpf_init(t3);
  mpf_init(t4);
  mpf_init(term);
  mpf_init(tmp);
  mpf_init_set_si(sum,0);
  mpf_init_set_si(mp1,1);
  mpf_init_set_si(mp2,2);
  mpf_init_set_si(mp4,4);
  mpf_init_set_si(mp5,5);
  mpf_init_set_si(mp6,6);
  mpf_init_set_si(mp8,8);
  mpf_init_set_si(mp16,16);
  mpf_init(mp8k);

  printf("pi= 3.1415926535897932384626433832795028841971 \n");

  mpf_set(r16k,mp1);

  for(k=0; k<30; k++)
  {
    mpf_mul_ui(mp8k,mp8,k);
    mpf_add_ui(tmp,mp8k,1);
    mpf_div(t1,mp4,tmp); /* t1 = 4/(8k+1) */
    mpf_add_ui(tmp,mp8k,4);
    mpf_div(t2,mp2,tmp); /* t2 = 2/(8k+4) */
    mpf_add_ui(tmp,mp8k,5);
    mpf_div(t3,mp1,tmp); /* t3 = 1/(8k+5) */
    mpf_add_ui(tmp,mp8k,6);
    mpf_div(t4,mp1,tmp); /* t4 = 1/(8k+6) */

    mpf_sub(term,t1,t2);
    mpf_sub(term,term,t3);
    mpf_sub(term,term,t4);
    mpf_mul(term,term,r16k); /* 1/16^k (t1-t2-t3-t4) */
    mpf_add(sum,sum,term);
    gmp_printf("pi=%75.72Ff \n", sum);

    mpf_div(r16k,r16k,mp16);
  }
  printf("pi= 3.1415926535897932384626433832795028841971 \n");

  return 0;
} /* end hex_pi.c */