Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/Q/log1pl.c
+++ new/usr/src/lib/libm/common/Q/log1pl.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
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.
|
↓ open down ↓ |
14 lines elided |
↑ open up ↑ |
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 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 +
25 26 /*
26 27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 28 * Use is subject to license terms.
28 29 */
29 30
30 31 #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
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
35 36 #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)
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)
40 41 #endif
41 42
42 43 /*
43 44 * log1pl(x)
44 45 * Table look-up algorithm by modifying logl.c
45 46 * By K.C. Ng, July 6, 1995
46 47 *
47 48 * (a). For 1+x in [31/33,33/31], using a special approximation:
48 49 * s = x/(2.0+x); ... here |s| <= 0.03125
49 50 * z = s*s;
50 51 * return x-s*(x-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
51 52 * (i.e., x is in [-2/33,2/31])
52 53 *
53 54 * (b). Otherwise, normalize 1+x = 2^n * 1.f.
54 - * Here we may need a correction term for 1+x rounded.
55 + * Here we may need a correction term for 1+x rounded.
55 56 * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
56 57 * then
57 58 * log(1+x) = n*ln2 + log(1.g) + log(1.f/1.g).
58 59 * Here the leading and trailing values of log(1.g) are obtained from
59 60 * a size-64 table.
60 61 * For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g). Note that
61 62 * 1.f = 2^-n(1+x)
62 63 *
63 64 * then
64 65 * log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +...
65 66 * Note that |s|<2**-8=0.00390625. We use an odd s-polynomial
66 67 * approximation to compute log(1.f/1.g):
|
↓ open down ↓ |
2 lines elided |
↑ open up ↑ |
67 68 * s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7))))))
68 69 * (Precision is 2**-136.91 bits, absolute error)
69 70 *
70 71 * CAUTION:
71 72 * For x>=1, compute 1+x will lost one bit (OK).
72 73 * For x in [-0.5,-1), 1+x is exact.
73 74 * For x in (-0.5,-2/33]U[2/31,1), up to 4 last bits of x will be lost
74 75 * in 1+x. Therefore, to recover the lost bits, one need to compute
75 76 * 1.f-1.g accurately.
76 77 *
77 - * Let hx = HI(x), m = (hx>>16)-0x3fff (=ilogbl(x)), note that
78 + * Let hx = HI(x), m = (hx>>16)-0x3fff (=ilogbl(x)), note that
78 79 * -2/33 = -0.0606...= 2^-5 * 1.939...,
79 80 * 2/31 = 0.09375 = 2^-4 * 1.500...,
80 81 * so for x in (-0.5,-2/33], -5<=m<=-2, n= -1, 1+f=2*(1+x)
81 82 * for x in [2/33,1), -4<=m<=-1, n= 0, f=x
82 83 *
83 84 * In short:
84 - * if x>0, let g: hg= ((hx + (0x200<<(-m)))>>(10-m))<<(10-m)
85 + * if x>0, let g: hg= ((hx + (0x200<<(-m)))>>(10-m))<<(10-m)
85 86 * then 1.f-1.g = x-g
86 87 * if x<0, let g': hg' =((ix-(0x200)<<(-m-1))>>(9-m))<<(9-m)
87 88 * (ix=hx&0x7fffffff)
88 89 * then 1.f-1.g = 2*(g'+x),
89 90 *
90 91 * (c). The final result is computed by
91 92 * (n*ln2_hi+_TBL_logl_hi[j]) +
92 93 * ( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) )
93 94 *
94 95 * Note.
95 96 * For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
96 97 * so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
97 98 * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
98 99 *
99 100 *
100 101 * Special cases:
101 102 * log(x) is NaN with signal if x < 0 (including -INF) ;
102 103 * log(+INF) is +INF; log(0) is -INF with signal;
103 104 * log(NaN) is that NaN with no signal.
104 105 *
105 106 * Constants:
106 107 * The hexadecimal values are the intended ones for the following constants.
|
↓ open down ↓ |
12 lines elided |
↑ open up ↑ |
107 108 * The decimal values may be used, provided that the compiler will convert
108 109 * from decimal to binary accurately enough to produce the hexadecimal values
109 110 * shown.
110 111 */
111 112
112 113 #pragma weak __log1pl = log1pl
113 114
114 115 #include "libm.h"
115 116
116 117 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;
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;
140 139
141 140 long double
142 -log1pl(long double x) {
141 +log1pl(long double x)
142 +{
143 143 long double f, s, z, qn, h, t, y, g;
144 144 int i, j, ix, iy, n, hx, m;
145 145
146 146 hx = H0(x);
147 147 ix = hx & 0x7fffffff;
148 - if (ix < 0x3ffaf07c) { /* |x|<2/33 */
148 +
149 + if (ix < 0x3ffaf07c) { /* |x|<2/33 */
149 150 if (ix <= 0x3f8d0000) { /* x <= 2**-114, return x */
150 - if ((int) x == 0)
151 + if ((int)x == 0)
151 152 return (x);
152 153 }
154 +
153 155 s = x / (two + x); /* |s|<2**-8 */
154 156 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 + return (x - s * (x - z * (B1 + z * (B2 + z * (B3 + z * (B4 + z *
158 + (B5 + z * (B6 + z * (B7 + z * (B8 + z * B9))))))))));
157 159 }
158 - if (ix >= 0x7fff0000) { /* x is +inf or NaN */
160 +
161 + if (ix >= 0x7fff0000) /* x is +inf or NaN */
159 162 return (x + fabsl(x));
160 - }
163 +
161 164 if (hx < 0 && ix >= 0x3fff0000) {
162 165 if (ix > 0x3fff0000 || (H1(x) | H2(x) | H3(x)) != 0)
163 166 x = zero;
167 +
164 168 return (x / zero); /* log1p(x) is NaN if x<-1 */
165 169 /* log1p(-1) is -inf */
166 170 }
171 +
167 172 if (ix >= 0x7ffeffff)
168 - y = x; /* avoid spurious overflow */
173 + y = x; /* avoid spurious overflow */
169 174 else
170 175 y = one + x;
176 +
171 177 iy = H0(y);
172 178 n = ((iy + 0x200) >> 16) - 0x3fff;
173 179 iy = (iy & 0x0000ffff) | 0x3fff0000; /* scale 1+x to [1,2] */
174 180 H0(y) = iy;
175 181 z = zero;
176 182 m = (ix >> 16) - 0x3fff;
183 +
177 184 /* HI(1+x) = (((hx&0xffff)|0x10000)>>(-m))|0x3fff0000 */
178 - if (n == 0) { /* x in [2/33,1) */
185 + if (n == 0) { /* x in [2/33,1) */
179 186 g = zero;
180 187 H0(g) = ((hx + (0x200 << (-m))) >> (10 - m)) << (10 - m);
181 188 t = x - g;
182 189 i = (((((hx & 0xffff) | 0x10000) >> (-m)) | 0x3fff0000) +
183 - 0x200) >> 10;
190 + 0x200) >> 10;
184 191 H0(z) = i << 10;
185 -
186 192 } else if ((1 + n) == 0 && (ix < 0x3ffe0000)) { /* x in (-0.5,-2/33] */
187 193 g = zero;
188 194 H0(g) = ((ix + (0x200 << (-m - 1))) >> (9 - m)) << (9 - m);
189 195 t = g + x;
190 196 t = t + t;
197 +
191 198 /*
192 199 * HI(2*(1+x)) =
193 200 * ((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000
194 201 */
202 +
195 203 /*
196 204 * i =
197 205 * ((((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000)+
198 206 * 0x200)>>10; H0(z)=i<<10;
199 207 */
200 208 z = two * (one - g);
201 209 i = H0(z) >> 10;
202 210 } else {
203 211 i = (iy + 0x200) >> 10;
204 212 H0(z) = i << 10;
205 213 t = y - z;
206 214 }
207 215
208 216 s = t / (y + z);
209 217 j = i & 0x3f;
210 218 z = s * s;
211 - qn = (long double) n;
219 + qn = (long double)n;
212 220 t = qn * ln2lo + _TBL_logl_lo[j];
213 221 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))))));
222 + f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + z * (A6 + z *
223 + A7))))));
216 224 return (h + f);
217 225 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX