diff options
author | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
---|---|---|
committer | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
commit | 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch) | |
tree | 3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/exp2.c | |
parent | c117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff) |
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD).
-Erik
Diffstat (limited to 'libm/double/exp2.c')
-rw-r--r-- | libm/double/exp2.c | 183 |
1 files changed, 0 insertions, 183 deletions
diff --git a/libm/double/exp2.c b/libm/double/exp2.c deleted file mode 100644 index be5bdfd0c..000000000 --- a/libm/double/exp2.c +++ /dev/null @@ -1,183 +0,0 @@ -/* 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 <math.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 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); -} |