/************************************************************************
 *                                                                      *
 * A replacement exp() routine for PML library.  An original algorithm  *
 * is not changed but a table of fractionals power of 2 was recalcul-   *
 * ated (believe it or not - this has an influence on a final result).  *
 * Also divisions by power of 2 were replaced by applications of ldexp  *
 * routine which is faster and better preserves precision               *
 *                                                                      *
 * Michal Jaegermann, December 1990                                     *    
 *                                                                      *
 ************************************************************************
 */
  /************************************************************************
 *									*
 *				N O T I C E				*
 *									*
 *			Copyright Abandoned, 1987, Fred Fish		*
 *									*
 *	This previously copyrighted work has been placed into the	*
 *	public domain by the author (Fred Fish) and may be freely used	*
 *	for any purpose, private or commercial.  I would appreciate	*
 *	it, as a courtesy, if this notice is left in all copies and	*
 *	derivative works.  Thank you, and enjoy...			*
 *									*
 *	The author makes no warranty of any kind with respect to this	*
 *	product and explicitly disclaims any implied warranties of	*
 *	merchantability or fitness for any particular purpose.		*
 *									*
 ************************************************************************
 */


/*
 *  FUNCTION
 *
 *	exp   double precision exponential
 *
 *  KEY WORDS
 *
 *	exp
 *	machine independent routines
 *	math libraries
 *
 *  DESCRIPTION
 *
 *	Returns double precision exponential of double precision
 *	floating point number.
 *
 *  USAGE
 *
 *	double exp (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 log(MAXDOUBLE) result in overflow.
 *	Inputs less than log(MINDOUBLE) 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
 *
 *  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.
 *
 */
 
#include <stdio.h>
#include <pmluser.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 computed on Atari ST     *
 *      using integer arithmetic and decimal values for 2**(1/2),       *
 *      2**(1/4) and 2**(1/8) taken from Appendix C of HART et al.      *
 *									*
 ************************************************************************
 */
 
static double fpof2tbl[] = {
    1.0/*0000000000000000000*/,		/*	2 ** 0/16	*/
    1.04427378242741384032,		/*	2 ** 1/16	*/
    1.09050773266525765920,		/*	2 ** 2/16	*/
    1.13878863475669165369,		/*	2 ** 3/16	*/
    1.18920711500272106671,		/*	2 ** 4/16	*/
    1.24185781207348404858,		/*	2 ** 5/16	*/
    1.29683955465100966592,		/*	2 ** 6/16	*/
    1.35425554693689272828,		/*	2 ** 7/16	*/
    1.41421356237309504880,		/*	2 ** 8/16	*/
    1.47682614593949931138,		/*	2 ** 9/16	*/
    1.54221082540794082360,		/*	2 ** 10/16	*/
    1.61049033194925430816,		/*	2 ** 11/16	*/
    1.68179283050742908605,		/*	2 ** 12/16	*/
    1.75625216037329948310,		/*	2 ** 13/16	*/
    1.83400808640934246346,		/*	2 ** 14/16	*/
    1.91520656139714729384		/*	2 ** 15/16	*/
};

static char funcname[] = "exp";


double exp (x)
double x;
{
    register int index;
    auto int y;
    auto double w;
    auto double a;
    auto double b;
    auto double temp;
    auto char * xcptstr;
    extern double dabs ();
    extern double ldexp ();
    extern double modf ();
    auto struct exception xcpt;
 
    DBUG_ENTER (funcname);
    DBUG_3 ("expin", "arg %le", x);
    if (x > LOGE_MAXDOUBLE) {
	xcpt.type = OVERFLOW;
	xcpt.retval = MAXDOUBLE;
	xcptstr = "OVER";
	goto eset;
    }
    if (x <= LOGE_MINDOUBLE) {
	xcpt.type = UNDERFLOW;
	xcpt.retval = 0.0;
	xcptstr = "UNDER";
eset:
	xcpt.name = funcname;
	xcpt.arg1 = x;
	if (!matherr (&xcpt)) {
            errno = ERANGE;
	    fprintf (stderr, "%s: %sFLOW error\n", funcname, xcptstr);
	}
    } else {
	x *= LOG2E;
	w = ldexp(modf(x,&a), 4);
	y = a;
	w = ldexp(modf(w, &temp), -4);
	index = temp;
	
	b = w * w;
	a = C + b;
	b = w * (E + (F * b));
	xcpt.retval = ldexp (((a + b) / (a - b)) *
	    (index < 0 ? ldexp(fpof2tbl[16 + index], -1) : fpof2tbl[index]), y);
    }
    DBUG_3 ("expout", "result %le", xcpt.retval);
    DBUG_RETURN (xcpt.retval);
}