diff options
Diffstat (limited to 'src/math/log.c')
| -rw-r--r-- | src/math/log.c | 84 |
1 files changed, 32 insertions, 52 deletions
diff --git a/src/math/log.c b/src/math/log.c index 98051205..e61e113d 100644 --- a/src/math/log.c +++ b/src/math/log.c @@ -10,7 +10,7 @@ * ==================================================== */ /* log(x) - * Return the logrithm of x + * Return the logarithm of x * * Method : * 1. Argument Reduction: find k and f such that @@ -60,12 +60,12 @@ * to produce the hexadecimal values shown. */ -#include "libm.h" +#include <math.h> +#include <stdint.h> static const double ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ -two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ @@ -76,63 +76,43 @@ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ double log(double x) { - double hfsq,f,s,z,R,w,t1,t2,dk; - int32_t k,hx,i,j; - uint32_t lx; - - EXTRACT_WORDS(hx, lx, x); + union {double f; uint64_t i;} u = {x}; + double_t hfsq,f,s,z,R,w,t1,t2,dk; + uint32_t hx; + int k; + hx = u.i>>32; k = 0; - if (hx < 0x00100000) { /* x < 2**-1022 */ - if (((hx&0x7fffffff)|lx) == 0) - return -two54/0.0; /* log(+-0)=-inf */ - if (hx < 0) - return (x-x)/0.0; /* log(-#) = NaN */ - /* subnormal number, scale up x */ + if (hx < 0x00100000 || hx>>31) { + if (u.i<<1 == 0) + return -1/(x*x); /* log(+-0)=-inf */ + if (hx>>31) + return (x-x)/0.0; /* log(-#) = NaN */ + /* subnormal number, scale x up */ k -= 54; - x *= two54; - GET_HIGH_WORD(hx,x); - } - if (hx >= 0x7ff00000) - return x+x; - k += (hx>>20) - 1023; - hx &= 0x000fffff; - i = (hx+0x95f64)&0x100000; - SET_HIGH_WORD(x, hx|(i^0x3ff00000)); /* normalize x or x/2 */ - k += i>>20; + x *= 0x1p54; + u.f = x; + hx = u.i>>32; + } else if (hx >= 0x7ff00000) { + return x; + } else if (hx == 0x3ff00000 && u.i<<32 == 0) + return 0; + + /* reduce x into [sqrt(2)/2, sqrt(2)] */ + hx += 0x3ff00000 - 0x3fe6a09e; + k += (int)(hx>>20) - 0x3ff; + hx = (hx&0x000fffff) + 0x3fe6a09e; + u.i = (uint64_t)hx<<32 | (u.i&0xffffffff); + x = u.f; + f = x - 1.0; - if ((0x000fffff&(2+hx)) < 3) { /* -2**-20 <= f < 2**-20 */ - if (f == 0.0) { - if (k == 0) { - return 0.0; - } - dk = (double)k; - return dk*ln2_hi + dk*ln2_lo; - } - R = f*f*(0.5-0.33333333333333333*f); - if (k == 0) - return f - R; - dk = (double)k; - return dk*ln2_hi - ((R-dk*ln2_lo)-f); - } + hfsq = 0.5*f*f; s = f/(2.0+f); - dk = (double)k; z = s*s; - i = hx - 0x6147a; w = z*z; - j = 0x6b851 - hx; t1 = w*(Lg2+w*(Lg4+w*Lg6)); t2 = z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); - i |= j; R = t2 + t1; - if (i > 0) { - hfsq = 0.5*f*f; - if (k == 0) - return f - (hfsq-s*(hfsq+R)); - return dk*ln2_hi - ((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f); - } else { - if (k == 0) - return f - s*(f-R); - return dk*ln2_hi - ((s*(f-R)-dk*ln2_lo)-f); - } + dk = k; + return s*(hfsq+R) + dk*ln2_lo - hfsq + f + dk*ln2_hi; } |