/*
 *	+++ NAME +++
 *
 *	 DEXP   Double precision exponential
 *
 *	+++ INDEX +++
 *
 *	 DEXP
 *	 machine independent routines
 *	 math libraries
 *
 *	+++ DESCRIPTION +++
 *
 *	Returns double precision exponential of double precision
 *	floating point number.
 *
 *	+++ USAGE +++
 *
 *	 double dexp(x)
 *	 double x;
 *
 *	+++ REFERENCES +++
 *
 *	Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
 *
 *	Computer Approximations, J.F. Hart et al, John Wiley & Sons,
 *	1968, pp. 96-104.
 *
 *	+++ RESTRICTIONS +++
 *
 *	Inputs greater than LN(MAX_POS_DBLF) result in overflow.
 *	Inputs less than LN(MIN_POS_DBLF) result in underflow.
 *
 *	The maximum relative error for the approximating polynomial
 *	is 10**(-16.4).  However, this assumes exact arithmetic
 *	in the polynomial evaluation.  Additional rounding and
 *	truncation errors may occur as the argument is reduced
 *	to the range over which the polynomial approximation
 *	is valid, and as the polynomial is evaluated using
 *	finite-precision arithmetic.
 *
 *	+++ PROGRAMMER +++
 *
 *	 Fred Fish
 *	 Goodyear Aerospace Corp, Arizona Div.
 *	 (602) 932-7000 work
 *	 (602) 894-6881 home
 *
 *	+++ INTERNALS +++
 *
 *	Computes exponential from:
 *
 *		EXP(X)	=	2**Y  *  2**Z  *  2**W
 *
 *	Where:
 *
 *		Y	=	INT ( X * LOG2(e) )
 *
 *		V	=	16 * FRAC ( X * LOG2(e))
 *
 *		Z	=	(1/16) * INT (V)
 *
 *		W	=	(1/16) * FRAC (V)
 *
 *	Note that:
 *
 *		0 =< V < 16
 *
 *		Z = {0, 1/16, 2/16, ...15/16}
 *
 *		0 =< W < 1/16
 *
 *	Then:
 *
 *		2**Z is looked up in a table of 2**0, 2**1/16, ...
 *
 *		2**W is computed from an approximation:
 *
 *			2**W	=  (A + B) / (A - B)
 *
 *			A	=  C + (D * W * W)
 *
 *			B	=  W * (E + (F * W * W))
 *
 *			C	=  20.8137711965230361973
 *
 *			D	=  1.0
 *
 *			E	=  7.2135034108448192083
 *
 *			F	=  0.057761135831801928
 *
 *		Coefficients are from HART, table #1121, pg 206.
 *
 *		Effective multiplication by 2**Y is done by a
 *		floating point scale with Y as scale argument.
 *
 *	---
 */

/*)LIBRARY
*/

#include <stdio.h>
#include "c:pmluse.h"
#include "pml.h"

# define C  20.8137711965230361973	/* Polynomial approx coeff.	*/
# define D  1.0				/* Polynomial approx coeff.	*/
# define E  7.2135034108448192083	/* Polynomial approx coeff.	*/
# define F  0.057761135831801928	/* Polynomial approx coeff.	*/


/************************************************************************
 *									*
 *	This table is fixed in size and reasonably hardware		*
 *	independent.  The given constants were generated on a 		*
 *	DECSYSTEM 20 using double precision FORTRAN.			*
 *									*
 ************************************************************************
 */

static double fpof2tbl [] = {
    1.00000000000000000000,		/*    2 ** 0/16		*/
    1.04427378242741384020,		/*    2 ** 1/16		*/
    1.09050773266525765930,		/*    2 ** 2/16		*/
    1.13878863475669165390,		/*    2 ** 3/16		*/
    1.18920711500272106640,		/*    2 ** 4/16		*/
    1.24185781207348404890,		/*    2 ** 5/16		*/
    1.29683955465100966610,		/*    2 ** 6/16		*/
    1.35425554693689272850,		/*    2 ** 7/16		*/
    1.41421356237309504880,		/*    2 ** 8/16		*/
    1.47682614593949931110,		/*    2 ** 9/16		*/
    1.54221082540794082350,		/*    2 ** 10/16	*/
    1.61049033194925430820,		/*    2 ** 11/16	*/
    1.68179283050742908600,		/*    2 ** 12/16	*/
    1.75625216037329948340,		/*    2 ** 13/16	*/
    1.83400808640934246360,		/*    2 ** 14/16	*/
    1.91520656139714729380		/*    2 ** 15/16	*/
};

double dexp(x)
double x;
{
    double w, v, a, b, t, wpof2, zpof2, dabs(), dscale(), dint(), dfrac();
    register int y, index;

    if (x > LN_MAXPOSDBL) {
	pmlerr(DEXP_OVERFLOW);
	return (MAX_POS_DBLF);
    } else if (x <= LN_MINPOSDBL) {
	pmlerr(DEXP_UNDERFLOW);
	return (MIN_POS_DBLF);
    } else {
	t = x * LOG2E;
	y = dint(t);
	v = 16.0 * dfrac(t);
	index = dint(dabs(v));
	if (x < 0.0) {
	    zpof2 = 1.0 / fpof2tbl[index];
	} else {
	    zpof2 = fpof2tbl[index];
	}
	w = dfrac(v) / 16.0;
	a = C + (D * w * w);
	b = w * (E + (F * w * w));
	wpof2 = (a + b) / (a - b);
	return ( dscale((wpof2 * zpof2),y) );
    }
}