/*
 * 				e x p . c
 */

/*)LIBRARY
*/

#ifdef	DOCUMENTATION
title	exp	exponential function
index	exponential function
usage
	.s
	double x, f, exp();
	.br
	f = exp(x);
	.s
description
	.s
	Returns exponential of x.
	.s
diagnostics
	.s
	If the argument is large enough to cause overflow in the result
	the message 'exp overflow', followed by the value of the argument,
	is written to stderr.  A value of HUGE is returned.
	.br
	If the argument would produce a result exponent
	too small to accomodate, a value of zero is returned. A warning 
	message 'exp underflow', followed by the value of the argument,
	is written to stderr, if warnings are enabled (normally not).
	.s
internal
	.s
	Algorithm from Cody and Waite pp. 60-83.
author
	.s
	Hamish Ross.
	.s
date
	.s
	28-Dec-84
#endif

#include <math.h>

#define EPS2 6.93889401e-18	/* exp(eps) = 1.0 to m/c precision	*/
#define C1 0.693359375		/* C1 + C2 should represent log 2 to	*/
#define C2 -2.12194440054690583e-4	/* more than machine precision	*/

static double p0 = 0.249999999999999993e0;
static double p1 = 0.694360001511792852e-2;
static double p2 = 0.165203300268279130e-4;
static double q0 = 0.5;
static double q1 = 0.555538666969001188e-1;
static double q2 = 0.495862884905441294e-3;
double exp(x)
double x;
{
    double g, gp, q, xn, z, ldexp();
    int n;

/* check for exponent too large, if so print argument and return large no */
    if (x >= LOG_HUGE){
	cmemsg(FP_LEXP, &x);
	return(HUGE);
    }
/* check for exponent underflow, if so call error system (set at warning
   level) and return zero */
    if (x < LOG_TINY){
	cmemsg(FP_SEXP, &x);
	return(0.0);
    }
    if (x < EPS2 && x > -EPS2)
	return(1.0);
    z = x * LOGB2E;
    if (z < 0.0)
	n = z - 0.5;
    else
	n = z + 0.5;
    xn = n;
    g = (x - xn * C1) - xn * C2;
    z = g * g;
    gp = ((p2 * z + p1) * z + p0) * g;
    q = (q2 * z + q1) * z + q0;
    return(ldexp(0.5 + gp / (q - gp), ++n));
}