/*							exp.c
 *
 *	Exponential function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, exp();
 *
 * y = exp( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns e (2.71828...) raised to the x power.
 *
 * Range reduction is accomplished by separating the argument
 * into an integer k and fraction f such that
 *
 *     x    k  f
 *    e  = 2  e.
 *
 * A Pade' form of degree 2/3 is used to approximate exp(f) - 1
 * in the basic range [-0.5 ln 2, 0.5 ln 2].
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       0, MAXLOG   38000       3.0e-17     6.2e-18
 *    IEEE      +- MAXLOG   35000       2.3e-16     5.7e-17
 *
 *
 * Error amplification in the exponential function can be
 * a serious matter.  The error propagation involves
 * exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
 * which shows that a 1 lsb error in representing X produces
 * a relative error of X times 1 lsb in the function.
 * While the routine gives an accurate result for arguments
 * that are exactly represented by a double precision
 * computer number, the result contains amplified roundoff
 * error for large arguments not exactly represented.
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * exp underflow    x < MINLOG         0.0
 * exp overflow     x > MAXLOG         MAXNUM
 *
 */

/*
Cephes Math Library Release 2.1:  January, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


/*	Exponential function	*/

#include "mconf.h"
static char fname[] = {"exp"};

#ifdef UNK

static double P[] = {
 1.26183092834458542160E-4,
 3.02996887658430129200E-2,
 1.00000000000000000000E0
};
static double Q[] = {
3.00227947279887615146E-6,
2.52453653553222894311E-3,
2.27266044198352679519E-1,
2.00000000000000000005E0
};
static double C1 = 6.9335937500000000000E-1;
static double C2 = 2.1219444005469058277E-4;
#endif

#ifdef DEC
static short P[] = {
0035004,0050004,0016315,0134545,
0036770,0033415,0105201,0034462,
0040200,0000000,0000000,0000000
};
static short Q[] = {
0033511,0075304,0141275,0061006,
0036045,0071261,0155620,0021143,
0037550,0134156,0006512,0174363,
0040400,0000000,0000000,0000000
};
static short sc1[] = {0040061,0100000,0000000,0000000};
#define C1 (*(double *)sc1)
static short sc2[] = {0035136,0100202,0161410,0062503};
#define C2 (*(double *)sc2)
#endif

#ifdef IBMPC
static short P[] = {
0xb72d,0x8399,0x8a00,0x3f20,
0x2726,0xb150,0x06e1,0x3f9f,
0x0000,0x0000,0x0000,0x3ff0
};
static short Q[] = {
0xac41,0x9857,0x2f58,0x3ec9,
0x044c,0x3b72,0xae56,0x3f64,
0x5f1e,0xc1a9,0x170d,0x3fcd,
0x0000,0x0000,0x0000,0x4000
};
static short sc1[] = {0x0000,0x0000,0x3000,0x3fe6};
#define C1 (*(double *)sc1)
static short sc2[] = {0x0ca8,0x5c61,0xd010,0x3f2b};
#define C2 (*(double *)sc2)
#endif

#ifdef MIEEE
static short P[] = {
0x3f20,0x8a00,0x8399,0xb72d,
0x3f9f,0x06e1,0xb150,0x2726,
0x3ff0,0x0000,0x0000,0x0000
};
static short Q[] = {
0x3ec9,0x2f58,0x9857,0xac41,
0x3f64,0xae56,0x3b72,0x044c,
0x3fcd,0x170d,0xc1a9,0x5f1e,
0x4000,0x0000,0x0000,0x0000
};
static short sc1[] = {
0x3fe6,0x3000,0x0000,0x0000
};
#define C1 (*(double *)sc1)
static short sc2[] = {
0x3f2b,0xd010,0x5c61,0x0ca8
};
#define C2 (*(double *)sc2)
#endif

extern double LOGE2, LOG2E, MAXLOG, MINLOG, MAXNUM;

double exp(x)
double x;
{
double px, qx, xx;
int n;
double polevl();
double floor(), ldexp();

if( x > MAXLOG)
	{
	mtherr( fname, OVERFLOW );
	return( MAXNUM );
	}

if( x < MINLOG )
	{
	mtherr( fname, UNDERFLOW );
	return(0.0);
	}

/* Express e**x = e**g 2**n
 *   = e**g e**( n loge(2) )
 *   = e**( g + n loge(2) )
 */
px = x * LOG2E;
qx = floor( px + 0.5 ); /* floor() truncates toward -infinity. */
n = qx;
x -= qx * C1;
x += qx * C2;


/* rational approximation for exponential
 * of the fractional part:
 * e**x - 1  =  2x P(x**2)/( Q(x**2) - P(x**2) )
 */
xx = x * x;
px = x * polevl( xx, P, 2 );
x =  px/( polevl( xx, Q, 3 ) - px );
x = ldexp( x, 1 );
x =  x + 1.0;
x = ldexp( x, n );
return(x);
}