From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/exp2.c | 183 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100644 libm/double/exp2.c (limited to 'libm/double/exp2.c') diff --git a/libm/double/exp2.c b/libm/double/exp2.c new file mode 100644 index 000000000..be5bdfd0c --- /dev/null +++ b/libm/double/exp2.c @@ -0,0 +1,183 @@ +/* exp2.c + * + * Base 2 exponential function + * + * + * + * SYNOPSIS: + * + * double x, y, exp2(); + * + * y = 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 exp.c for comments on error amplification. + * + * + * ERROR MESSAGES: + * + * message condition value returned + * exp underflow x < -MAXL2 0.0 + * 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 + +#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 floor ( double ); +extern double ldexp ( double, int ); +extern int isnan ( double ); +extern int isfinite ( double ); +#else +double polevl(), p1evl(), floor(), ldexp(); +int isnan(), isfinite(); +#endif +#ifdef INFINITIES +extern double INFINITY; +#endif +extern double MAXNUM; + +double 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( "exp2", OVERFLOW ); + return( MAXNUM ); +#endif + } + +if( x < MINL2 ) + { +#ifndef INFINITIES + mtherr( "exp2", UNDERFLOW ); +#endif + return(0.0); + } + +xx = x; /* save x */ +/* separate into integer and fractional parts */ +px = floor(x+0.5); +n = px; +x = x - px; + +/* rational approximation + * 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 + ldexp( x, 1 ); + +/* scale by power of 2 */ +x = ldexp( x, n ); +return(x); +} -- cgit v1.2.3