/***
*exp.c - exponential
*
*	Copyright (c) 1991-1991, Microsoft Corporation.	All rights reserved.
*
*Purpose:
*   Compute exp(x)
*
*Revision History:
*   8-15-91	GDP	written
*  12-21-91	GDP	support IEEE exceptions
*  02-03-92	GDP	added _exphlp for use by exp, sinh, and cosh
*
*******************************************************************************/
#include <math.h>
#include <trans.h>

double _exphlp(double, int *);

/*
 * Thresholds for over/underflow that results in an adjusted value
 * too big/small to be represented as a double.
 * OVFX: ln(XMAX * 2^IEEE_ADJ)
 * UFLX: ln(XIN * 2^(-IEEE_ADJ)
 */

static _dbl const ovfx = {SET_DBL(0x409bb9d3, 0xbeb8c86b)};
static _dbl const uflx = {SET_DBL(0xc09bb448, 0x2df909dd)};

#define OVFX	ovfx.dbl
#define UFLX	uflx.dbl


static double const  EPS    =  5.16987882845642297e-26;	  /* 2^(-53) / 2 */
static double const  LN2INV =  1.442695040889634074;	  /* 1/ln(2) */
static double const  C1	    =  0.693359375000000000;
static double const  C2	    = -2.1219444005469058277e-4;

/* constants for the rational approximation */
static double const p0 = 0.249999999999999993e+0;
static double const p1 = 0.694360001511792852e-2;
static double const p2 = 0.165203300268279130e-4;
static double const q0 = 0.500000000000000000e+0;
static double const q1 = 0.555538666969001188e-1;
static double const q2 = 0.495862884905441294e-3;

#define P(z)  ( (p2 * (z) + p1) * (z) + p0 )
#define Q(z)  ( (q2 * (z) + q1) * (z) + q0 )

/***
*double exp(double x) - exponential
*
*Purpose:
*   Compute the exponential of a number.
*   The algorithm (reduction / rational approximation) is
*   taken from Cody & Waite.
*
*Entry:
*
*Exit:
*
*Exceptions: O, U, P, I
*
*******************************************************************************/

double exp (double x)
{
    unsigned int savedcw;
    int n, newexp;
    double result;
    double xn;

    /* save user fp control word */
    savedcw = _maskfp();

    if (IS_D_SPECIAL(x)){
	switch (_sptype(x)) {
	case T_PINF:
	    RETURN(savedcw,x);
	case T_NINF:
	    RETURN(savedcw,0.0);
	case T_QNAN:
	    return _handle_qnan1(OP_EXP, x, savedcw);
	default: //T_SNAN
	    return _except1(FP_I, OP_EXP, x, _s2qnan(x), savedcw);
	}
    }

    if (x == 0.0) {
	RETURN(savedcw, 1.0);
    }

    if (x > OVFX) {

	// even after scaling the exponent of the result,
	// it is still too large.
	// Deliver infinity to the trap handler

	return _except1(FP_O | FP_P, OP_EXP, x, D_INF, savedcw);
    }

    if (x < UFLX) {

	// even after scaling the exponent of the result,
	// it is still too small.
	// Deliver 0 to the trap handler

	return _except1(FP_U | FP_P, OP_EXP, x, 0.0, savedcw);
    }

    if (ABS(x) < EPS) {
	result = 1.0;
    }

    else {
	result = _exphlp(x, &newexp);
	if (newexp > MAXEXP) {
	    result = _set_exp(result, newexp-IEEE_ADJUST);
	    return _except1(FP_O | FP_P, OP_EXP, x, result, savedcw);
	}
	else if (newexp < MINEXP) {
	    result = _set_exp(result, newexp+IEEE_ADJUST);
	    return _except1(FP_U | FP_P, OP_EXP, x, result, savedcw);
	}
	else
	    result = _set_exp(result, newexp);
    }

    RETURN_INEXACT1(OP_EXP, x, result, savedcw);
}




/***
*double _exphlp(double x, int * pnewexp) - exp helper routine
*
*Purpose:
*   Provide the mantissa and  the exponent of e^x
*
*Entry:
*   x : a (non special) double precision number
*
*Exit:
*   *newexp: the exponent of e^x
*   return value: the mantissa m of e^x scaled by a factor
*		  (the value of this factor has no significance.
*		   The mantissa can be obtained with _set_exp(m, 0).
*
*   _set_exp(m, *pnewexp) may be used for constructing the final
*   result, if it is within the representable range.
*
*Exceptions:
*   No exceptions are raised by this function
*
*******************************************************************************/



double _exphlp(double x, int * pnewexp)
{

    double xn;
    double g,z,gpz,qz,rg;
    int n;

    xn = _frnd(x * LN2INV);
    n = (int) xn;

    /* assume guard digit is present */
    g = (x - xn * C1) - xn * C2;
    z = g*g;
    gpz = g * P(z);
    qz = Q(z);
    rg = 0.5 + gpz/(qz-gpz);

    n++;

    *pnewexp = _get_exp(rg) + n;
    return rg;
}
