/*							md_exp2.c
 *
 *	Base 2 exponential function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, md_exp2();
 *
 * y = md_exp2( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns 2 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
 *    2  = 2  2.
 *
 * A Pade' form
 *
 *   1 + 2x P(x**2) / (Q(x**2) - x P(x**2) )
 *
 * approximates 2**x in the basic range [-0.5, 0.5].
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE    -1022,+1024   30000       1.8e-16     5.4e-17
 *
 *
 * See md_exp.c for comments on error amplification.
 *
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * md_exp underflow    x < -MAXL2        0.0
 * md_exp overflow     x > MAXL2         MAXNUM
 *
 * For DEC arithmetic, MAXL2 = 127.
 * For IEEE arithmetic, MAXL2 = 1024.
 */


/*
Cephes Math Library Release 2.8:  June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/



#include "mconf.h"

#ifdef UNK
static double P[] = {
 2.30933477057345225087E-2,
 2.02020656693165307700E1,
 1.51390680115615096133E3,
};
static double Q[] = {
/* 1.00000000000000000000E0,*/
 2.33184211722314911771E2,
 4.36821166879210612817E3,
};
#define MAXL2 1024.0
#define MINL2 -1024.0
#endif

#ifdef DEC
static unsigned short P[] = {
0036675,0027102,0122327,0053227,
0041241,0116724,0115412,0157355,
0042675,0036404,0101733,0132226,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0042151,0027450,0077732,0160744,
0043210,0100661,0077550,0056560,
};
#define MAXL2 127.0
#define MINL2 -127.0
#endif

#ifdef IBMPC
static unsigned short P[] = {
0xead3,0x549a,0xa5c8,0x3f97,
0x5bde,0x9361,0x33ba,0x4034,
0x7693,0x907b,0xa7a0,0x4097,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x5c3c,0x0ffb,0x25e5,0x406d,
0x0bae,0x2fed,0x1036,0x40b1,
};
#define MAXL2 1024.0
#define MINL2 -1022.0
#endif

#ifdef MIEEE
static unsigned short P[] = {
0x3f97,0xa5c8,0x549a,0xead3,
0x4034,0x33ba,0x9361,0x5bde,
0x4097,0xa7a0,0x907b,0x7693,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x406d,0x25e5,0x0ffb,0x5c3c,
0x40b1,0x1036,0x2fed,0x0bae,
};
#define MAXL2 1024.0
#define MINL2 -1022.0
#endif

#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double md_floor ( double );
extern double md_ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double polevl(), p1evl(), md_floor(), md_ldexp();
int isnan(), isfinite();
#endif
#ifdef INFINITIES
extern double INFINITY;
#endif
extern double MAXNUM;

double md_exp2(x)
double x;
{
double px, xx;
short n;

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

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

xx = x;	/* save x */
/* separate into integer and fractional parts */
px = md_floor(x+0.5);
n = px;
x = x - px;

/* rational approximation
 * md_exp2(x) = 1 +  2xP(xx)/(Q(xx) - P(xx))
 * where xx = x**2
 */
xx = x * x;
px = x * polevl( xx, P, 2 );
x =  px / ( p1evl( xx, Q, 2 ) - px );
x = 1.0 + md_ldexp( x, 1 );

/* scale by power of 2 */
x = md_ldexp( x, n );
return(x);
}