/* 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  1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
 * of degree 2/3 is used to approximate exp(f) in the basic
 * interval [-0.5, 0.5].
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       +- 88       50000       2.8e-17     7.0e-18
 *    IEEE      +- 708      40000       2.0e-16     5.6e-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         INFINITY
 *
 */

/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/


/* Exponential function */

#include "mconf.h"

#ifdef UNK

static double P[] = {
 1.26177193074810590878E-4,
 3.02994407707441961300E-2,
 9.99999999999999999910E-1,
};
static double Q[] = {
 3.00198505138664455042E-6,
 2.52448340349684104192E-3,
 2.27265548208155028766E-1,
 2.00000000000000000009E0,
};
static double C1 = 6.93145751953125E-1;
static double C2 = 1.42860682030941723212E-6;
#endif

#ifdef DEC
static unsigned short P[] = {
0035004,0047156,0127442,0057502,
0036770,0033210,0063121,0061764,
0040200,0000000,0000000,0000000,
};
static unsigned short Q[] = {
0033511,0072665,0160662,0176377,
0036045,0070715,0124105,0132777,
0037550,0134114,0142077,0001637,
0040400,0000000,0000000,0000000,
};
static unsigned short sc1[] = {0040061,0071000,0000000,0000000};
#define C1 (*(double *)sc1)
static unsigned short sc2[] = {0033277,0137216,0075715,0057117};
#define C2 (*(double *)sc2)
#endif

#ifdef IBMPC
static unsigned short P[] = {
0x4be8,0xd5e4,0x89cd,0x3f20,
0x2c7e,0x0cca,0x06d1,0x3f9f,
0x0000,0x0000,0x0000,0x3ff0,
};
static unsigned short Q[] = {
0x5fa0,0xbc36,0x2eb6,0x3ec9,
0xb6c0,0xb508,0xae39,0x3f64,
0xe074,0x9887,0x1709,0x3fcd,
0x0000,0x0000,0x0000,0x4000,
};
static unsigned short sc1[] = {0x0000,0x0000,0x2e40,0x3fe6};
#define C1 (*(double *)sc1)
static unsigned short sc2[] = {0xabca,0xcf79,0xf7d1,0x3eb7};
#define C2 (*(double *)sc2)
#endif

#ifdef MIEEE
static unsigned short P[] = {
0x3f20,0x89cd,0xd5e4,0x4be8,
0x3f9f,0x06d1,0x0cca,0x2c7e,
0x3ff0,0x0000,0x0000,0x0000,
};
static unsigned short Q[] = {
0x3ec9,0x2eb6,0xbc36,0x5fa0,
0x3f64,0xae39,0xb508,0xb6c0,
0x3fcd,0x1709,0x9887,0xe074,
0x4000,0x0000,0x0000,0x0000,
};
static unsigned short sc1[] = {0x3fe6,0x2e40,0x0000,0x0000};
#define C1 (*(double *)sc1)
static unsigned short sc2[] = {0x3eb7,0xf7d1,0xcf79,0xabca};
#define C2 (*(double *)sc2)
#endif

#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double polevl(), p1evl(), floor(), ldexp();
int isnan(), isfinite();
#endif
extern double LOGE2, LOG2E, MAXLOG, MINLOG, MAXNUM;
#ifdef INFINITIES
extern double INFINITY;
#endif

double exp(x)
double x;
{
double px, xx;
int n;

#ifdef NANS
if( isnan(x) )
 return(x);
#endif
if( x > MAXLOG)
 {
#ifdef INFINITIES
 return( INFINITY );
#else
 mtherr( "exp", OVERFLOW );
 return( MAXNUM );
#endif
 }

if( x < MINLOG )
 {
#ifndef INFINITIES
 mtherr( "exp", UNDERFLOW );
#endif
 return(0.0);
 }

/* Express e**x = e**g 2**n
 *   = e**g e**( n loge(2) )
 *   = e**( g + n loge(2) )
 */
px = floor( LOG2E * x + 0.5 ); /* floor() truncates toward -infinity. */
n = px;
x -= px * C1;
x -= px * 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 = 1.0 + 2.0 * x;

/* multiply by power of 2 */
x = ldexp( x, n );
return(x);
}