| 1 |
/* $MidnightBSD$ */ |
| 2 |
|
| 3 |
/* @(#)e_log10.c 1.3 95/01/18 */ |
| 4 |
/* |
| 5 |
* ==================================================== |
| 6 |
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 7 |
* |
| 8 |
* Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 9 |
* Permission to use, copy, modify, and distribute this |
| 10 |
* software is freely granted, provided that this notice |
| 11 |
* is preserved. |
| 12 |
* ==================================================== |
| 13 |
*/ |
| 14 |
|
| 15 |
#include <sys/cdefs.h> |
| 16 |
__FBSDID("$FreeBSD: stable/10/lib/msun/src/e_log10.c 251292 2013-06-03 09:14:31Z das $"); |
| 17 |
|
| 18 |
/* |
| 19 |
* Return the base 10 logarithm of x. See e_log.c and k_log.h for most |
| 20 |
* comments. |
| 21 |
* |
| 22 |
* log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2) |
| 23 |
* in not-quite-routine extra precision. |
| 24 |
*/ |
| 25 |
|
| 26 |
#include <float.h> |
| 27 |
|
| 28 |
#include "math.h" |
| 29 |
#include "math_private.h" |
| 30 |
#include "k_log.h" |
| 31 |
|
| 32 |
static const double |
| 33 |
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ |
| 34 |
ivln10hi = 4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */ |
| 35 |
ivln10lo = 2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */ |
| 36 |
log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */ |
| 37 |
log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ |
| 38 |
|
| 39 |
static const double zero = 0.0; |
| 40 |
static volatile double vzero = 0.0; |
| 41 |
|
| 42 |
double |
| 43 |
__ieee754_log10(double x) |
| 44 |
{ |
| 45 |
double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2; |
| 46 |
int32_t i,k,hx; |
| 47 |
u_int32_t lx; |
| 48 |
|
| 49 |
EXTRACT_WORDS(hx,lx,x); |
| 50 |
|
| 51 |
k=0; |
| 52 |
if (hx < 0x00100000) { /* x < 2**-1022 */ |
| 53 |
if (((hx&0x7fffffff)|lx)==0) |
| 54 |
return -two54/vzero; /* log(+-0)=-inf */ |
| 55 |
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */ |
| 56 |
k -= 54; x *= two54; /* subnormal number, scale up x */ |
| 57 |
GET_HIGH_WORD(hx,x); |
| 58 |
} |
| 59 |
if (hx >= 0x7ff00000) return x+x; |
| 60 |
if (hx == 0x3ff00000 && lx == 0) |
| 61 |
return zero; /* log(1) = +0 */ |
| 62 |
k += (hx>>20)-1023; |
| 63 |
hx &= 0x000fffff; |
| 64 |
i = (hx+0x95f64)&0x100000; |
| 65 |
SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */ |
| 66 |
k += (i>>20); |
| 67 |
y = (double)k; |
| 68 |
f = x - 1.0; |
| 69 |
hfsq = 0.5*f*f; |
| 70 |
r = k_log1p(f); |
| 71 |
|
| 72 |
/* See e_log2.c for most details. */ |
| 73 |
hi = f - hfsq; |
| 74 |
SET_LOW_WORD(hi,0); |
| 75 |
lo = (f - hi) - hfsq + r; |
| 76 |
val_hi = hi*ivln10hi; |
| 77 |
y2 = y*log10_2hi; |
| 78 |
val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi; |
| 79 |
|
| 80 |
/* |
| 81 |
* Extra precision in for adding y*log10_2hi is not strictly needed |
| 82 |
* since there is no very large cancellation near x = sqrt(2) or |
| 83 |
* x = 1/sqrt(2), but we do it anyway since it costs little on CPUs |
| 84 |
* with some parallelism and it reduces the error for many args. |
| 85 |
*/ |
| 86 |
w = y2 + val_hi; |
| 87 |
val_lo += (y2 - w) + val_hi; |
| 88 |
val_hi = w; |
| 89 |
|
| 90 |
return val_lo + val_hi; |
| 91 |
} |
| 92 |
|
| 93 |
#if (LDBL_MANT_DIG == 53) |
| 94 |
__weak_reference(log10, log10l); |
| 95 |
#endif |
