4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
29 #pragma weak __log10 = log10
33 * log10(x) = log(x)/log10
35 * Base on Table look-up algorithm with product polynomial
36 * approximation for log(x).
38 * By K.C. Ng, Nov 29, 2004
40 * (a). For x in [1-0.125, 1+0.125], from log.c we have
41 * log(x) = f + ((a1*f^2) *
42 * ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
43 * (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
44 * ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
46 * (i) modify a1 <- a1 / log10
47 * (ii) 1/log10 = 0.4342944819...
48 * = 0.4375 - 0.003205518... (7 bit shift)
49 * Let lgv = 0.4375 - 1/log10, then
50 * lgv = 0.003205518096748172348871081083395...,
51 * (iii) f*0.4375 is exact because f has 3 trailing zero.
52 * (iv) Thus, log10(x) = f*0.4375 - (lgv*f - PPoly)
54 * (b). For 0.09375 <= x < 24
55 * Let j = (ix - 0x3fb80000) >> 15. Look up Y[j], 1/Y[j], and log(Y[j])
56 * from _TBL_log.c. Then
57 * log10(x) = log10(Y[j]) + log10(1 + (x-Y[j])*(1/Y[j]))
58 * = log(Y[j])(1/log10) + log10(1 + s)
60 * s = (x-Y[j])*(1/Y[j])
61 * From log.c, we have log(1+s) =
63 * (b s) (b + b s + s ) [b + b s + s (b + s)] (b + b s + s )
66 * By setting b1 <- b1/log10, we have
67 * log10(x) = 0.4375 * T - (lgv * T - POLY(s))
69 * (c). Otherwise, get "n", the exponent of x, and then normalize x to
70 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
71 * significant bits. Then
72 * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
73 * log10(x) = n*(ln2/ln10) + log10(z).
76 * log10(x) is NaN with signal if x < 0 (including -INF) ;
77 * log10(+INF) is +INF; log10(0) is -INF with signal;
78 * log10(NaN) is that NaN with no signal.
80 * Maximum error observed: less than 0.89 ulp
83 * The hexadecimal values are the intended ones for the following constants.
84 * The decimal values may be used, provided that the compiler will convert
85 * from decimal to binary accurately enough to produce the hexadecimal values
92 extern const double _TBL_log
[];
94 static const double P
[] = {
96 /* TWO52 */ 4503599627370496.0,
97 /* LNAHI */ 3.01029995607677847147e-01, /* 3FD34413 50900000 */
98 /* LNALO */ 5.63033480667509769841e-11, /* 3DCEF3FD E623E256 */
99 /* A1 */ -2.9142521960136582507385480707044582802184e-02,
100 /* A2 */ 1.99628461483039965074226529395673424005508422852e+0000,
101 /* A3 */ 2.26812367662950720159642514772713184356689453125e+0000,
102 /* A4 */ -9.05030639084976384900471657601883634924888610840e-0001,
103 /* A5 */ -1.48275767132434044270894446526654064655303955078e+0000,
104 /* A6 */ 1.88158320939722756293122074566781520843505859375e+0000,
105 /* A7 */ 1.83309386046986411145098827546462416648864746094e+0000,
106 /* A8 */ 1.24847063988317086291601754055591300129890441895e+0000,
107 /* A9 */ 1.98372421445537705508854742220137268304824829102e+0000,
108 /* A10 */ -3.94711735767898475035764249696512706577777862549e-0001,
109 /* A11 */ 3.07890395362954372160402272129431366920471191406e+0000,
110 /* A12 */ -9.60099585275022149311041630426188930869102478027e-0001,
111 /* B1 */ -5.4304894950350052960838096752491540286689e-02,
112 /* B2 */ 1.87161713283355151891381127914642725337613123482e+0000,
113 /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000,
114 /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000,
115 /* B5 */ 1.64822828085258366037635369139024987816810607910e+0000,
116 /* B6 */ -1.24409107065868340669112512841820716857910156250e+0000,
117 /* B7 */ 1.70534231658220414296067701798165217041969299316e+0000,
118 /* B8 */ 1.99196833784655646937267192697618156671524047852e+0000,
120 /* LGL */ 0.003205518096748172348871081083395,
121 /* LG10V */ 0.43429448190325182765112891891660509576226,
154 double *tb
, dn
, dn1
, s
, z
, r
, w
;
155 int i
, hx
, ix
, n
, lx
;
158 hx
= ((int *)&x
)[HIWORD
];
159 ix
= hx
& 0x7fffffff;
160 lx
= ((int *)&x
)[LOWORD
];
162 /* subnormal,0,negative,inf,nan */
163 if ((hx
+ 0x100000) < 0x200000) {
164 if (ix
> 0x7ff00000 || (ix
== 0x7ff00000 && lx
!= 0)) /* nan */
166 if (((hx
<< 1) | lx
) == 0) /* zero */
167 return (_SVID_libm_err(x
, x
, 18));
168 if (hx
< 0) /* negative */
169 return (_SVID_libm_err(x
, x
, 19));
170 if (((hx
- 0x7ff00000) | lx
) == 0) /* +inf */
173 /* x must be positive and subnormal */
176 ix
= ((int *)&x
)[HIWORD
];
177 lx
= ((int *)&x
)[LOWORD
];
181 if (i
>= 0x7f7 && i
<= 0x806) {
182 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
183 if (ix
>= 0x3fec0000 && ix
< 0x3ff20000) {
184 /* 0.875 <= x < 1.125 */
187 if (((ix
- 0x3ff00000) | lx
) == 0) /* x = 1 */
189 r
= (A10
* s
) * (A11
+ s
);
191 return (LGH
* s
- (LGL
* s
- ((A1
* z
) *
192 ((A2
+ (A3
* s
) * (A4
+ s
)) + w
* (A5
+ s
))) *
193 (((A6
+ s
* (A7
+ s
)) + w
* (A8
+ s
)) *
194 ((A9
+ r
) + w
* (A12
+ s
)))));
196 i
= (ix
- 0x3fb80000) >> 15;
197 tb
= (double *)_TBL_log
+ (i
+ i
+ i
);
198 s
= (x
- tb
[0]) * tb
[1];
199 return (LGH
* tb
[2] - (LGL
* tb
[2] - ((B1
* s
) *
200 (B2
+ s
* (B3
+ s
))) *
201 (((B4
+ s
* B5
) + (s
* s
) * (B6
+ s
)) *
202 (B7
+ s
* (B8
+ s
)))));
205 dn
= (double)(n
+ ((ix
>> 20) - 0x3ff));
207 i
= (ix
& 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
208 ((int *)&x
)[HIWORD
] = i
;
209 i
= (i
- 0x3fb80000) >> 15;
210 tb
= (double *)_TBL_log
+ (i
+ i
+ i
);
211 s
= (x
- tb
[0]) * tb
[1];
212 dn
= dn
* LNALO
+ tb
[2] * LG10V
;
213 return (dn1
+ (dn
+ ((B1
* s
) *
214 (B2
+ s
* (B3
+ s
))) *
215 (((B4
+ s
* B5
) + (s
* s
) * (B6
+ s
)) *
216 (B7
+ s
* (B8
+ s
)))));