1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
|
/* 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 NPY_INFINITY
*
* For IEEE arithmetic, MAXL2 = 1024.
*/
/*
* Cephes Math Library Release 2.3: March, 1995
* Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
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
double exp2(double x)
{
double px, xx;
short n;
if (cephes_isnan(x))
return (x);
if (x > MAXL2) {
return (NPY_INFINITY);
}
if (x < MINL2) {
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);
}
|