Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/Q/powl.c
+++ new/usr/src/lib/libm/common/Q/powl.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 #pragma weak __powl = powl
31 32
32 33 #include "libm.h"
33 -#include "xpg6.h" /* __xpg6 */
34 -#define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int
34 +#include "xpg6.h" /* __xpg6 */
35 +#define _C99SUSv3_pow _C99SUSv3_pow_treats_Inf_as_an_even_int
35 36
36 37 #if defined(__sparc)
37 -#define i0 0
38 -#define i1 1
39 -#define i2 2
40 -#define i3 3
38 +#define i0 0
39 +#define i1 1
40 +#define i2 2
41 +#define i3 3
41 42
42 43 static const long double zero = 0.0L, one = 1.0L, two = 2.0L;
43 -
44 44 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
45 -
46 -static const long double
47 - two113 = 10384593717069655257060992658440192.0L,
45 +static const long double two113 = 10384593717069655257060992658440192.0L,
48 46 ln2hi = 6.931471805599453094172319547495844850203e-0001L,
49 47 ln2lo = 1.667085920830552208890449330400379754169e-0025L,
50 48 A2 = 6.666666666666666666666666666666091393804e-0001L,
51 49 A3 = 4.000000000000000000000000407167070220671e-0001L,
52 50 A4 = 2.857142857142857142730077490612903681164e-0001L,
53 51 A5 = 2.222222222222242577702836920812882605099e-0001L,
54 52 A6 = 1.818181816435493395985912667105885828356e-0001L,
55 53 A7 = 1.538537835211839751112067512805496931725e-0001L,
56 54 B1 = 6.666666666666666666666666666666666667787e-0001L,
57 55 B2 = 3.999999999999999999999999999999848524411e-0001L,
58 56 B3 = 2.857142857142857142857142865084581075070e-0001L,
59 57 B4 = 2.222222222222222222222010781800643808497e-0001L,
60 58 B5 = 1.818181818181818185051442171337036403674e-0001L,
61 59 B6 = 1.538461538461508363540720286292008207673e-0001L,
62 60 B7 = 1.333333333506731842033180638329317108428e-0001L,
63 61 B8 = 1.176469984587418890634302788283946761670e-0001L,
64 62 B9 = 1.053794891561452331722969901564862497132e-0001L;
65 63
66 64 static long double
67 -logl_x(long double x, long double *w) {
65 +logl_x(long double x, long double *w)
66 +{
68 67 long double f, f1, v, s, z, qn, h, t;
69 - int *px = (int *) &x;
70 - int *pz = (int *) &z;
68 + int *px = (int *)&x;
69 + int *pz = (int *)&z;
71 70 int i, j, ix, n;
72 71
73 72 n = 0;
74 73 ix = px[i0];
74 +
75 75 if (ix > 0x3ffef03f && ix < 0x3fff0820) { /* 65/63 > x > 63/65 */
76 76 f = x - one;
77 77 z = f * f;
78 +
78 79 if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) {
79 80 *w = zero;
80 81 return (zero); /* log(1)= +0 */
81 82 }
83 +
82 84 qn = one / (two + f);
83 - s = f * qn; /* |s|<2**-6 */
85 + s = f * qn; /* |s|<2**-6 */
84 86 v = s * s;
85 - h = (long double) (2.0 * (double) s);
86 - f1 = (long double) ((double) f);
87 - t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn +
88 - s * (v * (B1 + v * (B2 + v * (B3 + v * (B4 +
89 - v * (B5 + v * (B6 + v * (B7 + v * (B8 + v * B9)))))))));
90 - s = (long double) ((double) (h + t));
87 + h = (long double)(2.0 * (double)s);
88 + f1 = (long double)((double)f);
89 + t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn + s * (v *
90 + (B1 + v * (B2 + v * (B3 + v * (B4 + v * (B5 + v * (B6 + v *
91 + (B7 + v * (B8 + v * B9)))))))));
92 + s = (long double)((double)(h + t));
91 93 *w = t - (s - h);
92 94 return (s);
93 95 }
94 - if (ix < 0x00010000) { /* subnormal x */
96 +
97 + if (ix < 0x00010000) { /* subnormal x */
95 98 x *= two113;
96 99 n = -113;
97 100 ix = px[i0];
98 101 }
102 +
99 103 /* LARGE_N */
100 104 n += ((ix + 0x200) >> 16) - 0x3fff;
101 105 ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */
102 106 px[i0] = ix;
103 107 i = ix + 0x200;
104 108 pz[i0] = i & 0xfffffc00;
105 109 pz[i1] = pz[i2] = pz[i3] = 0;
106 110 qn = one / (x + z);
107 111 f = x - z;
108 112 s = f * qn;
109 - f1 = (long double) ((double) f);
110 - h = (long double) (2.0 * (double) s);
113 + f1 = (long double)((double)f);
114 + h = (long double)(2.0 * (double)s);
111 115 t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1));
112 116 j = (i >> 10) & 0x3f;
113 117 v = s * s;
114 - qn = (long double) n;
118 + qn = (long double)n;
115 119 t += qn * ln2lo + _TBL_logl_lo[j];
116 - t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 +
117 - v * A7))))));
120 + t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 + v *
121 + A7))))));
118 122 v = qn * ln2hi + _TBL_logl_hi[j];
119 123 s = h + v;
120 124 t += (h - (s - v));
121 - z = (long double) ((double) (s + t));
125 + z = (long double)((double)(s + t));
122 126 *w = t - (z - s);
123 127 return (z);
124 128 }
125 129
126 130 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
127 131 static const long double
128 132 invln2_32 = 4.616624130844682903551758979206054839765e+1L,
129 133 ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
130 134 ln2_32lo = 5.209643502595475652782654157501186731779e-27L,
131 135 ln2_64 = 1.083042469624914545964425189778400898568e-2L;
132 136
133 137 long double
134 -powl(long double x, long double y) {
138 +powl(long double x, long double y)
139 +{
135 140 long double z, ax;
136 141 long double y1, y2, w1, w2;
137 142 int sbx, sby, j, k, yisint, m;
138 143 int hx, lx, hy, ly, ahx, ahy;
139 - int *pz = (int *) &z;
140 - int *px = (int *) &x;
141 - int *py = (int *) &y;
144 + int *pz = (int *)&z;
145 + int *px = (int *)&x;
146 + int *py = (int *)&y;
142 147
143 148 hx = px[i0];
144 149 lx = px[i1] | px[i2] | px[i3];
145 150 hy = py[i0];
146 151 ly = py[i1] | py[i2] | py[i3];
147 152 ahx = hx & ~0x80000000;
148 153 ahy = hy & ~0x80000000;
149 154
150 155 if ((ahy | ly) == 0)
151 156 return (one); /* x**+-0 = 1 */
152 - else if (hx == 0x3fff0000 && lx == 0 &&
153 - (__xpg6 & _C99SUSv3_pow) != 0)
157 + else if (hx == 0x3fff0000 && lx == 0 && (__xpg6 & _C99SUSv3_pow) != 0)
154 158 return (one); /* C99: 1**anything = 1 */
155 - else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) ||
156 - ahy > 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0))
159 + else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) || ahy >
160 + 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0))
157 161 return (x + y); /* +-NaN return x+y */
158 - /* includes Sun: 1**NaN = NaN */
159 - sbx = (unsigned) hx >> 31;
160 - sby = (unsigned) hy >> 31;
162 +
163 + /* includes Sun: 1**NaN = NaN */
164 + sbx = (unsigned)hx >> 31;
165 + sby = (unsigned)hy >> 31;
161 166 ax = fabsl(x);
167 +
162 168 /*
163 169 * determine if y is an odd int when x < 0
164 170 * yisint = 0 ... y is not an integer
165 171 * yisint = 1 ... y is an odd int
166 172 * yisint = 2 ... y is an even int
167 173 */
168 174 yisint = 0;
175 +
169 176 if (sbx) {
170 - if (ahy >= 0x40700000) /* if |y|>=2**113 */
171 - yisint = 2; /* even integer y */
172 - else if (ahy >= 0x3fff0000) {
177 + if (ahy >= 0x40700000) { /* if |y|>=2**113 */
178 + yisint = 2; /* even integer y */
179 + } else if (ahy >= 0x3fff0000) {
173 180 k = (ahy >> 16) - 0x3fff; /* exponent */
181 +
174 182 if (k > 80) {
175 - j = ((unsigned) py[i3]) >> (112 - k);
183 + j = ((unsigned)py[i3]) >> (112 - k);
184 +
176 185 if ((j << (112 - k)) == py[i3])
177 186 yisint = 2 - (j & 1);
178 187 } else if (k > 48) {
179 - j = ((unsigned) py[i2]) >> (80 - k);
188 + j = ((unsigned)py[i2]) >> (80 - k);
189 +
180 190 if ((j << (80 - k)) == py[i2])
181 191 yisint = 2 - (j & 1);
182 192 } else if (k > 16) {
183 - j = ((unsigned) py[i1]) >> (48 - k);
193 + j = ((unsigned)py[i1]) >> (48 - k);
194 +
184 195 if ((j << (48 - k)) == py[i1])
185 196 yisint = 2 - (j & 1);
186 197 } else if (ly == 0) {
187 198 j = ahy >> (16 - k);
199 +
188 200 if ((j << (16 - k)) == ahy)
189 201 yisint = 2 - (j & 1);
190 202 }
191 203 }
192 204 }
193 205
194 206 /* special value of y */
195 207 if (ly == 0) {
196 208 if (ahy == 0x7fff0000) { /* y is +-inf */
197 209 if (((ahx - 0x3fff0000) | lx) == 0) {
198 210 if ((__xpg6 & _C99SUSv3_pow) != 0)
199 211 return (one);
200 - /* C99: (-1)**+-inf = 1 */
212 + /* C99: (-1)**+-inf = 1 */
201 213 else
202 214 return (y - y);
203 - /* Sun: (+-1)**+-inf = NaN */
204 - } else if (ahx >= 0x3fff0000)
205 - /* (|x|>1)**+,-inf = inf,0 */
215 +
216 + /* Sun: (+-1)**+-inf = NaN */
217 + } else if (ahx >= 0x3fff0000) {
218 + /* (|x|>1)**+,-inf = inf,0 */
206 219 return (sby == 0 ? y : zero);
207 - else /* (|x|<1)**-,+inf = inf,0 */
220 + } else { /* (|x|<1)**-,+inf = inf,0 */
208 221 return (sby != 0 ? -y : zero);
222 + }
209 223 } else if (ahy == 0x3fff0000) { /* y is +-1 */
210 224 if (sby != 0)
211 225 return (one / x);
212 226 else
213 227 return (x);
214 - } else if (hy == 0x40000000) /* y is 2 */
228 + } else if (hy == 0x40000000) { /* y is 2 */
215 229 return (x * x);
216 - else if (hy == 0x3ffe0000) { /* y is 0.5 */
230 + } else if (hy == 0x3ffe0000) { /* y is 0.5 */
217 231 if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) ==
218 - 0))
232 + 0))
219 233 return (sqrtl(x));
220 234 }
221 235 }
222 236
223 237 /* special value of x */
224 238 if (lx == 0) {
225 239 if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) {
226 - /* x is +-0,+-inf,+-1 */
240 + /* x is +-0,+-inf,+-1 */
227 241 z = ax;
242 +
228 243 if (sby == 1)
229 244 z = one / z; /* z = 1/|x| if y is negative */
245 +
230 246 if (sbx == 1) {
231 247 if (ahx == 0x3fff0000 && yisint == 0)
232 248 z = zero / zero;
233 - /* (-1)**non-int is NaN */
249 + /* (-1)**non-int is NaN */
234 250 else if (yisint == 1)
235 251 z = -z; /* (x<0)**odd = -(|x|**odd) */
236 252 }
253 +
237 254 return (z);
238 255 }
239 256 }
240 257
241 258 /* (x<0)**(non-int) is NaN */
242 259 if (sbx == 1 && yisint == 0)
243 260 return (zero / zero); /* should be volatile */
244 261
245 - /* Now ax is finite, y is finite */
246 - /* first compute log(ax) = w1+w2, with 53 bits w1 */
262 + /*
263 + * Now ax is finite, y is finite
264 + * first compute log(ax) = w1+w2, with 53 bits w1
265 + */
247 266 w1 = logl_x(ax, &w2);
248 267
249 268 /* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
250 269 if (ly == 0 || ahy >= 0x43fe0000) {
251 270 y1 = y * w1;
252 271 y2 = y * w2;
253 272 } else {
254 - y1 = (long double) ((double) y);
273 + y1 = (long double)((double)y);
255 274 y2 = (y - y1) * w1 + y * w2;
256 275 y1 *= w1;
257 276 }
277 +
258 278 z = y1 + y2;
259 279 j = pz[i0];
260 - if ((unsigned) j >= 0xffff0000) { /* NaN or -inf */
280 +
281 + if ((unsigned)j >= 0xffff0000) { /* NaN or -inf */
261 282 if (sbx == 1 && yisint == 1)
262 283 return (one / z);
263 284 else
264 285 return (-one / z);
265 286 } else if ((j & ~0x80000000) < 0x3fc30000) { /* |x|<2^-60 */
266 287 if (sbx == 1 && yisint == 1)
267 288 return (-one - z);
268 289 else
269 290 return (one + z);
270 291 } else if (j > 0) {
271 292 if (j > 0x400d0000) {
272 293 if (sbx == 1 && yisint == 1)
273 294 return (scalbnl(-one, 20000));
274 295 else
275 296 return (scalbnl(one, 20000));
276 297 }
277 - k = (int) (invln2_32 * (z + ln2_64));
298 +
299 + k = (int)(invln2_32 * (z + ln2_64));
278 300 } else {
279 - if ((unsigned) j > 0xc00d0000) {
301 + if ((unsigned)j > 0xc00d0000) {
280 302 if (sbx == 1 && yisint == 1)
281 303 return (scalbnl(-one, -20000));
282 304 else
283 305 return (scalbnl(one, -20000));
284 306 }
285 - k = (int) (invln2_32 * (z - ln2_64));
307 +
308 + k = (int)(invln2_32 * (z - ln2_64));
286 309 }
310 +
287 311 j = k & 0x1f;
288 312 m = k >> 5;
289 313 {
290 314 /* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
291 - long double
292 - t1 = 1.666666666666666666666666666660876387437e-1L,
293 - t2 = -2.777777777777777777777707812093173478756e-3L,
294 - t3 = 6.613756613756613482074280932874221202424e-5L,
295 - t4 = -1.653439153392139954169609822742235851120e-6L,
296 - t5 = 4.175314851769539751387852116610973796053e-8L;
297 - long double t = (long double) k;
315 + long double t1 = 1.666666666666666666666666666660876387437e-1L,
316 + t2 = -2.777777777777777777777707812093173478756e-3L,
317 + t3 = 6.613756613756613482074280932874221202424e-5L,
318 + t4 = -1.653439153392139954169609822742235851120e-6L,
319 + t5 = 4.175314851769539751387852116610973796053e-8L;
320 + long double t = (long double)k;
298 321
299 322 w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo;
300 323 t = w1 * w1;
301 324 w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) -
302 - two;
325 + two;
303 326 z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 -
304 - _TBL_expl_lo[j]);
327 + _TBL_expl_lo[j]);
305 328 }
306 329 j = m + (pz[i0] >> 16);
307 - if (j && (unsigned) j < 0x7fff)
330 +
331 + if (j && (unsigned)j < 0x7fff)
308 332 pz[i0] += m << 16;
309 333 else
310 334 z = scalbnl(z, m);
311 335
312 336 if (sbx == 1 && yisint == 1)
313 - z = -z; /* (-ve)**(odd int) */
337 + z = -z; /* (-ve)**(odd int) */
338 +
314 339 return (z);
315 340 }
316 341 #else
317 342 #error Unsupported Architecture
318 -#endif /* defined(__sparc) */
343 +#endif /* defined(__sparc) */
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX