1*05a0b428SJohn Marino /* @(#)e_log10.c 5.1 93/09/24 */
2*05a0b428SJohn Marino /*
3*05a0b428SJohn Marino * ====================================================
4*05a0b428SJohn Marino * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5*05a0b428SJohn Marino *
6*05a0b428SJohn Marino * Developed at SunPro, a Sun Microsystems, Inc. business.
7*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this
8*05a0b428SJohn Marino * software is freely granted, provided that this notice
9*05a0b428SJohn Marino * is preserved.
10*05a0b428SJohn Marino * ====================================================
11*05a0b428SJohn Marino */
12*05a0b428SJohn Marino
13*05a0b428SJohn Marino /* log10(x)
14*05a0b428SJohn Marino * Return the base 10 logarithm of x
15*05a0b428SJohn Marino *
16*05a0b428SJohn Marino * Method :
17*05a0b428SJohn Marino * Let log10_2hi = leading 40 bits of log10(2) and
18*05a0b428SJohn Marino * log10_2lo = log10(2) - log10_2hi,
19*05a0b428SJohn Marino * ivln10 = 1/log(10) rounded.
20*05a0b428SJohn Marino * Then
21*05a0b428SJohn Marino * n = ilogb(x),
22*05a0b428SJohn Marino * if(n<0) n = n+1;
23*05a0b428SJohn Marino * x = scalbn(x,-n);
24*05a0b428SJohn Marino * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
25*05a0b428SJohn Marino *
26*05a0b428SJohn Marino * Note 1:
27*05a0b428SJohn Marino * To guarantee log10(10**n)=n, where 10**n is normal, the rounding
28*05a0b428SJohn Marino * mode must set to Round-to-Nearest.
29*05a0b428SJohn Marino * Note 2:
30*05a0b428SJohn Marino * [1/log(10)] rounded to 53 bits has error .198 ulps;
31*05a0b428SJohn Marino * log10 is monotonic at all binary break points.
32*05a0b428SJohn Marino *
33*05a0b428SJohn Marino * Special cases:
34*05a0b428SJohn Marino * log10(x) is NaN with signal if x < 0;
35*05a0b428SJohn Marino * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
36*05a0b428SJohn Marino * log10(NaN) is that NaN with no signal;
37*05a0b428SJohn Marino * log10(10**N) = N for N=0,1,...,22.
38*05a0b428SJohn Marino *
39*05a0b428SJohn Marino * Constants:
40*05a0b428SJohn Marino * The hexadecimal values are the intended ones for the following constants.
41*05a0b428SJohn Marino * The decimal values may be used, provided that the compiler will convert
42*05a0b428SJohn Marino * from decimal to binary accurately enough to produce the hexadecimal values
43*05a0b428SJohn Marino * shown.
44*05a0b428SJohn Marino */
45*05a0b428SJohn Marino
46*05a0b428SJohn Marino #include <float.h>
47*05a0b428SJohn Marino #include <math.h>
48*05a0b428SJohn Marino
49*05a0b428SJohn Marino #include "math_private.h"
50*05a0b428SJohn Marino
51*05a0b428SJohn Marino static const double
52*05a0b428SJohn Marino two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
53*05a0b428SJohn Marino ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */
54*05a0b428SJohn Marino log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
55*05a0b428SJohn Marino log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
56*05a0b428SJohn Marino
57*05a0b428SJohn Marino static const double zero = 0.0;
58*05a0b428SJohn Marino
59*05a0b428SJohn Marino double
log10(double x)60*05a0b428SJohn Marino log10(double x)
61*05a0b428SJohn Marino {
62*05a0b428SJohn Marino double y,z;
63*05a0b428SJohn Marino int32_t i,k,hx;
64*05a0b428SJohn Marino u_int32_t lx;
65*05a0b428SJohn Marino
66*05a0b428SJohn Marino EXTRACT_WORDS(hx,lx,x);
67*05a0b428SJohn Marino
68*05a0b428SJohn Marino k=0;
69*05a0b428SJohn Marino if (hx < 0x00100000) { /* x < 2**-1022 */
70*05a0b428SJohn Marino if (((hx&0x7fffffff)|lx)==0)
71*05a0b428SJohn Marino return -two54/zero; /* log(+-0)=-inf */
72*05a0b428SJohn Marino if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
73*05a0b428SJohn Marino k -= 54; x *= two54; /* subnormal number, scale up x */
74*05a0b428SJohn Marino GET_HIGH_WORD(hx,x);
75*05a0b428SJohn Marino }
76*05a0b428SJohn Marino if (hx >= 0x7ff00000) return x+x;
77*05a0b428SJohn Marino k += (hx>>20)-1023;
78*05a0b428SJohn Marino i = ((u_int32_t)k&0x80000000)>>31;
79*05a0b428SJohn Marino hx = (hx&0x000fffff)|((0x3ff-i)<<20);
80*05a0b428SJohn Marino y = (double)(k+i);
81*05a0b428SJohn Marino SET_HIGH_WORD(x,hx);
82*05a0b428SJohn Marino z = y*log10_2lo + ivln10*log(x);
83*05a0b428SJohn Marino return z+y*log10_2hi;
84*05a0b428SJohn Marino }
85*05a0b428SJohn Marino
86*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG
87*05a0b428SJohn Marino __strong_alias(log10l, log10);
88*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
89