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 #ifdef __LITTLE_ENDIAN
31 #define H0(x) *(3 + (int *) &x)
32 #define H1(x) *(2 + (int *) &x)
33 #define H2(x) *(1 + (int *) &x)
34 #define H3(x) *(int *) &x
35 #else
36 #define H0(x) *(int *) &x
37 #define H1(x) *(1 + (int *) &x)
38 #define H2(x) *(2 + (int *) &x)
39 #define H3(x) *(3 + (int *) &x)
40 #endif
41
42 /*
43 * log1pl(x)
44 * Table look-up algorithm by modifying logl.c
45 * By K.C. Ng, July 6, 1995
46 *
47 * (a). For 1+x in [31/33,33/31], using a special approximation:
48 * s = x/(2.0+x); ... here |s| <= 0.03125
49 * z = s*s;
50 * return x-s*(x-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
51 * (i.e., x is in [-2/33,2/31])
52 *
53 * (b). Otherwise, normalize 1+x = 2^n * 1.f.
54 * Here we may need a correction term for 1+x rounded.
55 * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
56 * then
57 * log(1+x) = n*ln2 + log(1.g) + log(1.f/1.g).
58 * Here the leading and trailing values of log(1.g) are obtained from
59 * a size-64 table.
97 * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
98 *
99 *
100 * Special cases:
101 * log(x) is NaN with signal if x < 0 (including -INF) ;
102 * log(+INF) is +INF; log(0) is -INF with signal;
103 * log(NaN) is that NaN with no signal.
104 *
105 * Constants:
106 * The hexadecimal values are the intended ones for the following constants.
107 * The decimal values may be used, provided that the compiler will convert
108 * from decimal to binary accurately enough to produce the hexadecimal values
109 * shown.
110 */
111
112 #pragma weak __log1pl = log1pl
113
114 #include "libm.h"
115
116 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
117
118 static const long double
119 zero = 0.0L,
120 one = 1.0L,
121 two = 2.0L,
122 ln2hi = 6.931471805599453094172319547495844850203e-0001L,
123 ln2lo = 1.667085920830552208890449330400379754169e-0025L,
124 A1 = 2.000000000000000000000000000000000000024e+0000L,
125 A2 = 6.666666666666666666666666666666091393804e-0001L,
126 A3 = 4.000000000000000000000000407167070220671e-0001L,
127 A4 = 2.857142857142857142730077490612903681164e-0001L,
128 A5 = 2.222222222222242577702836920812882605099e-0001L,
129 A6 = 1.818181816435493395985912667105885828356e-0001L,
130 A7 = 1.538537835211839751112067512805496931725e-0001L,
131 B1 = 6.666666666666666666666666666666961498329e-0001L,
132 B2 = 3.999999999999999999999999990037655042358e-0001L,
133 B3 = 2.857142857142857142857273426428347457918e-0001L,
134 B4 = 2.222222222222222221353229049747910109566e-0001L,
135 B5 = 1.818181818181821503532559306309070138046e-0001L,
136 B6 = 1.538461538453809210486356084587356788556e-0001L,
137 B7 = 1.333333344463358756121456892645178795480e-0001L,
138 B8 = 1.176460904783899064854645174603360383792e-0001L,
139 B9 = 1.057293869956598995326368602518056990746e-0001L;
140
141 long double
142 log1pl(long double x) {
143 long double f, s, z, qn, h, t, y, g;
144 int i, j, ix, iy, n, hx, m;
145
146 hx = H0(x);
147 ix = hx & 0x7fffffff;
148 if (ix < 0x3ffaf07c) { /* |x|<2/33 */
149 if (ix <= 0x3f8d0000) { /* x <= 2**-114, return x */
150 if ((int) x == 0)
151 return (x);
152 }
153 s = x / (two + x); /* |s|<2**-8 */
154 z = s * s;
155 return (x - s * (x - z * (B1 + z * (B2 + z * (B3 + z * (B4 +
156 z * (B5 + z * (B6 + z * (B7 + z * (B8 + z * B9))))))))));
157 }
158 if (ix >= 0x7fff0000) { /* x is +inf or NaN */
159 return (x + fabsl(x));
160 }
161 if (hx < 0 && ix >= 0x3fff0000) {
162 if (ix > 0x3fff0000 || (H1(x) | H2(x) | H3(x)) != 0)
163 x = zero;
164 return (x / zero); /* log1p(x) is NaN if x<-1 */
165 /* log1p(-1) is -inf */
166 }
167 if (ix >= 0x7ffeffff)
168 y = x; /* avoid spurious overflow */
169 else
170 y = one + x;
171 iy = H0(y);
172 n = ((iy + 0x200) >> 16) - 0x3fff;
173 iy = (iy & 0x0000ffff) | 0x3fff0000; /* scale 1+x to [1,2] */
174 H0(y) = iy;
175 z = zero;
176 m = (ix >> 16) - 0x3fff;
177 /* HI(1+x) = (((hx&0xffff)|0x10000)>>(-m))|0x3fff0000 */
178 if (n == 0) { /* x in [2/33,1) */
179 g = zero;
180 H0(g) = ((hx + (0x200 << (-m))) >> (10 - m)) << (10 - m);
181 t = x - g;
182 i = (((((hx & 0xffff) | 0x10000) >> (-m)) | 0x3fff0000) +
183 0x200) >> 10;
184 H0(z) = i << 10;
185
186 } else if ((1 + n) == 0 && (ix < 0x3ffe0000)) { /* x in (-0.5,-2/33] */
187 g = zero;
188 H0(g) = ((ix + (0x200 << (-m - 1))) >> (9 - m)) << (9 - m);
189 t = g + x;
190 t = t + t;
191 /*
192 * HI(2*(1+x)) =
193 * ((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000
194 */
195 /*
196 * i =
197 * ((((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000)+
198 * 0x200)>>10; H0(z)=i<<10;
199 */
200 z = two * (one - g);
201 i = H0(z) >> 10;
202 } else {
203 i = (iy + 0x200) >> 10;
204 H0(z) = i << 10;
205 t = y - z;
206 }
207
208 s = t / (y + z);
209 j = i & 0x3f;
210 z = s * s;
211 qn = (long double) n;
212 t = qn * ln2lo + _TBL_logl_lo[j];
213 h = qn * ln2hi + _TBL_logl_hi[j];
214 f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + z * (A6 +
215 z * A7))))));
216 return (h + f);
217 }
|
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 #ifdef __LITTLE_ENDIAN
32 #define H0(x) *(3 + (int *)&x)
33 #define H1(x) *(2 + (int *)&x)
34 #define H2(x) *(1 + (int *)&x)
35 #define H3(x) *(int *)&x
36 #else
37 #define H0(x) *(int *)&x
38 #define H1(x) *(1 + (int *)&x)
39 #define H2(x) *(2 + (int *)&x)
40 #define H3(x) *(3 + (int *)&x)
41 #endif
42
43 /*
44 * log1pl(x)
45 * Table look-up algorithm by modifying logl.c
46 * By K.C. Ng, July 6, 1995
47 *
48 * (a). For 1+x in [31/33,33/31], using a special approximation:
49 * s = x/(2.0+x); ... here |s| <= 0.03125
50 * z = s*s;
51 * return x-s*(x-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
52 * (i.e., x is in [-2/33,2/31])
53 *
54 * (b). Otherwise, normalize 1+x = 2^n * 1.f.
55 * Here we may need a correction term for 1+x rounded.
56 * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
57 * then
58 * log(1+x) = n*ln2 + log(1.g) + log(1.f/1.g).
59 * Here the leading and trailing values of log(1.g) are obtained from
60 * a size-64 table.
98 * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
99 *
100 *
101 * Special cases:
102 * log(x) is NaN with signal if x < 0 (including -INF) ;
103 * log(+INF) is +INF; log(0) is -INF with signal;
104 * log(NaN) is that NaN with no signal.
105 *
106 * Constants:
107 * The hexadecimal values are the intended ones for the following constants.
108 * The decimal values may be used, provided that the compiler will convert
109 * from decimal to binary accurately enough to produce the hexadecimal values
110 * shown.
111 */
112
113 #pragma weak __log1pl = log1pl
114
115 #include "libm.h"
116
117 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
118 static const long double zero = 0.0L,
119 one = 1.0L,
120 two = 2.0L,
121 ln2hi = 6.931471805599453094172319547495844850203e-0001L,
122 ln2lo = 1.667085920830552208890449330400379754169e-0025L,
123 A1 = 2.000000000000000000000000000000000000024e+0000L,
124 A2 = 6.666666666666666666666666666666091393804e-0001L,
125 A3 = 4.000000000000000000000000407167070220671e-0001L,
126 A4 = 2.857142857142857142730077490612903681164e-0001L,
127 A5 = 2.222222222222242577702836920812882605099e-0001L,
128 A6 = 1.818181816435493395985912667105885828356e-0001L,
129 A7 = 1.538537835211839751112067512805496931725e-0001L,
130 B1 = 6.666666666666666666666666666666961498329e-0001L,
131 B2 = 3.999999999999999999999999990037655042358e-0001L,
132 B3 = 2.857142857142857142857273426428347457918e-0001L,
133 B4 = 2.222222222222222221353229049747910109566e-0001L,
134 B5 = 1.818181818181821503532559306309070138046e-0001L,
135 B6 = 1.538461538453809210486356084587356788556e-0001L,
136 B7 = 1.333333344463358756121456892645178795480e-0001L,
137 B8 = 1.176460904783899064854645174603360383792e-0001L,
138 B9 = 1.057293869956598995326368602518056990746e-0001L;
139
140 long double
141 log1pl(long double x)
142 {
143 long double f, s, z, qn, h, t, y, g;
144 int i, j, ix, iy, n, hx, m;
145
146 hx = H0(x);
147 ix = hx & 0x7fffffff;
148
149 if (ix < 0x3ffaf07c) { /* |x|<2/33 */
150 if (ix <= 0x3f8d0000) { /* x <= 2**-114, return x */
151 if ((int)x == 0)
152 return (x);
153 }
154
155 s = x / (two + x); /* |s|<2**-8 */
156 z = s * s;
157 return (x - s * (x - z * (B1 + z * (B2 + z * (B3 + z * (B4 + z *
158 (B5 + z * (B6 + z * (B7 + z * (B8 + z * B9))))))))));
159 }
160
161 if (ix >= 0x7fff0000) /* x is +inf or NaN */
162 return (x + fabsl(x));
163
164 if (hx < 0 && ix >= 0x3fff0000) {
165 if (ix > 0x3fff0000 || (H1(x) | H2(x) | H3(x)) != 0)
166 x = zero;
167
168 return (x / zero); /* log1p(x) is NaN if x<-1 */
169 /* log1p(-1) is -inf */
170 }
171
172 if (ix >= 0x7ffeffff)
173 y = x; /* avoid spurious overflow */
174 else
175 y = one + x;
176
177 iy = H0(y);
178 n = ((iy + 0x200) >> 16) - 0x3fff;
179 iy = (iy & 0x0000ffff) | 0x3fff0000; /* scale 1+x to [1,2] */
180 H0(y) = iy;
181 z = zero;
182 m = (ix >> 16) - 0x3fff;
183
184 /* HI(1+x) = (((hx&0xffff)|0x10000)>>(-m))|0x3fff0000 */
185 if (n == 0) { /* x in [2/33,1) */
186 g = zero;
187 H0(g) = ((hx + (0x200 << (-m))) >> (10 - m)) << (10 - m);
188 t = x - g;
189 i = (((((hx & 0xffff) | 0x10000) >> (-m)) | 0x3fff0000) +
190 0x200) >> 10;
191 H0(z) = i << 10;
192 } else if ((1 + n) == 0 && (ix < 0x3ffe0000)) { /* x in (-0.5,-2/33] */
193 g = zero;
194 H0(g) = ((ix + (0x200 << (-m - 1))) >> (9 - m)) << (9 - m);
195 t = g + x;
196 t = t + t;
197
198 /*
199 * HI(2*(1+x)) =
200 * ((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000
201 */
202
203 /*
204 * i =
205 * ((((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000)+
206 * 0x200)>>10; H0(z)=i<<10;
207 */
208 z = two * (one - g);
209 i = H0(z) >> 10;
210 } else {
211 i = (iy + 0x200) >> 10;
212 H0(z) = i << 10;
213 t = y - z;
214 }
215
216 s = t / (y + z);
217 j = i & 0x3f;
218 z = s * s;
219 qn = (long double)n;
220 t = qn * ln2lo + _TBL_logl_lo[j];
221 h = qn * ln2hi + _TBL_logl_hi[j];
222 f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + z * (A6 + z *
223 A7))))));
224 return (h + f);
225 }
|