123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171 |
- /* 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 "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 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);
- }
|