xref: /illumos-gate/usr/src/lib/libm/common/C/log10.c (revision 32991bedc3a6475f1401855c2318ae5b15f8a16b)
1 /*
2  * CDDL HEADER START
3  *
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.
7  *
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.
12  *
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]
18  *
19  * CDDL HEADER END
20  */
21 /*
22  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
23  */
24 /*
25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
26  * Use is subject to license terms.
27  */
28 
29 #pragma weak __log10 = log10
30 
31 /* INDENT OFF */
32 /*
33  * log10(x) = log(x)/log10
34  *
35  * Base on Table look-up algorithm with product polynomial
36  * approximation for log(x).
37  *
38  * By K.C. Ng, Nov 29, 2004
39  *
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)))
45  *	where f = x - 1.
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)
53  *
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)
59  *	where
60  *		s = (x-Y[j])*(1/Y[j])
61  *	From log.c, we have log(1+s) =
62  *				  2              2                     2
63  *		(b s) (b + b s + s ) [b + b s + s (b + s)] (b + b s + s )
64  *		  1     2   3          4   5        6        7   8
65  *
66  *	By setting b1 <- b1/log10, we have
67  *		log10(x) = 0.4375 * T - (lgv * T - POLY(s))
68  *
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).
74  *
75  * Special cases:
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.
79  *
80  * Maximum error observed: less than 0.89 ulp
81  *
82  * Constants:
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
86  * shown.
87  */
88 /* INDENT ON */
89 
90 #include "libm.h"
91 
92 extern const double _TBL_log[];
93 
94 static const double P[] = {
95 /* ONE   */  1.0,
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,
119 /* LGH   */  0.4375,
120 /* LGL   */  0.003205518096748172348871081083395,
121 /* LG10V */  0.43429448190325182765112891891660509576226,
122 };
123 
124 #define	ONE   P[0]
125 #define	TWO52 P[1]
126 #define	LNAHI P[2]
127 #define	LNALO P[3]
128 #define	A1    P[4]
129 #define	A2    P[5]
130 #define	A3    P[6]
131 #define	A4    P[7]
132 #define	A5    P[8]
133 #define	A6    P[9]
134 #define	A7    P[10]
135 #define	A8    P[11]
136 #define	A9    P[12]
137 #define	A10   P[13]
138 #define	A11   P[14]
139 #define	A12   P[15]
140 #define	B1    P[16]
141 #define	B2    P[17]
142 #define	B3    P[18]
143 #define	B4    P[19]
144 #define	B5    P[20]
145 #define	B6    P[21]
146 #define	B7    P[22]
147 #define	B8    P[23]
148 #define	LGH   P[24]
149 #define	LGL   P[25]
150 #define	LG10V P[26]
151 
152 double
153 log10(double x) {
154 	double	*tb, dn, dn1, s, z, r, w;
155 	int	i, hx, ix, n, lx;
156 
157 	n = 0;
158 	hx = ((int *)&x)[HIWORD];
159 	ix = hx & 0x7fffffff;
160 	lx = ((int *)&x)[LOWORD];
161 
162 	/* subnormal,0,negative,inf,nan */
163 	if ((hx + 0x100000) < 0x200000) {
164 		if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
165 			return (x * x);
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 */
171 			return (x);
172 
173 		/* x must be positive and subnormal */
174 		x *= TWO52;
175 		n = -52;
176 		ix = ((int *)&x)[HIWORD];
177 		lx = ((int *)&x)[LOWORD];
178 	}
179 
180 	i = ix >> 19;
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 */
185 			s = x - ONE;
186 			z = s * s;
187 			if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
188 				return (z);
189 			r = (A10 * s) * (A11 + s);
190 			w = z * 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)))));
195 		} else {
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)))));
203 		}
204 	} else {
205 		dn = (double)(n + ((ix >> 20) - 0x3ff));
206 		dn1 = dn * LNAHI;
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)))));
217 	}
218 }
219