/* mpf_exp.c  normalize and use series                      */
/* e^x = e^j * e^xc   integer j, xc = x-j                   */
/* xc<1/4 series(xc), xc<1/2 series(xc/2)^2, series(xc/4)^4 */
/* 1 + x + x^2/2! + x^3/3! + ...  use rfact = 1/factorial   */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpf_exp.h" /* includes gmp and digits*/

static int n = 50;  /* need  0.125^n/n! <  10^(-digits) */
static mpf_t zero;
static mpf_t one;
static mpf_t half;
static mpf_t forth;
static mpf_t e;
static mpf_t temp;
static mpf_t sum;
static mpf_t prod;
static mpf_t xp;
static mpf_t rfact[51]; /* reciprocal factorial */
static init = 0;
static int debug = 0;

#define e1check  2.7182818284590452353602874713526624977572

void series_exp(mpf_t sum, mpf_t xc) /* |xc|<=0.125 */
{
  int i;

  mpf_set(xp, xc);
  mpf_add(sum, one, xc); /* 1 + xc/1! */
  for(i=2; i<n; i++)
  {
    mpf_mul(xp, xp, xc); /* first term is xc^2/2! */
    mpf_mul(temp, xp, rfact[i]);
    mpf_add(sum, sum, temp);
  }
} /* end series_exp */

void mpf_exp(mpf_t result, mpf_t x)
{
  mpf_t xc;   /* reduced x */
  mpf_t ep;   /* power of e, e^j */
  mpf_t ej;   /* mpf j */
  int i, j, k;

  if(init==0) /* one time only */
  {
    mpf_set_default_prec(digits*3.32); /* in bits */
    mpf_init(zero);  /* must init every variable */
    mpf_init(one);
    mpf_init(half);
    mpf_init(forth);
    mpf_init(e);
    mpf_init(temp);
    mpf_init(sum);
    mpf_init(prod);
    mpf_init(xp);
    for(i=0; i<=n; i++)  mpf_init(rfact[i]);
    mpf_set_si(zero, 0);
    mpf_set_str(one, "1.0", 10);
    mpf_set_str(half, "0.5", 10);
    mpf_mul(forth, half, half);
    mpf_set(rfact[0], one);
    mpf_set(rfact[1], one);
    mpf_set(sum, one);
    mpf_set(e, one);
    mpf_add(e, e, forth); /* e = (e^(1/4))^4 */
    mpf_set(prod, forth);
    for(i=2; i<=n; i++)
    {
      mpf_add(sum, sum, one);
      mpf_div(rfact[i], rfact[i-1], sum);
      mpf_mul(prod, prod, forth);
      mpf_mul(temp, prod, rfact[i]);
      mpf_add(e, e, temp);
    }
    init = 1;
    mpf_mul(e, e, e); /* (e^(1/4))^2 */
    mpf_mul(e, e, e); /* (e^(1/2))^2 */
    if(debug) gmp_printf("e=%50.47Ff \n", e);
    if(debug) gmp_printf("rfact[%d]=%Fg \n", n, rfact[n]);
    if(debug) gmp_printf("0.125^%d=%Fg \n", n, prod);
    if(debug) gmp_printf("0.125^%d/factorial(%d)=%Fg \n", n, n, temp);
  }
  else
  {
    mpf_set_default_prec(digits*3.32); 
  }
  mpf_init(xc);
  mpf_init(ep);
  mpf_init(ej);

  if(mpf_cmp(x, zero)==0)
  {
    mpf_set(result, one);
    return;
  }
  mpf_abs(temp, x);
  if(mpf_cmp(temp, forth)<=0)
  {
    series_exp(result, x);
    return;
  } 
  if(mpf_cmp(temp, half)<=0)
  {
    mpf_mul(xc, x, half);
    series_exp(result, xc);
    mpf_mul(result, result, result);
    return;
  } 
  if(mpf_cmp(temp, one)<=0)
  {
    mpf_mul(xc, x, forth);
    series_exp(result, xc);
    mpf_mul(result, result, result);
    mpf_mul(result, result, result);
    return;
  } 
  /* get j = (int)(x - floor(x))   result=e^j*series(xc) */
  j = (int)mpf_get_d(x); /* j should not be zero */
  mpf_set_si(ej, j);
  mpf_sub(xc, x, ej);
  series_exp(result, xc);
  if(debug) printf("mpf_exp j=%d \n", j);
  if(debug) gmp_printf("xc=%50.45Ff \n", xc);
  if(debug) gmp_printf("result=%50.45Ff \n", result);
  if(mpf_cmp(x, zero)>=0)
  {
    for(k=0; k<j; k++) mpf_mul(result, result, e);
    if(debug) gmp_printf("+result=%50.45Ff \n", result);
  }
  else
  {
    for(k=0; k<(-j); k++) mpf_div(result, result, e);
    if(debug) gmp_printf("-result=%50.45Ff \n", result);
  }
} /* end mpf_exp.c */