/*							exp10f.c
 *
 *	Base 10 exponential function
 *      (Common antilogarithm)
 *
 *
 *
 * SYNOPSIS:
 *
 * float x, y, exp10f();
 *
 * y = exp10f( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns 10 raised to the x power.
 *
 * Range reduction is accomplished by expressing the argument
 * as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
 * A polynomial approximates 10**f.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -38,+38     100000      9.8e-8      2.8e-8
 *
 * ERROR MESSAGES:
 *
 *   message         condition      value returned
 * exp10 underflow    x < -MAXL10        0.0
 * exp10 overflow     x > MAXL10       MAXNUM
 *
 * IEEE single arithmetic: MAXL10 = 38.230809449325611792.
 *
 */

/*
Cephes Math Library Release 2.2:  June, 1992
Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/


#include <math.h>

static float P[] = {
 2.063216740311022E-001,
 5.420251702225484E-001,
 1.171292686296281E+000,
 2.034649854009453E+000,
 2.650948748208892E+000,
 2.302585167056758E+000
};

/*static float LOG102 = 3.01029995663981195214e-1;*/
static float LOG210 = 3.32192809488736234787e0;
static float LG102A = 3.00781250000000000000E-1;
static float LG102B = 2.48745663981195213739E-4;
static float MAXL10 = 38.230809449325611792;




extern float MAXNUMF;

float floorf(float), ldexpf(float, int), polevlf(float, float *, int);

float exp10f(float xx)
{
float x, px, qx;
short n;

x = xx;
if( x > MAXL10 )
	{
	mtherr( "exp10f", OVERFLOW );
	return( MAXNUMF );
	}

if( x < -MAXL10 )	/* Would like to use MINLOG but can't */
	{
	mtherr( "exp10f", UNDERFLOW );
	return(0.0);
	}

/* The following is necessary because range reduction blows up: */
if( x == 0 )
	return(1.0);

/* Express 10**x = 10**g 2**n
 *   = 10**g 10**( n log10(2) )
 *   = 10**( g + n log10(2) )
 */
px = x * LOG210;
qx = floorf( px + 0.5 );
n = qx;
x -= qx * LG102A;
x -= qx * LG102B;

/* rational approximation for exponential
 * of the fractional part:
 * 10**x - 1  =  2x P(x**2)/( Q(x**2) - P(x**2) )
 */
px = 1.0 + x * polevlf( x, P, 5 );

/* multiply by power of 2 */
x = ldexpf( px, n );

return(x);
}