diff options
Diffstat (limited to 'libm/ldouble/exp10l.c')
| -rw-r--r-- | libm/ldouble/exp10l.c | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/libm/ldouble/exp10l.c b/libm/ldouble/exp10l.c new file mode 100644 index 000000000..b837571b4 --- /dev/null +++ b/libm/ldouble/exp10l.c @@ -0,0 +1,192 @@ +/* exp10l.c + * + * Base 10 exponential function, long double precision + * (Common antilogarithm) + * + * + * + * SYNOPSIS: + * + * long double x, y, exp10l() + * + * y = exp10l( 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). + * The Pade' form + * + * 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) + * + * is used to approximate 10**f. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-4900 30000 1.0e-19 2.7e-20 + * + * ERROR MESSAGES: + * + * message condition value returned + * exp10l underflow x < -MAXL10 0.0 + * exp10l overflow x > MAXL10 MAXNUM + * + * IEEE arithmetic: MAXL10 = 4932.0754489586679023819 + * + */ + +/* +Cephes Math Library Release 2.2: January, 1991 +Copyright 1984, 1991 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +#include <math.h> + +#ifdef UNK +static long double P[] = { + 3.1341179396892496811523E1L, + 4.5618283154904699073999E3L, + 1.3433113468542797218610E5L, + 7.6025447914440301593592E5L, +}; +static long double Q[] = { +/* 1.0000000000000000000000E0,*/ + 4.7705440288425157637739E2L, + 2.9732606548049614870598E4L, + 4.0843697951001026189583E5L, + 6.6034865026929015925608E5L, +}; +/*static long double LOG102 = 3.0102999566398119521373889e-1L;*/ +static long double LOG210 = 3.3219280948873623478703L; +static long double LG102A = 3.01025390625e-1L; +static long double LG102B = 4.6050389811952137388947e-6L; +#endif + + +#ifdef IBMPC +static short P[] = { +0x399a,0x7dc7,0xbc43,0xfaba,0x4003, XPD +0xb526,0xdf32,0xa063,0x8e8e,0x400b, XPD +0x18da,0xafa1,0xc89e,0x832e,0x4010, XPD +0x503d,0x9352,0xe7aa,0xb99b,0x4012, XPD +}; +static short Q[] = { +/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/ +0x947d,0x7855,0xf6ac,0xee86,0x4007, XPD +0x18cf,0x7749,0x368d,0xe849,0x400d, XPD +0x85be,0x2560,0x9f58,0xc76e,0x4011, XPD +0x6d3c,0x80c5,0xca67,0xa137,0x4012, XPD +}; +/* +static short L102[] = {0xf799,0xfbcf,0x9a84,0x9a20,0x3ffd, XPD}; +#define LOG102 *(long double *)L102 +*/ +static short L210[] = {0x8afe,0xcd1b,0x784b,0xd49a,0x4000, XPD}; +#define LOG210 *(long double *)L210 +static short L102A[] = {0x0000,0x0000,0x0000,0x9a20,0x3ffd, XPD}; +#define LG102A *(long double *)L102A +static short L102B[] = {0x8f89,0xf798,0xfbcf,0x9a84,0x3fed, XPD}; +#define LG102B *(long double *)L102B +#endif + +#ifdef MIEEE +static long P[] = { +0x40030000,0xfababc43,0x7dc7399a, +0x400b0000,0x8e8ea063,0xdf32b526, +0x40100000,0x832ec89e,0xafa118da, +0x40120000,0xb99be7aa,0x9352503d, +}; +static long Q[] = { +/* 0x3fff0000,0x80000000,0x00000000, */ +0x40070000,0xee86f6ac,0x7855947d, +0x400d0000,0xe849368d,0x774918cf, +0x40110000,0xc76e9f58,0x256085be, +0x40120000,0xa137ca67,0x80c56d3c, +}; +/* +static long L102[] = {0x3ffd0000,0x9a209a84,0xfbcff799}; +#define LOG102 *(long double *)L102 +*/ +static long L210[] = {0x40000000,0xd49a784b,0xcd1b8afe}; +#define LOG210 *(long double *)L210 +static long L102A[] = {0x3ffd0000,0x9a200000,0x00000000}; +#define LG102A *(long double *)L102A +static long L102B[] = {0x3fed0000,0x9a84fbcf,0xf7988f89}; +#define LG102B *(long double *)L102B +#endif + +static long double MAXL10 = 4.9320754489586679023819e3L; +extern long double MAXNUML; +#ifdef ANSIPROT +extern long double floorl ( long double ); +extern long double ldexpl ( long double, int ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern int isnanl ( long double ); +#else +long double floorl(), ldexpl(), polevll(), p1evll(), isnanl(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif + +long double exp10l(x) +long double x; +{ +long double px, xx; +short n; + +#ifdef NANS +if( isnanl(x) ) + return(x); +#endif +if( x > MAXL10 ) + { +#ifdef INFINITIES + return( INFINITYL ); +#else + mtherr( "exp10l", OVERFLOW ); + return( MAXNUML ); +#endif + } + +if( x < -MAXL10 ) /* Would like to use MINLOG but can't */ + { +#ifndef INFINITIES + mtherr( "exp10l", UNDERFLOW ); +#endif + return(0.0L); + } + +/* Express 10**x = 10**g 2**n + * = 10**g 10**( n log10(2) ) + * = 10**( g + n log10(2) ) + */ +px = floorl( LOG210 * x + 0.5L ); +n = px; +x -= px * LG102A; +x -= px * LG102B; + +/* rational approximation for exponential + * of the fractional part: + * 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) + */ +xx = x * x; +px = x * polevll( xx, P, 3 ); +x = px/( p1evll( xx, Q, 4 ) - px ); +x = 1.0L + ldexpl( x, 1 ); + +/* multiply by power of 2 */ +x = ldexpl( x, n ); +return(x); +} |