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