/* mpf_base.c  one time computation of constants to "digits" precision */
#include "mpf_math.h"
#include <stdio.h>
#include <stdlib.h>

static int initial = 1;
static mpf_t twopi;
static mpf_t pi;
static mpf_t halfpi;
static mpf_t fourthpi;
static mpf_t e;
static mpf_t *brfact;
static void mpf_base(); /* only once */

void mpf_twopi(mpf_t result)
{
  if(initial==1) mpf_base();
  mpf_set(result,twopi);
}

void mpf_pi(mpf_t result)
{
  if(initial==1) mpf_base();
  mpf_set(result,pi);
}

void mpf_halfpi(mpf_t result)
{
  if(initial==1) mpf_base();
  mpf_set(result,halfpi);
}

void mpf_fourthpi(mpf_t result)
{
  if(initial==1) mpf_base();
  mpf_set(result,fourthpi);
}

void mpf_e(mpf_t result)
{
  if(initial==1) mpf_base();
  mpf_set(result,e);
}

void mpf_rfact(mpf_t rfact[])
{
  int i; 

  if(initial==1) mpf_base();
  for(i=0; i<=digits; i++) mpf_init_set(rfact[i],brfact[i]);
}

static void mpf_base()
{
  int i, n;
  mpf_t mpf_cnt, tmp;
  mpf_t b1, b2, a2, sum, mp1, mp2, mp3;

  mpf_set_default_prec(digits*3.32);
  mpf_init(mpf_cnt);
  mpf_init(tmp);
  mpf_init(twopi);
  mpf_init(pi);
  mpf_init(halfpi);
  mpf_init(fourthpi);
  mpf_init(e);
  mpf_init(mpf_cnt);
  mpf_init(tmp);
  /* other variable initialized as used */

  /* precompute reciprocal factorials for for future use */
  brfact = (mpf_t *)malloc((digits+1)*sizeof(mpf_t));
  for(i=0; i<=digits; i++) mpf_init(brfact[i]);
  mpf_set_si(brfact[0],1);
  mpf_set_si(brfact[1],1);
  for(i=2; i<=digits; i++)
  {
    mpf_set_si(mpf_cnt,i);
    mpf_div(brfact[i],brfact[i-1],mpf_cnt);
  }

  /* precompute pi/4, pi/2, pi, 2pi for future use */
  mpf_init(b1);
  mpf_init(b2);
  mpf_init(a2);
  mpf_init_set_si(sum,0);
  mpf_init_set_si(mp1,1);
  mpf_init_set_si(mp2,2);
  mpf_init_set_si(mp3,3);
  
  mpf_sqrt(b1,mp3); 
  mpf_sub(b1,mp2,b1); /* b1 = 2-sqrt(3) */
  mpf_sqrt(b2,b1);
  mpf_mul(b2,mp2,b2);
  mpf_sub_ui(b2,b2,1);
  mpf_div(b2,b2,b1); /* b2 = (2*sqrt(b1)-1)/b1 */
  /* compute pi/24 = atan(b2) to  digits precision */
  n = ((digits-1)/4)*4+1; /* last term */
  mpf_mul(a2,b2,b2);
  for(i=n; i>1; i=i-4)
  {
    mpf_div_ui(tmp,mp1,i);
    mpf_add(sum,sum,tmp);
    mpf_mul(sum,sum,a2);
    mpf_div_ui(tmp,mp1,i-2);
    mpf_sub(sum,sum,tmp);
    mpf_mul(sum,sum,a2);
  }
  mpf_add_ui(sum,sum,1);
  mpf_mul(sum,sum,b2);
  mpf_mul_ui(fourthpi,sum,6);
  mpf_mul_ui(halfpi,sum,12);
  mpf_mul_ui(pi,sum,24);
  mpf_mul_ui(twopi,pi,2);

  initial = 0;
} /* end mpf_base */
/* end mpf_base.c */