fluencelabs/musl
1
0
Fork 0
mirror of https://github.com/fluencelabs/musl synced 2025-05-29 07:31:53 +00:00
Code Issues Projects Releases Wiki Activity

math: expl.c cleanup

Browse Source 
raise overflow and underflow when necessary, fix various comments.
This commit is contained in:
Szabolcs Nagy 2012-11-18 03:49:16 +01:00
parent ab1772c597
commit e93a0fe49d
1 changed files with 19 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
src/math/expl.c
View File

@ -35,7 +35,7 @@
* x k f * x k f
* e = 2 e. * e = 2 e.
* *
* A Pade' form of degree 2/3 is used to approximate exp(f) - 1 * A Pade' form of degree 5/6 is used to approximate exp(f) - 1
* in the basic range [-0.5 ln 2, 0.5 ln 2]. * in the basic range [-0.5 ln 2, 0.5 ln 2].
* *
* *
@ -86,42 +86,37 @@ static const long double Q[4] = {
2.0000000000000000000897E0L, 2.0000000000000000000897E0L,
}; };
static const long double static const long double
C1 = 6.9314575195312500000000E-1L, LN2HI = 6.9314575195312500000000E-1L,
C2 = 1.4286068203094172321215E-6L, LN2LO = 1.4286068203094172321215E-6L,
MAXLOGL = 1.1356523406294143949492E4L, LOG2E = 1.4426950408889634073599E0L;
MINLOGL = -1.13994985314888605586758E4L,
LOG2EL = 1.4426950408889634073599E0L;
long double expl(long double x) long double expl(long double x)
{ {
long double px, xx; long double px, xx;
int n; int k;
if (isnan(x)) if (isnan(x))
return x; return x;
if (x > MAXLOGL) if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */
return INFINITY; return x * 0x1p16383L;
if (x < MINLOGL) if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */
return 0.0; return 0x1p-10000L * 0x1p-10000L;
/* Express e**x = e**g 2**n /* Express e**x = e**f 2**k
* = e**g e**(n loge(2)) * = e**(f + k ln(2))
* = e**(g + n loge(2))
*/ */
px = floorl(LOG2EL * x + 0.5); /* floor() truncates toward -infinity. */ px = floorl(LOG2E * x + 0.5);
n = px; k = px;
x -= px * C1; x -= px * LN2HI;
x -= px * C2; x -= px * LN2LO;
/* rational approximation for exponential /* rational approximation of the fractional part:
* of the fractional part: * e**x = 1 + 2x P(x**2)/(Q(x**2) - x P(x**2))
* e**x = 1 + 2x P(x**2)/(Q(x**2) - P(x**2))
*/ */
xx = x * x; xx = x * x;
px = x * __polevll(xx, P, 2); px = x * __polevll(xx, P, 2);
x = px/(__polevll(xx, Q, 3) - px); x = px/(__polevll(xx, Q, 3) - px);
x = 1.0 + 2.0 * x; x = 1.0 + 2.0 * x;
x = scalbnl(x, n); return scalbnl(x, k);
return x;
} }
#endif #endif