/*
 *      e x p 
 */

#define __LIBRARY__

#include <math.h>
#include <limits.h>
#include <errno.h>
#include <debug.h>

static double expcoeff[] =
       {1.0000000000000000000E0,	/* 1/0! */
		1.0000000000000000000E0,	/* 1/1! */
		5.0000000000000000000E-1,	/* 1/2! and so on */
		1.6666666666666666667E-1,
		4.1666666666666666667E-2,
		8.3333333333333333333E-3,
		1.3888888888888888889E-3,
		1.9841269841269841270E-4,
		2.4801587301587301587E-5,
		2.7557319223985890652E-6,
		2.7557319223985890652E-7,
		2.5052108385441718775E-8,
		2.0876756987868098979E-9,	/* coeffs corrected to 20 s.f. 3/92 */
		1.6059043836821614599E-10,  /* and extended S E Peterson */ 
		1.1470745597729724713E-11,
		7.6471637318198164759E-13,
		4.7794773323873852974E-14,
		2.8114572543455207631E-15,
		1.5619206968586226462E-16,
		8.2206352466243297170E-18,
		4.1103176233121648585E-19,
       };

#define	NCOEFF	((sizeof expcoeff) / sizeof(double))

double exp(x)
double x;
{
    double ipart, frac;

    frac = modf(x * M_LOG2E, &ipart);
    if(frac > 0.5) {
        DBG(("EXP",0x21,"fraction > +0.5"));
        ipart += 1.0;
        frac -= 1.0;
    }
    else if(frac < -0.5) {
        DBG(("EXP",0x21,"fraction < -0.5"));
        ipart -= 1.0;
        frac += 1.0;
    }

    if(ipart < (double) INT_MIN) {		/* really small */
        DBG(("EXP",0x21,"ipart < INT_MIN"));
        return(0);
    }
    else if(ipart > (double) INT_MAX) {	/* really big */
        DBG(("EXP",0x21,"ipart > INT_MAX"));
        ipart = (double) INT_MAX;
    }
    return (ldexp(_poly(frac * M_LN2, expcoeff, NCOEFF), (int) ipart));
}