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 __log = log
30
31 /* INDENT OFF */
32 /*
33 * log(x)
34 * Table look-up algorithm with product polynomial approximation.
35 * By K.C. Ng, Oct 23, 2004. Updated Oct 18, 2005.
36 *
37 * (a). For x in [1-0.125, 1+0.1328125], using a special approximation:
38 * Let f = x - 1 and z = f*f.
39 * return f + ((a1*z) *
40 * ((a2 + (a3*f)*(a4+f)) + (f*z)*(a5+f))) *
41 * (((a6 + f*(a7+f)) + (f*z)*(a8+f)) *
42 * ((a9 + (a10*f)*(a11+f)) + (f*z)*(a12+f)))
43 * a1 -6.88821452420390473170286327331268694251775741577e-0002,
44 * a2 1.97493380704769294631262255279580131173133850098e+0000,
45 * a3 2.24963218866067560242072431719861924648284912109e+0000,
46 * a4 -9.02975906958474405783476868236903101205825805664e-0001,
47 * a5 -1.47391630715542865104339398385491222143173217773e+0000,
48 * a6 1.86846544648220058704168877738993614912033081055e+0000,
49 * a7 1.82277370459347465292410106485476717352867126465e+0000,
50 * a8 1.25295479915214102994980294170090928673744201660e+0000,
51 * a9 1.96709676945198275177517643896862864494323730469e+0000,
78 * with remez error |log(1+s) - P(s)| <= 2**-63.5
79 *
80 * (c). Otherwise, get "n", the exponent of x, and then normalize x to
81 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
82 * significant bits. Then
83 * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
84 *
85 * Special cases:
86 * log(x) is NaN with signal if x < 0 (including -INF) ;
87 * log(+INF) is +INF; log(0) is -INF with signal;
88 * log(NaN) is that NaN with no signal.
89 *
90 * Maximum error observed: less than 0.90 ulp
91 *
92 * Constants:
93 * The hexadecimal values are the intended ones for the following constants.
94 * The decimal values may be used, provided that the compiler will convert
95 * from decimal to binary accurately enough to produce the hexadecimal values
96 * shown.
97 */
98 /* INDENT ON */
99
100 #include "libm.h"
101
102 extern const double _TBL_log[];
103
104 static const double P[] = {
105 /* ONE */ 1.0,
106 /* TWO52 */ 4503599627370496.0,
107 /* LN2HI */ 6.93147180369123816490e-01, /* 3fe62e42, fee00000 */
108 /* LN2LO */ 1.90821492927058770002e-10, /* 3dea39ef, 35793c76 */
109 /* A1 */ -6.88821452420390473170286327331268694251775741577e-0002,
110 /* A2 */ 1.97493380704769294631262255279580131173133850098e+0000,
111 /* A3 */ 2.24963218866067560242072431719861924648284912109e+0000,
112 /* A4 */ -9.02975906958474405783476868236903101205825805664e-0001,
113 /* A5 */ -1.47391630715542865104339398385491222143173217773e+0000,
114 /* A6 */ 1.86846544648220058704168877738993614912033081055e+0000,
115 /* A7 */ 1.82277370459347465292410106485476717352867126465e+0000,
116 /* A8 */ 1.25295479915214102994980294170090928673744201660e+0000,
117 /* A9 */ 1.96709676945198275177517643896862864494323730469e+0000,
118 /* A10 */ -4.00127989749189894030934055990655906498432159424e-0001,
119 /* A11 */ 3.01675528558798333733648178167641162872314453125e+0000,
120 /* A12 */ -9.52325445049240770778453679668018594384193420410e-0001,
121 /* B1 */ -1.25041641589283658575482149899471551179885864258e-0001,
122 /* B2 */ 1.87161713283355151891381127914642725337613123482e+0000,
123 /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000,
124 /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000,
125 /* B5 */ 1.64822828085258366037635369139024987816810607910e+0000,
137 #define A3 P[6]
138 #define A4 P[7]
139 #define A5 P[8]
140 #define A6 P[9]
141 #define A7 P[10]
142 #define A8 P[11]
143 #define A9 P[12]
144 #define A10 P[13]
145 #define A11 P[14]
146 #define A12 P[15]
147 #define B1 P[16]
148 #define B2 P[17]
149 #define B3 P[18]
150 #define B4 P[19]
151 #define B5 P[20]
152 #define B6 P[21]
153 #define B7 P[22]
154 #define B8 P[23]
155
156 double
157 log(double x) {
158 double *tb, dn, dn1, s, z, r, w;
159 int i, hx, ix, n, lx;
160
161 n = 0;
162 hx = ((int *)&x)[HIWORD];
163 ix = hx & 0x7fffffff;
164 lx = ((int *)&x)[LOWORD];
165
166 /* subnormal,0,negative,inf,nan */
167 if ((hx + 0x100000) < 0x200000) {
168 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
169 return (x * x);
170 if (((hx << 1) | lx) == 0) /* zero */
171 return (_SVID_libm_err(x, x, 16));
172 if (hx < 0) /* negative */
173 return (_SVID_libm_err(x, x, 17));
174 if (((hx - 0x7ff00000) | lx) == 0) /* +inf */
175 return (x);
176
177 /* x must be positive and subnormal */
178 x *= TWO52;
179 n = -52;
180 ix = ((int *)&x)[HIWORD];
181 lx = ((int *)&x)[LOWORD];
182 }
183
184 i = ix >> 19;
185 if (i >= 0x7f7 && i <= 0x806) {
186 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
187 if (ix >= 0x3fec0000 && ix < 0x3ff22000) {
188 /* 0.875 <= x < 1.125 */
189 s = x - ONE;
190 z = s * s;
191 if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
192 return (z);
193 r = (A10 * s) * (A11 + s);
194 w = z * s;
195 return (s + ((A1 * z) *
196 (A2 + ((A3 * s) * (A4 + s) + w * (A5 + s)))) *
197 ((A6 + (s * (A7 + s) + w * (A8 + s))) *
198 (A9 + (r + w * (A12 + s)))));
199 } else {
200 i = (ix - 0x3fb80000) >> 15;
201 tb = (double *)_TBL_log + (i + i + i);
202 s = (x - tb[0]) * tb[1];
203 return (tb[2] + ((B1 * s) * (B2 + s * (B3 + s))) *
204 (((B4 + s * B5) + (s * s) * (B6 + s)) *
205 (B7 + s * (B8 + s))));
206 }
207 } else {
208 dn = (double)(n + ((ix >> 20) - 0x3ff));
209 dn1 = dn * LN2HI;
210 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
211 ((int *)&x)[HIWORD] = i;
212 i = (i - 0x3fb80000) >> 15;
213 tb = (double *)_TBL_log + (i + i + i);
214 s = (x - tb[0]) * tb[1];
215 dn = dn * LN2LO + tb[2];
216 return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *
217 (((B4 + s * B5) + (s * s) * (B6 + s)) *
218 (B7 + s * (B8 + s)))));
|
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 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __log = log
32
33
34 /*
35 * log(x)
36 * Table look-up algorithm with product polynomial approximation.
37 * By K.C. Ng, Oct 23, 2004. Updated Oct 18, 2005.
38 *
39 * (a). For x in [1-0.125, 1+0.1328125], using a special approximation:
40 * Let f = x - 1 and z = f*f.
41 * return f + ((a1*z) *
42 * ((a2 + (a3*f)*(a4+f)) + (f*z)*(a5+f))) *
43 * (((a6 + f*(a7+f)) + (f*z)*(a8+f)) *
44 * ((a9 + (a10*f)*(a11+f)) + (f*z)*(a12+f)))
45 * a1 -6.88821452420390473170286327331268694251775741577e-0002,
46 * a2 1.97493380704769294631262255279580131173133850098e+0000,
47 * a3 2.24963218866067560242072431719861924648284912109e+0000,
48 * a4 -9.02975906958474405783476868236903101205825805664e-0001,
49 * a5 -1.47391630715542865104339398385491222143173217773e+0000,
50 * a6 1.86846544648220058704168877738993614912033081055e+0000,
51 * a7 1.82277370459347465292410106485476717352867126465e+0000,
52 * a8 1.25295479915214102994980294170090928673744201660e+0000,
53 * a9 1.96709676945198275177517643896862864494323730469e+0000,
80 * with remez error |log(1+s) - P(s)| <= 2**-63.5
81 *
82 * (c). Otherwise, get "n", the exponent of x, and then normalize x to
83 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
84 * significant bits. Then
85 * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
86 *
87 * Special cases:
88 * log(x) is NaN with signal if x < 0 (including -INF) ;
89 * log(+INF) is +INF; log(0) is -INF with signal;
90 * log(NaN) is that NaN with no signal.
91 *
92 * Maximum error observed: less than 0.90 ulp
93 *
94 * Constants:
95 * The hexadecimal values are the intended ones for the following constants.
96 * The decimal values may be used, provided that the compiler will convert
97 * from decimal to binary accurately enough to produce the hexadecimal values
98 * shown.
99 */
100
101 #include "libm.h"
102
103 extern const double _TBL_log[];
104
105 static const double P[] = {
106 /* ONE */
107 1.0,
108 /* TWO52 */ 4503599627370496.0,
109 /* LN2HI */ 6.93147180369123816490e-01, /* 3fe62e42, fee00000 */
110 /* LN2LO */ 1.90821492927058770002e-10, /* 3dea39ef, 35793c76 */
111 /* A1 */ -6.88821452420390473170286327331268694251775741577e-0002,
112 /* A2 */ 1.97493380704769294631262255279580131173133850098e+0000,
113 /* A3 */ 2.24963218866067560242072431719861924648284912109e+0000,
114 /* A4 */ -9.02975906958474405783476868236903101205825805664e-0001,
115 /* A5 */ -1.47391630715542865104339398385491222143173217773e+0000,
116 /* A6 */ 1.86846544648220058704168877738993614912033081055e+0000,
117 /* A7 */ 1.82277370459347465292410106485476717352867126465e+0000,
118 /* A8 */ 1.25295479915214102994980294170090928673744201660e+0000,
119 /* A9 */ 1.96709676945198275177517643896862864494323730469e+0000,
120 /* A10 */ -4.00127989749189894030934055990655906498432159424e-0001,
121 /* A11 */ 3.01675528558798333733648178167641162872314453125e+0000,
122 /* A12 */ -9.52325445049240770778453679668018594384193420410e-0001,
123 /* B1 */ -1.25041641589283658575482149899471551179885864258e-0001,
124 /* B2 */ 1.87161713283355151891381127914642725337613123482e+0000,
125 /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000,
126 /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000,
127 /* B5 */ 1.64822828085258366037635369139024987816810607910e+0000,
139 #define A3 P[6]
140 #define A4 P[7]
141 #define A5 P[8]
142 #define A6 P[9]
143 #define A7 P[10]
144 #define A8 P[11]
145 #define A9 P[12]
146 #define A10 P[13]
147 #define A11 P[14]
148 #define A12 P[15]
149 #define B1 P[16]
150 #define B2 P[17]
151 #define B3 P[18]
152 #define B4 P[19]
153 #define B5 P[20]
154 #define B6 P[21]
155 #define B7 P[22]
156 #define B8 P[23]
157
158 double
159 log(double x)
160 {
161 double *tb, dn, dn1, s, z, r, w;
162 int i, hx, ix, n, lx;
163
164 n = 0;
165 hx = ((int *)&x)[HIWORD];
166 ix = hx & 0x7fffffff;
167 lx = ((int *)&x)[LOWORD];
168
169 /* subnormal,0,negative,inf,nan */
170 if ((hx + 0x100000) < 0x200000) {
171 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
172 return (x * x);
173
174 if (((hx << 1) | lx) == 0) /* zero */
175 return (_SVID_libm_err(x, x, 16));
176
177 if (hx < 0) /* negative */
178 return (_SVID_libm_err(x, x, 17));
179
180 if (((hx - 0x7ff00000) | lx) == 0) /* +inf */
181 return (x);
182
183 /* x must be positive and subnormal */
184 x *= TWO52;
185 n = -52;
186 ix = ((int *)&x)[HIWORD];
187 lx = ((int *)&x)[LOWORD];
188 }
189
190 i = ix >> 19;
191
192 if (i >= 0x7f7 && i <= 0x806) {
193 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
194 if (ix >= 0x3fec0000 && ix < 0x3ff22000) {
195 /* 0.875 <= x < 1.125 */
196 s = x - ONE;
197 z = s * s;
198
199 if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
200 return (z);
201
202 r = (A10 * s) * (A11 + s);
203 w = z * s;
204 return (s + ((A1 * z) * (A2 + ((A3 * s) * (A4 + s) + w *
205 (A5 + s)))) * ((A6 + (s * (A7 + s) + w *
206 (A8 + s))) * (A9 + (r + w * (A12 + s)))));
207 } else {
208 i = (ix - 0x3fb80000) >> 15;
209 tb = (double *)_TBL_log + (i + i + i);
210 s = (x - tb[0]) * tb[1];
211 return (tb[2] + ((B1 * s) * (B2 + s * (B3 + s))) *
212 (((B4 + s * B5) + (s * s) * (B6 + s)) *
213 (B7 + s * (B8 + s))));
214 }
215 } else {
216 dn = (double)(n + ((ix >> 20) - 0x3ff));
217 dn1 = dn * LN2HI;
218 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
219 ((int *)&x)[HIWORD] = i;
220 i = (i - 0x3fb80000) >> 15;
221 tb = (double *)_TBL_log + (i + i + i);
222 s = (x - tb[0]) * tb[1];
223 dn = dn * LN2LO + tb[2];
224 return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *
225 (((B4 + s * B5) + (s * s) * (B6 + s)) *
226 (B7 + s * (B8 + s)))));
|