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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 #pragma weak __logl = logl
31
32 /*
33 * logl(x)
34 * Table look-up algorithm
35 * By K.C. Ng, March 6, 1989
36 *
37 * (a). For x in [31/33,33/31], using a special approximation:
38 * f = x - 1;
39 * s = f/(2.0+f); ... here |s| <= 0.03125
40 * z = s*s;
41 * return f-s*(f-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
42 *
43 * (b). Otherwise, normalize x = 2^n * 1.f.
44 * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
61 * For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
62 * so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
63 * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
64 *
65 *
66 * Special cases:
67 * log(x) is NaN with signal if x < 0 (including -INF) ;
68 * log(+INF) is +INF; log(0) is -INF with signal;
69 * log(NaN) is that NaN with no signal.
70 *
71 * Constants:
72 * The hexadecimal values are the intended ones for the following constants.
73 * The decimal values may be used, provided that the compiler will convert
74 * from decimal to binary accurately enough to produce the hexadecimal values
75 * shown.
76 */
77
78 #include "libm.h"
79
80 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
81
82 static const long double
83 zero = 0.0L,
84 one = 1.0L,
85 two = 2.0L,
86 two113 = 10384593717069655257060992658440192.0L,
87 ln2hi = 6.931471805599453094172319547495844850203e-0001L,
88 ln2lo = 1.667085920830552208890449330400379754169e-0025L,
89 A1 = 2.000000000000000000000000000000000000024e+0000L,
90 A2 = 6.666666666666666666666666666666091393804e-0001L,
91 A3 = 4.000000000000000000000000407167070220671e-0001L,
92 A4 = 2.857142857142857142730077490612903681164e-0001L,
93 A5 = 2.222222222222242577702836920812882605099e-0001L,
94 A6 = 1.818181816435493395985912667105885828356e-0001L,
95 A7 = 1.538537835211839751112067512805496931725e-0001L,
96 B1 = 6.666666666666666666666666666666961498329e-0001L,
97 B2 = 3.999999999999999999999999990037655042358e-0001L,
98 B3 = 2.857142857142857142857273426428347457918e-0001L,
99 B4 = 2.222222222222222221353229049747910109566e-0001L,
100 B5 = 1.818181818181821503532559306309070138046e-0001L,
101 B6 = 1.538461538453809210486356084587356788556e-0001L,
102 B7 = 1.333333344463358756121456892645178795480e-0001L,
103 B8 = 1.176460904783899064854645174603360383792e-0001L,
104 B9 = 1.057293869956598995326368602518056990746e-0001L;
105
106 long double
107 logl(long double x) {
108 long double f, s, z, qn, h, t;
109 int *px = (int *) &x;
110 int *pz = (int *) &z;
111 int i, j, ix, i0, i1, n;
112
113 /* get long double precision word ordering */
114 if (*(int *) &one == 0) {
115 i0 = 3;
116 i1 = 0;
117 } else {
118 i0 = 0;
119 i1 = 3;
120 }
121
122 n = 0;
123 ix = px[i0];
124 if (ix > 0x3ffee0f8) { /* if x > 31/33 */
125 if (ix < 0x3fff1084) { /* if x < 33/31 */
126 f = x - one;
127 z = f * f;
128 if (((ix - 0x3fff0000) | px[i1] | px[2] | px[1]) == 0) {
129 return (zero); /* log(1)= +0 */
130 }
131 s = f / (two + f); /* |s|<2**-8 */
132 z = s * s;
133 return (f - s * (f - z * (B1 + z * (B2 + z * (B3 +
134 z * (B4 + z * (B5 + z * (B6 + z * (B7 +
135 z * (B8 + z * B9))))))))));
136 }
137 if (ix >= 0x7fff0000)
138 return (x + x); /* x is +inf or NaN */
139 goto LARGE_N;
140 }
141 if (ix >= 0x00010000)
142 goto LARGE_N;
143 i = ix & 0x7fffffff;
144 if ((i | px[i1] | px[2] | px[1]) == 0) {
145 px[i0] |= 0x80000000;
146 return (one / x); /* log(0.0) = -inf */
147 }
148 if (ix < 0) {
149 if ((unsigned) ix >= 0xffff0000)
150 return (x - x); /* x is -inf or NaN */
151 return (zero / zero); /* log(x<0) is NaN */
152 }
153 /* subnormal x */
154 x *= two113;
155 n = -113;
156 ix = px[i0];
157 LARGE_N:
158 n += ((ix + 0x200) >> 16) - 0x3fff;
159 ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */
160 px[i0] = ix;
161 i = ix + 0x200;
162 pz[i0] = i & 0xfffffc00;
163 pz[i1] = pz[1] = pz[2] = 0;
164 s = (x - z) / (x + z);
165 j = (i >> 10) & 0x3f;
166 z = s * s;
167 qn = (long double) n;
168 t = qn * ln2lo + _TBL_logl_lo[j];
169 h = qn * ln2hi + _TBL_logl_hi[j];
170 f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 +
171 z * (A6 + z * A7))))));
172 return (h + f);
173 }
|
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 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __logl = logl
32
33 /*
34 * logl(x)
35 * Table look-up algorithm
36 * By K.C. Ng, March 6, 1989
37 *
38 * (a). For x in [31/33,33/31], using a special approximation:
39 * f = x - 1;
40 * s = f/(2.0+f); ... here |s| <= 0.03125
41 * z = s*s;
42 * return f-s*(f-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
43 *
44 * (b). Otherwise, normalize x = 2^n * 1.f.
45 * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
62 * For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
63 * so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
64 * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
65 *
66 *
67 * Special cases:
68 * log(x) is NaN with signal if x < 0 (including -INF) ;
69 * log(+INF) is +INF; log(0) is -INF with signal;
70 * log(NaN) is that NaN with no signal.
71 *
72 * Constants:
73 * The hexadecimal values are the intended ones for the following constants.
74 * The decimal values may be used, provided that the compiler will convert
75 * from decimal to binary accurately enough to produce the hexadecimal values
76 * shown.
77 */
78
79 #include "libm.h"
80
81 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
82 static const long double zero = 0.0L,
83 one = 1.0L,
84 two = 2.0L,
85 two113 = 10384593717069655257060992658440192.0L,
86 ln2hi = 6.931471805599453094172319547495844850203e-0001L,
87 ln2lo = 1.667085920830552208890449330400379754169e-0025L,
88 A1 = 2.000000000000000000000000000000000000024e+0000L,
89 A2 = 6.666666666666666666666666666666091393804e-0001L,
90 A3 = 4.000000000000000000000000407167070220671e-0001L,
91 A4 = 2.857142857142857142730077490612903681164e-0001L,
92 A5 = 2.222222222222242577702836920812882605099e-0001L,
93 A6 = 1.818181816435493395985912667105885828356e-0001L,
94 A7 = 1.538537835211839751112067512805496931725e-0001L,
95 B1 = 6.666666666666666666666666666666961498329e-0001L,
96 B2 = 3.999999999999999999999999990037655042358e-0001L,
97 B3 = 2.857142857142857142857273426428347457918e-0001L,
98 B4 = 2.222222222222222221353229049747910109566e-0001L,
99 B5 = 1.818181818181821503532559306309070138046e-0001L,
100 B6 = 1.538461538453809210486356084587356788556e-0001L,
101 B7 = 1.333333344463358756121456892645178795480e-0001L,
102 B8 = 1.176460904783899064854645174603360383792e-0001L,
103 B9 = 1.057293869956598995326368602518056990746e-0001L;
104
105 long double
106 logl(long double x)
107 {
108 long double f, s, z, qn, h, t;
109 int *px = (int *)&x;
110 int *pz = (int *)&z;
111 int i, j, ix, i0, i1, n;
112
113 /* get long double precision word ordering */
114 if (*(int *)&one == 0) {
115 i0 = 3;
116 i1 = 0;
117 } else {
118 i0 = 0;
119 i1 = 3;
120 }
121
122 n = 0;
123 ix = px[i0];
124
125 if (ix > 0x3ffee0f8) { /* if x > 31/33 */
126 if (ix < 0x3fff1084) { /* if x < 33/31 */
127 f = x - one;
128 z = f * f;
129
130 if (((ix - 0x3fff0000) | px[i1] | px[2] | px[1]) == 0)
131 return (zero); /* log(1)= +0 */
132
133 s = f / (two + f); /* |s|<2**-8 */
134 z = s * s;
135 return (f - s * (f - z * (B1 + z * (B2 + z * (B3 + z *
136 (B4 + z * (B5 + z * (B6 + z * (B7 + z * (B8 + z *
137 B9))))))))));
138 }
139
140 if (ix >= 0x7fff0000)
141 return (x + x); /* x is +inf or NaN */
142
143 goto LARGE_N;
144 }
145
146 if (ix >= 0x00010000)
147 goto LARGE_N;
148
149 i = ix & 0x7fffffff;
150
151 if ((i | px[i1] | px[2] | px[1]) == 0) {
152 px[i0] |= 0x80000000;
153 return (one / x); /* log(0.0) = -inf */
154 }
155
156 if (ix < 0) {
157 if ((unsigned)ix >= 0xffff0000)
158 return (x - x); /* x is -inf or NaN */
159
160 return (zero / zero); /* log(x<0) is NaN */
161 }
162
163 /* subnormal x */
164 x *= two113;
165 n = -113;
166 ix = px[i0];
167 LARGE_N:
168 n += ((ix + 0x200) >> 16) - 0x3fff;
169 ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */
170 px[i0] = ix;
171 i = ix + 0x200;
172 pz[i0] = i & 0xfffffc00;
173 pz[i1] = pz[1] = pz[2] = 0;
174 s = (x - z) / (x + z);
175 j = (i >> 10) & 0x3f;
176 z = s * s;
177 qn = (long double)n;
178 t = qn * ln2lo + _TBL_logl_lo[j];
179 h = qn * ln2hi + _TBL_logl_hi[j];
180 f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + z * (A6 + z *
181 A7))))));
182 return (h + f);
183 }
|