Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/Q/sqrtl.c
+++ new/usr/src/lib/libm/common/Q/sqrtl.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 __sqrtl = sqrtl
31 32
32 33 #include "libm.h"
33 34 #include "longdouble.h"
34 35
35 36 extern int __swapTE(int);
36 37 extern int __swapEX(int);
|
↓ open down ↓ |
2 lines elided |
↑ open up ↑ |
37 38 extern enum fp_direction_type __swapRD(enum fp_direction_type);
38 39
39 40 /*
40 41 * in struct longdouble, msw consists of
41 42 * unsigned short sgn:1;
42 43 * unsigned short exp:15;
43 44 * unsigned short frac1:16;
44 45 */
45 46
46 47 #ifdef __LITTLE_ENDIAN
47 -
48 48 /* array indices used to access words within a double */
49 -#define HIWORD 1
50 -#define LOWORD 0
49 +#define HIWORD 1
50 +#define LOWORD 0
51 51
52 52 /* structure used to access words within a quad */
53 53 union longdouble {
54 54 struct {
55 - unsigned int frac4;
56 - unsigned int frac3;
57 - unsigned int frac2;
58 - unsigned int msw;
55 + unsigned int frac4;
56 + unsigned int frac3;
57 + unsigned int frac2;
58 + unsigned int msw;
59 59 } l;
60 - long double d;
60 + long double d;
61 61 };
62 62
63 63 /* default NaN returned for sqrt(neg) */
64 -static const union longdouble
65 - qnan = { 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff };
64 +static const union longdouble qnan = {
65 + 0xffffffff, 0xffffffff, 0xffffffff, 0x7fffffff
66 +};
66 67
67 68 /* signalling NaN used to raise invalid */
68 69 static const union {
69 70 unsigned u[2];
70 71 double d;
71 72 } snan = { 0, 0x7ff00001 };
72 -
73 73 #else
74 -
75 74 /* array indices used to access words within a double */
76 -#define HIWORD 0
77 -#define LOWORD 1
75 +#define HIWORD 0
76 +#define LOWORD 1
78 77
79 78 /* structure used to access words within a quad */
80 79 union longdouble {
81 80 struct {
82 - unsigned int msw;
83 - unsigned int frac2;
84 - unsigned int frac3;
85 - unsigned int frac4;
81 + unsigned int msw;
82 + unsigned int frac2;
83 + unsigned int frac3;
84 + unsigned int frac4;
86 85 } l;
87 - long double d;
86 + long double d;
88 87 };
89 88
90 89 /* default NaN returned for sqrt(neg) */
91 -static const union longdouble
92 - qnan = { 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff };
90 +static const union longdouble qnan = {
91 + 0x7fffffff, 0xffffffff, 0xffffffff, 0xffffffff
92 +};
93 93
94 94 /* signalling NaN used to raise invalid */
95 95 static const union {
96 96 unsigned u[2];
97 97 double d;
98 98 } snan = { 0x7ff00001, 0 };
99 -
100 99 #endif /* __LITTLE_ENDIAN */
101 100
102 -
103 -static const double
104 - zero = 0.0,
101 +static const double zero = 0.0,
105 102 half = 0.5,
106 103 one = 1.0,
107 104 huge = 1.0e300,
108 105 tiny = 1.0e-300,
109 106 two36 = 6.87194767360000000000e+10,
110 107 two30 = 1.07374182400000000000e+09,
111 108 two6 = 6.40000000000000000000e+01,
112 109 two4 = 1.60000000000000000000e+01,
113 110 twom18 = 3.81469726562500000000e-06,
114 111 twom28 = 3.72529029846191406250e-09,
115 112 twom42 = 2.27373675443232059479e-13,
116 113 twom60 = 8.67361737988403547206e-19,
117 114 twom62 = 2.16840434497100886801e-19,
118 115 twom66 = 1.35525271560688054251e-20,
119 116 twom90 = 8.07793566946316088742e-28,
120 117 twom113 = 9.62964972193617926528e-35,
121 118 twom124 = 4.70197740328915003187e-38;
122 119
123 -
124 120 /*
125 -* Extract the exponent and normalized significand (represented as
126 -* an array of five doubles) from a finite, nonzero quad.
127 -*/
121 + * Extract the exponent and normalized significand (represented as
122 + * an array of five doubles) from a finite, nonzero quad.
123 + */
128 124 static int
129 125 __q_unpack(const union longdouble *x, double *s)
130 126 {
131 127 union {
132 - double d;
133 - unsigned int l[2];
128 + double d;
129 + unsigned int l[2];
134 130 } u;
135 - double b;
136 - unsigned int lx, w[3];
137 - int ex;
131 +
132 + double b;
133 + unsigned int lx, w[3];
134 + int ex;
138 135
139 136 /* get the normalized significand and exponent */
140 - ex = (int) ((x->l.msw & 0x7fffffff) >> 16);
137 + ex = (int)((x->l.msw & 0x7fffffff) >> 16);
141 138 lx = x->l.msw & 0xffff;
142 - if (ex)
143 - {
139 +
140 + if (ex) {
144 141 lx |= 0x10000;
145 142 w[0] = x->l.frac2;
146 143 w[1] = x->l.frac3;
147 144 w[2] = x->l.frac4;
148 - }
149 - else
150 - {
151 - if (lx | (x->l.frac2 & 0xfffe0000))
152 - {
145 + } else {
146 + if (lx | (x->l.frac2 & 0xfffe0000)) {
153 147 w[0] = x->l.frac2;
154 148 w[1] = x->l.frac3;
155 149 w[2] = x->l.frac4;
156 150 ex = 1;
157 - }
158 - else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000))
159 - {
151 + } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
160 152 lx = x->l.frac2;
161 153 w[0] = x->l.frac3;
162 154 w[1] = x->l.frac4;
163 155 w[2] = 0;
164 156 ex = -31;
165 - }
166 - else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000))
167 - {
157 + } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
168 158 lx = x->l.frac3;
169 159 w[0] = x->l.frac4;
170 160 w[1] = w[2] = 0;
171 161 ex = -63;
172 - }
173 - else
174 - {
162 + } else {
175 163 lx = x->l.frac4;
176 164 w[0] = w[1] = w[2] = 0;
177 165 ex = -95;
178 166 }
179 - while ((lx & 0x10000) == 0)
180 - {
167 +
168 + while ((lx & 0x10000) == 0) {
181 169 lx = (lx << 1) | (w[0] >> 31);
182 170 w[0] = (w[0] << 1) | (w[1] >> 31);
183 171 w[1] = (w[1] << 1) | (w[2] >> 31);
184 172 w[2] <<= 1;
185 173 ex--;
186 174 }
187 175 }
188 176
189 177 /* extract the significand into five doubles */
190 178 u.l[HIWORD] = 0x42300000;
191 179 u.l[LOWORD] = 0;
192 180 b = u.d;
193 181 u.l[LOWORD] = lx;
194 182 s[0] = u.d - b;
195 183
196 184 u.l[HIWORD] = 0x40300000;
197 185 u.l[LOWORD] = 0;
198 186 b = u.d;
199 187 u.l[LOWORD] = w[0] & 0xffffff00;
200 188 s[1] = u.d - b;
201 189
202 190 u.l[HIWORD] = 0x3e300000;
203 191 u.l[LOWORD] = 0;
204 192 b = u.d;
205 193 u.l[HIWORD] |= w[0] & 0xff;
206 194 u.l[LOWORD] = w[1] & 0xffff0000;
207 195 s[2] = u.d - b;
208 196
209 197 u.l[HIWORD] = 0x3c300000;
210 198 u.l[LOWORD] = 0;
211 199 b = u.d;
|
↓ open down ↓ |
21 lines elided |
↑ open up ↑ |
212 200 u.l[HIWORD] |= w[1] & 0xffff;
213 201 u.l[LOWORD] = w[2] & 0xff000000;
214 202 s[3] = u.d - b;
215 203
216 204 u.l[HIWORD] = 0x3c300000;
217 205 u.l[LOWORD] = 0;
218 206 b = u.d;
219 207 u.l[LOWORD] = w[2] & 0xffffff;
220 208 s[4] = u.d - b;
221 209
222 - return ex - 0x3fff;
210 + return (ex - 0x3fff);
223 211 }
224 212
225 -
226 213 /*
227 -* Pack an exponent and array of three doubles representing a finite,
228 -* nonzero number into a quad. Assume the sign is already there and
229 -* the rounding mode has been fudged accordingly.
230 -*/
214 + * Pack an exponent and array of three doubles representing a finite,
215 + * nonzero number into a quad. Assume the sign is already there and
216 + * the rounding mode has been fudged accordingly.
217 + */
231 218 static void
232 219 __q_pack(const double *z, int exp, enum fp_direction_type rm,
233 - union longdouble *x, int *inexact)
220 + union longdouble *x, int *inexact)
234 221 {
235 222 union {
236 - double d;
237 - unsigned int l[2];
223 + double d;
224 + unsigned int l[2];
238 225 } u;
239 - double s[3], t, t2;
240 - unsigned int msw, frac2, frac3, frac4;
226 +
227 + double s[3], t, t2;
228 + unsigned int msw, frac2, frac3, frac4;
241 229
242 230 /* bias exponent and strip off integer bit */
243 231 exp += 0x3fff;
244 232 s[0] = z[0] - one;
245 233 s[1] = z[1];
246 234 s[2] = z[2];
247 235
248 236 /*
249 237 * chop the significand to obtain the fraction;
250 238 * use round-to-minus-infinity to ensure chopping
251 239 */
252 240 (void) __swapRD(fp_negative);
253 241
254 242 /* extract the first eighty bits of fraction */
255 243 t = s[1] + s[2];
256 244 u.d = two36 + (s[0] + t);
257 245 msw = u.l[LOWORD];
258 246 s[0] -= (u.d - two36);
259 247
260 248 u.d = two4 + (s[0] + t);
261 249 frac2 = u.l[LOWORD];
262 250 s[0] -= (u.d - two4);
263 251
264 252 u.d = twom28 + (s[0] + t);
265 253 frac3 = u.l[LOWORD];
266 254 s[0] -= (u.d - twom28);
267 255
268 256 /* condense the remaining fraction; errors here won't matter */
269 257 t = s[0] + s[1];
270 258 s[1] = ((s[0] - t) + s[1]) + s[2];
271 259 s[0] = t;
272 260
273 261 /* get the last word of fraction */
274 262 u.d = twom60 + (s[0] + s[1]);
|
↓ open down ↓ |
24 lines elided |
↑ open up ↑ |
275 263 frac4 = u.l[LOWORD];
276 264 s[0] -= (u.d - twom60);
277 265
278 266 /*
279 267 * keep track of what's left for rounding; note that
280 268 * t2 will be non-negative due to rounding mode
281 269 */
282 270 t = s[0] + s[1];
283 271 t2 = (s[0] - t) + s[1];
284 272
285 - if (t != zero)
286 - {
273 + if (t != zero) {
287 274 *inexact = 1;
288 275
289 276 /* decide whether to round the fraction up */
290 277 if (rm == fp_positive || (rm == fp_nearest && (t > twom113 ||
291 - (t == twom113 && (t2 != zero || frac4 & 1)))))
292 - {
278 + (t == twom113 && (t2 != zero || frac4 & 1))))) {
293 279 /* round up and renormalize if necessary */
294 - if (++frac4 == 0)
295 - if (++frac3 == 0)
296 - if (++frac2 == 0)
297 - if (++msw == 0x10000)
298 - {
280 + if (++frac4 == 0) {
281 + if (++frac3 == 0) {
282 + if (++frac2 == 0) {
283 + if (++msw == 0x10000) {
299 284 msw = 0;
300 285 exp++;
301 286 }
287 + }
288 + }
289 + }
302 290 }
303 291 }
304 292
305 293 /* assemble the result */
306 294 x->l.msw |= msw | (exp << 16);
307 295 x->l.frac2 = frac2;
308 296 x->l.frac3 = frac3;
309 297 x->l.frac4 = frac4;
310 298 }
311 299
312 -
313 300 /*
314 -* Compute the square root of x and place the TP result in s.
315 -*/
301 + * Compute the square root of x and place the TP result in s.
302 + */
316 303 static void
317 304 __q_tp_sqrt(const double *x, double *s)
318 305 {
319 - double c, rr, r[3], tt[3], t[5];
306 + double c, rr, r[3], tt[3], t[5];
320 307
321 308 /* approximate the divisor for the Newton iteration */
322 309 c = sqrt((x[0] + x[1]) + x[2]);
323 310 rr = half / c;
324 311
325 312 /* compute the first five "digits" of the square root */
326 313 t[0] = (c + two30) - two30;
327 314 tt[0] = t[0] + t[0];
328 315 r[0] = ((x[0] - t[0] * t[0]) + x[1]) + x[2];
329 316
330 317 t[1] = (rr * (r[0] + x[3]) + two6) - two6;
331 318 tt[1] = t[1] + t[1];
332 319 r[0] -= tt[0] * t[1];
333 320 r[1] = x[3] - t[1] * t[1];
334 321 c = (r[1] + twom18) - twom18;
335 322 r[0] += c;
336 323 r[1] = (r[1] - c) + x[4];
337 324
338 325 t[2] = (rr * (r[0] + r[1]) + twom18) - twom18;
339 326 tt[2] = t[2] + t[2];
340 327 r[0] -= tt[0] * t[2];
341 328 r[1] -= tt[1] * t[2];
342 329 c = (r[1] + twom42) - twom42;
343 330 r[0] += c;
344 331 r[1] = (r[1] - c) - t[2] * t[2];
345 332
|
↓ open down ↓ |
16 lines elided |
↑ open up ↑ |
346 333 t[3] = (rr * (r[0] + r[1]) + twom42) - twom42;
347 334 r[0] = ((r[0] - tt[0] * t[3]) + r[1]) - tt[1] * t[3];
348 335 r[1] = -tt[2] * t[3];
349 336 c = (r[1] + twom90) - twom90;
350 337 r[0] += c;
351 338 r[1] = (r[1] - c) - t[3] * t[3];
352 339
353 340 t[4] = (rr * (r[0] + r[1]) + twom66) - twom66;
354 341
355 342 /* here we just need to get the sign of the remainder */
356 - c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1])
357 - - tt[2] * t[4]) - (t[3] + t[3]) * t[4]) - t[4] * t[4];
343 + c = (((((r[0] - tt[0] * t[4]) - tt[1] * t[4]) + r[1]) - tt[2] * t[4]) -
344 + (t[3] + t[3]) * t[4]) - t[4] * t[4];
358 345
359 346 /* reduce to three doubles */
360 347 t[0] += t[1];
361 348 t[1] = t[2] + t[3];
362 349 t[2] = t[4];
363 350
364 351 /* if the third term might lie on a rounding boundary, perturb it */
365 - if (c != zero && t[2] == (twom62 + t[2]) - twom62)
366 - {
352 + if (c != zero && t[2] == (twom62 + t[2]) - twom62) {
367 353 if (c < zero)
368 354 t[2] -= twom124;
369 355 else
370 356 t[2] += twom124;
371 357 }
372 358
373 359 /* condense the square root */
374 360 c = t[1] + t[2];
375 361 t[2] += (t[1] - c);
376 362 t[1] = c;
377 363 c = t[0] + t[1];
378 364 s[1] = t[1] + (t[0] - c);
379 365 s[0] = c;
380 - if (s[1] == zero)
381 - {
366 +
367 + if (s[1] == zero) {
382 368 c = s[0] + t[2];
383 369 s[1] = t[2] + (s[0] - c);
384 370 s[0] = c;
385 371 s[2] = zero;
386 - }
387 - else
388 - {
372 + } else {
389 373 c = s[1] + t[2];
390 374 s[2] = t[2] + (s[1] - c);
391 375 s[1] = c;
392 376 }
393 377 }
394 378
395 -
396 379 long double
397 380 sqrtl(long double ldx)
398 381 {
399 - union longdouble x;
400 - volatile double t;
401 - double xx[5], zz[3];
402 - enum fp_direction_type rm;
403 - int ex, inexact, exc, traps;
382 + union longdouble x;
383 + volatile double t;
384 + double xx[5], zz[3];
385 + enum fp_direction_type rm;
386 + int ex, inexact, exc, traps;
404 387
405 388 /* clear cexc */
406 389 t = zero;
407 390 t -= zero;
408 391
409 392 /* check for zero operand */
410 393 x.d = ldx;
394 +
411 395 if (!((x.l.msw & 0x7fffffff) | x.l.frac2 | x.l.frac3 | x.l.frac4))
412 - return ldx;
396 + return (ldx);
413 397
414 398 /* handle nan and inf cases */
415 - if ((x.l.msw & 0x7fffffff) >= 0x7fff0000)
416 - {
417 - if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4)
418 - {
419 - if (!(x.l.msw & 0x8000))
420 - {
399 + if ((x.l.msw & 0x7fffffff) >= 0x7fff0000) {
400 + if ((x.l.msw & 0xffff) | x.l.frac2 | x.l.frac3 | x.l.frac4) {
401 + if (!(x.l.msw & 0x8000)) {
421 402 /* snan, signal invalid */
422 403 t += snan.d;
423 404 }
405 +
424 406 x.l.msw |= 0x8000;
425 - return x.d;
407 + return (x.d);
426 408 }
427 - if (x.l.msw & 0x80000000)
428 - {
409 +
410 + if (x.l.msw & 0x80000000) {
429 411 /* sqrt(-inf), signal invalid */
430 412 t = -one;
431 413 t = sqrt(t);
432 - return qnan.d;
414 + return (qnan.d);
433 415 }
416 +
434 417 /* sqrt(inf), return inf */
435 - return x.d;
418 + return (x.d);
436 419 }
437 420
438 421 /* handle negative numbers */
439 - if (x.l.msw & 0x80000000)
440 - {
422 + if (x.l.msw & 0x80000000) {
441 423 t = -one;
442 424 t = sqrt(t);
443 - return qnan.d;
425 + return (qnan.d);
444 426 }
445 427
446 428 /* now x is finite, positive */
447 429
448 430 traps = __swapTE(0);
449 431 exc = __swapEX(0);
450 432 rm = __swapRD(fp_nearest);
451 433
452 434 ex = __q_unpack(&x, xx);
453 - if (ex & 1)
454 - {
435 +
436 + if (ex & 1) {
455 437 /* make exponent even */
456 438 xx[0] += xx[0];
457 439 xx[1] += xx[1];
458 440 xx[2] += xx[2];
459 441 xx[3] += xx[3];
460 442 xx[4] += xx[4];
461 443 ex--;
462 444 }
445 +
463 446 __q_tp_sqrt(xx, zz);
464 447
465 448 /* put everything together */
466 449 x.l.msw = 0;
467 450 inexact = 0;
468 451 __q_pack(zz, ex >> 1, rm, &x, &inexact);
469 452
470 453 (void) __swapRD(rm);
471 454 (void) __swapEX(exc);
472 455 (void) __swapTE(traps);
473 - if (inexact)
474 - {
456 +
457 + if (inexact) {
475 458 t = huge;
476 459 t += tiny;
477 460 }
478 - return x.d;
461 +
462 + return (x.d);
479 463 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX