Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/C/log10.c
+++ new/usr/src/lib/libm/common/C/log10.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
|
↓ open down ↓ |
10 lines elided |
↑ open up ↑ |
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 +
21 22 /*
22 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 24 */
25 +
24 26 /*
25 27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 28 * Use is subject to license terms.
27 29 */
28 30
29 31 #pragma weak __log10 = log10
30 32
31 -/* INDENT OFF */
33 +
32 34 /*
33 35 * log10(x) = log(x)/log10
34 36 *
35 37 * Base on Table look-up algorithm with product polynomial
36 38 * approximation for log(x).
37 39 *
38 40 * By K.C. Ng, Nov 29, 2004
39 41 *
40 42 * (a). For x in [1-0.125, 1+0.125], from log.c we have
41 43 * log(x) = f + ((a1*f^2) *
42 44 * ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
43 45 * (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
44 46 * ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
45 47 * where f = x - 1.
46 48 * (i) modify a1 <- a1 / log10
47 49 * (ii) 1/log10 = 0.4342944819...
48 50 * = 0.4375 - 0.003205518... (7 bit shift)
49 51 * Let lgv = 0.4375 - 1/log10, then
50 52 * lgv = 0.003205518096748172348871081083395...,
51 53 * (iii) f*0.4375 is exact because f has 3 trailing zero.
52 54 * (iv) Thus, log10(x) = f*0.4375 - (lgv*f - PPoly)
53 55 *
54 56 * (b). For 0.09375 <= x < 24
55 57 * Let j = (ix - 0x3fb80000) >> 15. Look up Y[j], 1/Y[j], and log(Y[j])
56 58 * from _TBL_log.c. Then
57 59 * log10(x) = log10(Y[j]) + log10(1 + (x-Y[j])*(1/Y[j]))
58 60 * = log(Y[j])(1/log10) + log10(1 + s)
59 61 * where
60 62 * s = (x-Y[j])*(1/Y[j])
61 63 * From log.c, we have log(1+s) =
62 64 * 2 2 2
63 65 * (b s) (b + b s + s ) [b + b s + s (b + s)] (b + b s + s )
64 66 * 1 2 3 4 5 6 7 8
65 67 *
66 68 * By setting b1 <- b1/log10, we have
67 69 * log10(x) = 0.4375 * T - (lgv * T - POLY(s))
68 70 *
69 71 * (c). Otherwise, get "n", the exponent of x, and then normalize x to
70 72 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
71 73 * significant bits. Then
72 74 * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
73 75 * log10(x) = n*(ln2/ln10) + log10(z).
74 76 *
75 77 * Special cases:
76 78 * log10(x) is NaN with signal if x < 0 (including -INF) ;
77 79 * log10(+INF) is +INF; log10(0) is -INF with signal;
|
↓ open down ↓ |
36 lines elided |
↑ open up ↑ |
78 80 * log10(NaN) is that NaN with no signal.
79 81 *
80 82 * Maximum error observed: less than 0.89 ulp
81 83 *
82 84 * Constants:
83 85 * The hexadecimal values are the intended ones for the following constants.
84 86 * The decimal values may be used, provided that the compiler will convert
85 87 * from decimal to binary accurately enough to produce the hexadecimal values
86 88 * shown.
87 89 */
88 -/* INDENT ON */
89 90
90 91 #include "libm.h"
91 92
92 93 extern const double _TBL_log[];
93 94
94 95 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 */
96 +/* ONE */
97 + 1.0,
98 +/* TWO52 */ 4503599627370496.0,
99 +/* LNAHI */ 3.01029995607677847147e-01, /* 3FD34413 50900000 */
100 +/* LNALO */ 5.63033480667509769841e-11, /* 3DCEF3FD E623E256 */
99 101 /* A1 */ -2.9142521960136582507385480707044582802184e-02,
100 -/* A2 */ 1.99628461483039965074226529395673424005508422852e+0000,
101 -/* A3 */ 2.26812367662950720159642514772713184356689453125e+0000,
102 +/* A2 */ 1.99628461483039965074226529395673424005508422852e+0000,
103 +/* A3 */ 2.26812367662950720159642514772713184356689453125e+0000,
102 104 /* A4 */ -9.05030639084976384900471657601883634924888610840e-0001,
103 105 /* A5 */ -1.48275767132434044270894446526654064655303955078e+0000,
104 -/* A6 */ 1.88158320939722756293122074566781520843505859375e+0000,
105 -/* A7 */ 1.83309386046986411145098827546462416648864746094e+0000,
106 -/* A8 */ 1.24847063988317086291601754055591300129890441895e+0000,
107 -/* A9 */ 1.98372421445537705508854742220137268304824829102e+0000,
106 +/* A6 */ 1.88158320939722756293122074566781520843505859375e+0000,
107 +/* A7 */ 1.83309386046986411145098827546462416648864746094e+0000,
108 +/* A8 */ 1.24847063988317086291601754055591300129890441895e+0000,
109 +/* A9 */ 1.98372421445537705508854742220137268304824829102e+0000,
108 110 /* A10 */ -3.94711735767898475035764249696512706577777862549e-0001,
109 -/* A11 */ 3.07890395362954372160402272129431366920471191406e+0000,
111 +/* A11 */ 3.07890395362954372160402272129431366920471191406e+0000,
110 112 /* A12 */ -9.60099585275022149311041630426188930869102478027e-0001,
111 113 /* B1 */ -5.4304894950350052960838096752491540286689e-02,
112 -/* B2 */ 1.87161713283355151891381127914642725337613123482e+0000,
114 +/* B2 */ 1.87161713283355151891381127914642725337613123482e+0000,
113 115 /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000,
114 116 /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000,
115 -/* B5 */ 1.64822828085258366037635369139024987816810607910e+0000,
117 +/* B5 */ 1.64822828085258366037635369139024987816810607910e+0000,
116 118 /* 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,
119 +/* B7 */ 1.70534231658220414296067701798165217041969299316e+0000,
120 +/* B8 */ 1.99196833784655646937267192697618156671524047852e+0000,
121 +/* LGH */ 0.4375,
122 +/* LGL */ 0.003205518096748172348871081083395,
123 +/* LG10V */ 0.43429448190325182765112891891660509576226,
122 124 };
123 125
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]
126 +#define ONE P[0]
127 +#define TWO52 P[1]
128 +#define LNAHI P[2]
129 +#define LNALO P[3]
130 +#define A1 P[4]
131 +#define A2 P[5]
132 +#define A3 P[6]
133 +#define A4 P[7]
134 +#define A5 P[8]
135 +#define A6 P[9]
136 +#define A7 P[10]
137 +#define A8 P[11]
138 +#define A9 P[12]
139 +#define A10 P[13]
140 +#define A11 P[14]
141 +#define A12 P[15]
142 +#define B1 P[16]
143 +#define B2 P[17]
144 +#define B3 P[18]
145 +#define B4 P[19]
146 +#define B5 P[20]
147 +#define B6 P[21]
148 +#define B7 P[22]
149 +#define B8 P[23]
150 +#define LGH P[24]
151 +#define LGL P[25]
152 +#define LG10V P[26]
151 153
152 154 double
153 -log10(double x) {
154 - double *tb, dn, dn1, s, z, r, w;
155 - int i, hx, ix, n, lx;
155 +log10(double x)
156 +{
157 + double *tb, dn, dn1, s, z, r, w;
158 + int i, hx, ix, n, lx;
156 159
157 160 n = 0;
158 161 hx = ((int *)&x)[HIWORD];
159 162 ix = hx & 0x7fffffff;
160 163 lx = ((int *)&x)[LOWORD];
161 164
162 165 /* subnormal,0,negative,inf,nan */
163 166 if ((hx + 0x100000) < 0x200000) {
164 167 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
165 168 return (x * x);
166 - if (((hx << 1) | lx) == 0) /* zero */
169 +
170 + if (((hx << 1) | lx) == 0) /* zero */
167 171 return (_SVID_libm_err(x, x, 18));
168 - if (hx < 0) /* negative */
172 +
173 + if (hx < 0) /* negative */
169 174 return (_SVID_libm_err(x, x, 19));
170 - if (((hx - 0x7ff00000) | lx) == 0) /* +inf */
175 +
176 + if (((hx - 0x7ff00000) | lx) == 0) /* +inf */
171 177 return (x);
172 178
173 179 /* x must be positive and subnormal */
174 180 x *= TWO52;
175 181 n = -52;
176 182 ix = ((int *)&x)[HIWORD];
177 183 lx = ((int *)&x)[LOWORD];
178 184 }
179 185
180 186 i = ix >> 19;
187 +
181 188 if (i >= 0x7f7 && i <= 0x806) {
182 189 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
183 190 if (ix >= 0x3fec0000 && ix < 0x3ff20000) {
184 191 /* 0.875 <= x < 1.125 */
185 192 s = x - ONE;
186 193 z = s * s;
187 - if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
194 +
195 + if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
188 196 return (z);
197 +
189 198 r = (A10 * s) * (A11 + s);
190 199 w = z * s;
191 200 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)))));
201 + ((A2 + (A3 * s) * (A4 + s)) + w * (A5 + s))) *
202 + (((A6 + s * (A7 + s)) + w * (A8 + s)) *
203 + ((A9 + r) + w * (A12 + s)))));
195 204 } else {
196 205 i = (ix - 0x3fb80000) >> 15;
197 206 tb = (double *)_TBL_log + (i + i + i);
198 207 s = (x - tb[0]) * tb[1];
199 208 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)))));
209 + (B2 + s * (B3 + s))) * (((B4 + s * B5) + (s * s) *
210 + (B6 + s)) * (B7 + s * (B8 + s)))));
203 211 }
204 212 } else {
205 213 dn = (double)(n + ((ix >> 20) - 0x3ff));
206 214 dn1 = dn * LNAHI;
207 215 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
208 216 ((int *)&x)[HIWORD] = i;
209 217 i = (i - 0x3fb80000) >> 15;
210 218 tb = (double *)_TBL_log + (i + i + i);
211 219 s = (x - tb[0]) * tb[1];
212 220 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)))));
221 + return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *
222 + (((B4 + s * B5) + (s * s) * (B6 + s)) *
223 + (B7 + s * (B8 + s)))));
217 224 }
218 225 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX