/* my_exp.c  normalize and use series */
/* e^x = e^j * e^xc   xc = x-j if |xc|>1/2 sqrt(e) */
#include <stdio.h>
#include <math.h>

#define e  2.7182818284590452353602874713526624977572
#define abs(x) ((x)<0.0?(-(x)):(x))

double ser_exp(double x)
{
  double y;    /* sum of exp series on xc */
  double xp;   /* power of x */
  double fact; /* factorial */
  int i;
  int n = 28;  /* need 0.5^n / n! < 0.5^52 */

  if(abs(x)>0.5) printf("series_exp x too big = %g \n", x);
  xp = 1.0;
  fact = 1.0;
  y = 1.0;
  for(i=1; i<n; i++)
  {
    fact=fact*(double)i;
    xp = xp * x;
    y = y + xp/fact;
  }
  return y;
} /* end ser_exp */

double my_exp(double x)
{
  double y;    /* sum of exp series on xc */
  double xc;   /* reduced x */
  double xp;   /* power of xc */
  double fact; /* factorial */
  double ep;   /* power of e, e^2^j */
  double p2;   /* power of 2 */
  int i, j, k;

  if(abs(x)<=0.5) return ser_exp(x); 
  if(x>0.0)
  {
    ep = 1.0;
    j=0;
    while(x>(double)(j+1))
    { 
      ep = ep * e;
      j++;
    }
    xc = x-(double)j;
    printf("ep=%g, j=%d, xc=%g \n", ep, j, xc);
    if(xc<=0.5) return ep*ser_exp(xc); 
    else       return ep*sqrt(e)*ser_exp(xc-0.5);
  }
  else /* x<0.0 */
  {
    ep = 1.0;
    j=0;
    while(x<(double)(-(j+1)))
    { 
      ep = ep / e;
      j++;
    }
    xc = x+(double)j;
    printf("ep=%g, j=%d, xc=%g \n", ep, j, xc);
    if(xc>=-0.5) return ep*ser_exp(xc); 
    else       return ep*ser_exp(xc+0.5)/sqrt(e);
  }

} /* end exp */

void test_exp(double x)
{
  double y, ym, diff;

  y  = my_exp(x);
  ym = exp(x);
  diff = (y - ym)/ym; /* relative error */
  printf("                            x=%g \n",x);
  printf("%30.3f  my_exp \n", y);
  printf("%30.3f  exp \n", ym);
  printf("                            rel err %g \n\n", diff);
} /* end test_exp */

int main(int argc, char *argv[])
{
  double x, y, ym, diff;
  int i;

  printf("my_exp.c compare to math.h exp \n");
  for(i=0; i<52; i++)
  {
    x  = (double)i;
    y  = my_exp(x);
    ym = exp(x);
    diff = (y - ym)/ym; /* relative error */
    printf("%30.3f  my_exp(%d) \n", y, i);
    printf("%30.3f  exp \n", ym);
    printf("                            rel err %g \n\n", diff);
  }
  test_exp(0.1);
  test_exp(0.4);
  test_exp(0.6);
  test_exp(0.99);
  test_exp(1.01);
  test_exp(2.8);
  test_exp(9.0);
  test_exp(-0.1);
  test_exp(-0.4);
  test_exp(-0.6);
  test_exp(-0.99);
  test_exp(-1.01);
  test_exp(-2.8);
  test_exp(-9.0);
  return 0;
} /* end main of my_exp.c */