1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 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 #if defined(ELFOBJ) 31 #pragma weak fmal = __fmal 32 #endif 33 34 #include "libm.h" 35 #include "fma.h" 36 #include "fenv_inlines.h" 37 38 #if defined(__sparc) 39 40 static const union { 41 unsigned i[2]; 42 double d; 43 } C[] = { 44 { 0x3fe00000u, 0 }, 45 { 0x40000000u, 0 }, 46 { 0x3ef00000u, 0 }, 47 { 0x3e700000u, 0 }, 48 { 0x41300000u, 0 }, 49 { 0x3e300000u, 0 }, 50 { 0x3b300000u, 0 }, 51 { 0x38300000u, 0 }, 52 { 0x42300000u, 0 }, 53 { 0x3df00000u, 0 }, 54 { 0x7fe00000u, 0 }, 55 { 0x00100000u, 0 }, 56 { 0x00100001u, 0 }, 57 { 0, 0 }, 58 { 0x7ff00000u, 0 }, 59 { 0x7ff00001u, 0 } 60 }; 61 62 #define half C[0].d 63 #define two C[1].d 64 #define twom16 C[2].d 65 #define twom24 C[3].d 66 #define two20 C[4].d 67 #define twom28 C[5].d 68 #define twom76 C[6].d 69 #define twom124 C[7].d 70 #define two36 C[8].d 71 #define twom32 C[9].d 72 #define huge C[10].d 73 #define tiny C[11].d 74 #define tiny2 C[12].d 75 #define zero C[13].d 76 #define inf C[14].d 77 #define snan C[15].d 78 79 static const unsigned int fsr_rm = 0xc0000000u; 80 81 /* 82 * fmal for SPARC: 128-bit quad precision, big-endian 83 */ 84 long double 85 __fmal(long double x, long double y, long double z) { 86 union { 87 unsigned int i[4]; 88 long double q; 89 } xx, yy, zz; 90 union { 91 unsigned int i[2]; 92 double d; 93 } u; 94 double dx[5], dy[5], dxy[9], c, s; 95 unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7; 96 unsigned int z0, z1, z2, z3, z4, z5, z6, z7; 97 unsigned int rm, sticky; 98 unsigned int fsr; 99 int hx, hy, hz, ex, ey, ez, exy, sxy, sz, e, ibit; 100 int cx, cy, cz; 101 volatile double dummy; 102 103 /* extract the high order words of the arguments */ 104 xx.q = x; 105 yy.q = y; 106 zz.q = z; 107 hx = xx.i[0] & ~0x80000000; 108 hy = yy.i[0] & ~0x80000000; 109 hz = zz.i[0] & ~0x80000000; 110 111 /* 112 * distinguish zero, finite nonzero, infinite, and quiet nan 113 * arguments; raise invalid and return for signaling nans 114 */ 115 if (hx >= 0x7fff0000) { 116 if ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3]) { 117 if (!(hx & 0x8000)) { 118 /* signaling nan, raise invalid */ 119 dummy = snan; 120 dummy += snan; 121 xx.i[0] |= 0x8000; 122 return (xx.q); 123 } 124 cx = 3; /* quiet nan */ 125 } else 126 cx = 2; /* inf */ 127 } else if (hx == 0) { 128 cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0; 129 /* subnormal or zero */ 130 } else 131 cx = 1; /* finite nonzero */ 132 133 if (hy >= 0x7fff0000) { 134 if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) { 135 if (!(hy & 0x8000)) { 136 dummy = snan; 137 dummy += snan; 138 yy.i[0] |= 0x8000; 139 return (yy.q); 140 } 141 cy = 3; 142 } else 143 cy = 2; 144 } else if (hy == 0) { 145 cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0; 146 } else 147 cy = 1; 148 149 if (hz >= 0x7fff0000) { 150 if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) { 151 if (!(hz & 0x8000)) { 152 dummy = snan; 153 dummy += snan; 154 zz.i[0] |= 0x8000; 155 return (zz.q); 156 } 157 cz = 3; 158 } else 159 cz = 2; 160 } else if (hz == 0) { 161 cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0; 162 } else 163 cz = 1; 164 165 /* get the fsr and clear current exceptions */ 166 __fenv_getfsr32(&fsr); 167 fsr &= ~FSR_CEXC; 168 169 /* handle all other zero, inf, and nan cases */ 170 if (cx != 1 || cy != 1 || cz != 1) { 171 /* if x or y is a quiet nan, return it */ 172 if (cx == 3) { 173 __fenv_setfsr32(&fsr); 174 return (x); 175 } 176 if (cy == 3) { 177 __fenv_setfsr32(&fsr); 178 return (y); 179 } 180 181 /* if x*y is 0*inf, raise invalid and return the default nan */ 182 if ((cx == 0 && cy == 2) || (cx == 2 && cy == 0)) { 183 dummy = zero; 184 dummy *= inf; 185 zz.i[0] = 0x7fffffff; 186 zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff; 187 return (zz.q); 188 } 189 190 /* if z is a quiet nan, return it */ 191 if (cz == 3) { 192 __fenv_setfsr32(&fsr); 193 return (z); 194 } 195 196 /* 197 * now none of x, y, or z is nan; handle cases where x or y 198 * is inf 199 */ 200 if (cx == 2 || cy == 2) { 201 /* 202 * if z is also inf, either we have inf-inf or 203 * the result is the same as z depending on signs 204 */ 205 if (cz == 2) { 206 if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) { 207 dummy = inf; 208 dummy -= inf; 209 zz.i[0] = 0x7fffffff; 210 zz.i[1] = zz.i[2] = zz.i[3] = 211 0xffffffff; 212 return (zz.q); 213 } 214 __fenv_setfsr32(&fsr); 215 return (z); 216 } 217 218 /* otherwise the result is inf with appropriate sign */ 219 zz.i[0] = ((xx.i[0] ^ yy.i[0]) & 0x80000000) | 220 0x7fff0000; 221 zz.i[1] = zz.i[2] = zz.i[3] = 0; 222 __fenv_setfsr32(&fsr); 223 return (zz.q); 224 } 225 226 /* if z is inf, return it */ 227 if (cz == 2) { 228 __fenv_setfsr32(&fsr); 229 return (z); 230 } 231 232 /* 233 * now x, y, and z are all finite; handle cases where x or y 234 * is zero 235 */ 236 if (cx == 0 || cy == 0) { 237 /* either we have 0-0 or the result is the same as z */ 238 if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 239 0) { 240 zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 : 241 0; 242 __fenv_setfsr32(&fsr); 243 return (zz.q); 244 } 245 __fenv_setfsr32(&fsr); 246 return (z); 247 } 248 249 /* if we get here, x and y are nonzero finite, z must be zero */ 250 return (x * y); 251 } 252 253 /* 254 * now x, y, and z are all finite and nonzero; set round-to- 255 * negative-infinity mode 256 */ 257 __fenv_setfsr32(&fsr_rm); 258 259 /* 260 * get the signs and exponents and normalize the significands 261 * of x and y 262 */ 263 sxy = (xx.i[0] ^ yy.i[0]) & 0x80000000; 264 ex = hx >> 16; 265 hx &= 0xffff; 266 if (!ex) { 267 if (hx | (xx.i[1] & 0xfffe0000)) { 268 ex = 1; 269 } else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) { 270 hx = xx.i[1]; 271 xx.i[1] = xx.i[2]; 272 xx.i[2] = xx.i[3]; 273 xx.i[3] = 0; 274 ex = -31; 275 } else if (xx.i[2] | (xx.i[3] & 0xfffe0000)) { 276 hx = xx.i[2]; 277 xx.i[1] = xx.i[3]; 278 xx.i[2] = xx.i[3] = 0; 279 ex = -63; 280 } else { 281 hx = xx.i[3]; 282 xx.i[1] = xx.i[2] = xx.i[3] = 0; 283 ex = -95; 284 } 285 while ((hx & 0x10000) == 0) { 286 hx = (hx << 1) | (xx.i[1] >> 31); 287 xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31); 288 xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31); 289 xx.i[3] <<= 1; 290 ex--; 291 } 292 } else 293 hx |= 0x10000; 294 ey = hy >> 16; 295 hy &= 0xffff; 296 if (!ey) { 297 if (hy | (yy.i[1] & 0xfffe0000)) { 298 ey = 1; 299 } else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) { 300 hy = yy.i[1]; 301 yy.i[1] = yy.i[2]; 302 yy.i[2] = yy.i[3]; 303 yy.i[3] = 0; 304 ey = -31; 305 } else if (yy.i[2] | (yy.i[3] & 0xfffe0000)) { 306 hy = yy.i[2]; 307 yy.i[1] = yy.i[3]; 308 yy.i[2] = yy.i[3] = 0; 309 ey = -63; 310 } else { 311 hy = yy.i[3]; 312 yy.i[1] = yy.i[2] = yy.i[3] = 0; 313 ey = -95; 314 } 315 while ((hy & 0x10000) == 0) { 316 hy = (hy << 1) | (yy.i[1] >> 31); 317 yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31); 318 yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31); 319 yy.i[3] <<= 1; 320 ey--; 321 } 322 } else 323 hy |= 0x10000; 324 exy = ex + ey - 0x3fff; 325 326 /* convert the significands of x and y to doubles */ 327 c = twom16; 328 dx[0] = (double) ((int) hx) * c; 329 dy[0] = (double) ((int) hy) * c; 330 331 c *= twom24; 332 dx[1] = (double) ((int) (xx.i[1] >> 8)) * c; 333 dy[1] = (double) ((int) (yy.i[1] >> 8)) * c; 334 335 c *= twom24; 336 dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) & 337 0xffffff)) * c; 338 dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) & 339 0xffffff)) * c; 340 341 c *= twom24; 342 dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) & 343 0xffffff)) * c; 344 dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) & 345 0xffffff)) * c; 346 347 c *= twom24; 348 dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c; 349 dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c; 350 351 /* form the "digits" of the product */ 352 dxy[0] = dx[0] * dy[0]; 353 dxy[1] = dx[0] * dy[1] + dx[1] * dy[0]; 354 dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0]; 355 dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] + 356 dx[3] * dy[0]; 357 dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] + 358 dx[3] * dy[1] + dx[4] * dy[0]; 359 dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] + 360 dx[4] * dy[1]; 361 dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2]; 362 dxy[7] = dx[3] * dy[4] + dx[4] * dy[3]; 363 dxy[8] = dx[4] * dy[4]; 364 365 /* split odd-numbered terms and combine into even-numbered terms */ 366 c = (dxy[1] + two20) - two20; 367 dxy[0] += c; 368 dxy[1] -= c; 369 c = (dxy[3] + twom28) - twom28; 370 dxy[2] += c + dxy[1]; 371 dxy[3] -= c; 372 c = (dxy[5] + twom76) - twom76; 373 dxy[4] += c + dxy[3]; 374 dxy[5] -= c; 375 c = (dxy[7] + twom124) - twom124; 376 dxy[6] += c + dxy[5]; 377 dxy[8] += (dxy[7] - c); 378 379 /* propagate carries, adjusting the exponent if need be */ 380 dxy[7] = dxy[6] + dxy[8]; 381 dxy[5] = dxy[4] + dxy[7]; 382 dxy[3] = dxy[2] + dxy[5]; 383 dxy[1] = dxy[0] + dxy[3]; 384 if (dxy[1] >= two) { 385 dxy[0] *= half; 386 dxy[1] *= half; 387 dxy[2] *= half; 388 dxy[3] *= half; 389 dxy[4] *= half; 390 dxy[5] *= half; 391 dxy[6] *= half; 392 dxy[7] *= half; 393 dxy[8] *= half; 394 exy++; 395 } 396 397 /* extract the significand of x*y */ 398 s = two36; 399 u.d = c = dxy[1] + s; 400 xy0 = u.i[1]; 401 c -= s; 402 dxy[1] -= c; 403 dxy[0] -= c; 404 405 s *= twom32; 406 u.d = c = dxy[1] + s; 407 xy1 = u.i[1]; 408 c -= s; 409 dxy[2] += (dxy[0] - c); 410 dxy[3] = dxy[2] + dxy[5]; 411 412 s *= twom32; 413 u.d = c = dxy[3] + s; 414 xy2 = u.i[1]; 415 c -= s; 416 dxy[4] += (dxy[2] - c); 417 dxy[5] = dxy[4] + dxy[7]; 418 419 s *= twom32; 420 u.d = c = dxy[5] + s; 421 xy3 = u.i[1]; 422 c -= s; 423 dxy[4] -= c; 424 dxy[5] = dxy[4] + dxy[7]; 425 426 s *= twom32; 427 u.d = c = dxy[5] + s; 428 xy4 = u.i[1]; 429 c -= s; 430 dxy[6] += (dxy[4] - c); 431 dxy[7] = dxy[6] + dxy[8]; 432 433 s *= twom32; 434 u.d = c = dxy[7] + s; 435 xy5 = u.i[1]; 436 c -= s; 437 dxy[8] += (dxy[6] - c); 438 439 s *= twom32; 440 u.d = c = dxy[8] + s; 441 xy6 = u.i[1]; 442 c -= s; 443 dxy[8] -= c; 444 445 s *= twom32; 446 u.d = c = dxy[8] + s; 447 xy7 = u.i[1]; 448 449 /* extract the sign, exponent, and significand of z */ 450 sz = zz.i[0] & 0x80000000; 451 ez = hz >> 16; 452 z0 = hz & 0xffff; 453 if (!ez) { 454 if (z0 | (zz.i[1] & 0xfffe0000)) { 455 z1 = zz.i[1]; 456 z2 = zz.i[2]; 457 z3 = zz.i[3]; 458 ez = 1; 459 } else if (zz.i[1] | (zz.i[2] & 0xfffe0000)) { 460 z0 = zz.i[1]; 461 z1 = zz.i[2]; 462 z2 = zz.i[3]; 463 z3 = 0; 464 ez = -31; 465 } else if (zz.i[2] | (zz.i[3] & 0xfffe0000)) { 466 z0 = zz.i[2]; 467 z1 = zz.i[3]; 468 z2 = z3 = 0; 469 ez = -63; 470 } else { 471 z0 = zz.i[3]; 472 z1 = z2 = z3 = 0; 473 ez = -95; 474 } 475 while ((z0 & 0x10000) == 0) { 476 z0 = (z0 << 1) | (z1 >> 31); 477 z1 = (z1 << 1) | (z2 >> 31); 478 z2 = (z2 << 1) | (z3 >> 31); 479 z3 <<= 1; 480 ez--; 481 } 482 } else { 483 z0 |= 0x10000; 484 z1 = zz.i[1]; 485 z2 = zz.i[2]; 486 z3 = zz.i[3]; 487 } 488 z4 = z5 = z6 = z7 = 0; 489 490 /* 491 * now x*y is represented by sxy, exy, and xy[0-7], and z is 492 * represented likewise; swap if need be so |xy| <= |z| 493 */ 494 if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 || 495 (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 || 496 (xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) { 497 e = sxy; sxy = sz; sz = e; 498 e = exy; exy = ez; ez = e; 499 e = xy0; xy0 = z0; z0 = e; 500 e = xy1; xy1 = z1; z1 = e; 501 e = xy2; xy2 = z2; z2 = e; 502 e = xy3; xy3 = z3; z3 = e; 503 z4 = xy4; xy4 = 0; 504 z5 = xy5; xy5 = 0; 505 z6 = xy6; xy6 = 0; 506 z7 = xy7; xy7 = 0; 507 } 508 509 /* shift the significand of xy keeping a sticky bit */ 510 e = ez - exy; 511 if (e > 236) { 512 xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0; 513 xy7 = 1; 514 } else if (e >= 224) { 515 sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 | 516 ((xy0 << 1) << (255 - e)); 517 xy7 = xy0 >> (e - 224); 518 if (sticky) 519 xy7 |= 1; 520 xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0; 521 } else if (e >= 192) { 522 sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | 523 ((xy1 << 1) << (223 - e)); 524 xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e)); 525 if (sticky) 526 xy7 |= 1; 527 xy6 = xy0 >> (e - 192); 528 xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0; 529 } else if (e >= 160) { 530 sticky = xy7 | xy6 | xy5 | xy4 | xy3 | 531 ((xy2 << 1) << (191 - e)); 532 xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e)); 533 if (sticky) 534 xy7 |= 1; 535 xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e)); 536 xy5 = xy0 >> (e - 160); 537 xy0 = xy1 = xy2 = xy3 = xy4 = 0; 538 } else if (e >= 128) { 539 sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e)); 540 xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e)); 541 if (sticky) 542 xy7 |= 1; 543 xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e)); 544 xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e)); 545 xy4 = xy0 >> (e - 128); 546 xy0 = xy1 = xy2 = xy3 = 0; 547 } else if (e >= 96) { 548 sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e)); 549 xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e)); 550 if (sticky) 551 xy7 |= 1; 552 xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e)); 553 xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e)); 554 xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e)); 555 xy3 = xy0 >> (e - 96); 556 xy0 = xy1 = xy2 = 0; 557 } else if (e >= 64) { 558 sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e)); 559 xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e)); 560 if (sticky) 561 xy7 |= 1; 562 xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e)); 563 xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e)); 564 xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e)); 565 xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e)); 566 xy2 = xy0 >> (e - 64); 567 xy0 = xy1 = 0; 568 } else if (e >= 32) { 569 sticky = xy7 | ((xy6 << 1) << (63 - e)); 570 xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e)); 571 if (sticky) 572 xy7 |= 1; 573 xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e)); 574 xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e)); 575 xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e)); 576 xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e)); 577 xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e)); 578 xy1 = xy0 >> (e - 32); 579 xy0 = 0; 580 } else if (e) { 581 sticky = (xy7 << 1) << (31 - e); 582 xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e)); 583 if (sticky) 584 xy7 |= 1; 585 xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e)); 586 xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e)); 587 xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e)); 588 xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e)); 589 xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e)); 590 xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e)); 591 xy0 >>= e; 592 } 593 594 /* if this is a magnitude subtract, negate the significand of xy */ 595 if (sxy ^ sz) { 596 xy0 = ~xy0; 597 xy1 = ~xy1; 598 xy2 = ~xy2; 599 xy3 = ~xy3; 600 xy4 = ~xy4; 601 xy5 = ~xy5; 602 xy6 = ~xy6; 603 xy7 = -xy7; 604 if (xy7 == 0) 605 if (++xy6 == 0) 606 if (++xy5 == 0) 607 if (++xy4 == 0) 608 if (++xy3 == 0) 609 if (++xy2 == 0) 610 if (++xy1 == 0) 611 xy0++; 612 } 613 614 /* add, propagating carries */ 615 z7 += xy7; 616 e = (z7 < xy7); 617 z6 += xy6; 618 if (e) { 619 z6++; 620 e = (z6 <= xy6); 621 } else 622 e = (z6 < xy6); 623 z5 += xy5; 624 if (e) { 625 z5++; 626 e = (z5 <= xy5); 627 } else 628 e = (z5 < xy5); 629 z4 += xy4; 630 if (e) { 631 z4++; 632 e = (z4 <= xy4); 633 } else 634 e = (z4 < xy4); 635 z3 += xy3; 636 if (e) { 637 z3++; 638 e = (z3 <= xy3); 639 } else 640 e = (z3 < xy3); 641 z2 += xy2; 642 if (e) { 643 z2++; 644 e = (z2 <= xy2); 645 } else 646 e = (z2 < xy2); 647 z1 += xy1; 648 if (e) { 649 z1++; 650 e = (z1 <= xy1); 651 } else 652 e = (z1 < xy1); 653 z0 += xy0; 654 if (e) 655 z0++; 656 657 /* postnormalize and collect rounding information into z4 */ 658 if (ez < 1) { 659 /* result is tiny; shift right until exponent is within range */ 660 e = 1 - ez; 661 if (e > 116) { 662 z4 = 1; /* result can't be exactly zero */ 663 z0 = z1 = z2 = z3 = 0; 664 } else if (e >= 96) { 665 sticky = z7 | z6 | z5 | z4 | z3 | z2 | 666 ((z1 << 1) << (127 - e)); 667 z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e)); 668 if (sticky) 669 z4 |= 1; 670 z3 = z0 >> (e - 96); 671 z0 = z1 = z2 = 0; 672 } else if (e >= 64) { 673 sticky = z7 | z6 | z5 | z4 | z3 | 674 ((z2 << 1) << (95 - e)); 675 z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e)); 676 if (sticky) 677 z4 |= 1; 678 z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e)); 679 z2 = z0 >> (e - 64); 680 z0 = z1 = 0; 681 } else if (e >= 32) { 682 sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e)); 683 z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e)); 684 if (sticky) 685 z4 |= 1; 686 z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e)); 687 z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e)); 688 z1 = z0 >> (e - 32); 689 z0 = 0; 690 } else { 691 sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e); 692 z4 = (z4 >> e) | ((z3 << 1) << (31 - e)); 693 if (sticky) 694 z4 |= 1; 695 z3 = (z3 >> e) | ((z2 << 1) << (31 - e)); 696 z2 = (z2 >> e) | ((z1 << 1) << (31 - e)); 697 z1 = (z1 >> e) | ((z0 << 1) << (31 - e)); 698 z0 >>= e; 699 } 700 ez = 1; 701 } else if (z0 >= 0x20000) { 702 /* carry out; shift right by one */ 703 sticky = (z4 & 1) | z5 | z6 | z7; 704 z4 = (z4 >> 1) | (z3 << 31); 705 if (sticky) 706 z4 |= 1; 707 z3 = (z3 >> 1) | (z2 << 31); 708 z2 = (z2 >> 1) | (z1 << 31); 709 z1 = (z1 >> 1) | (z0 << 31); 710 z0 >>= 1; 711 ez++; 712 } else { 713 if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7) 714 != 0) { 715 /* 716 * borrow/cancellation; shift left as much as 717 * exponent allows 718 */ 719 while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) { 720 z0 = z1; 721 z1 = z2; 722 z2 = z3; 723 z3 = z4; 724 z4 = z5; 725 z5 = z6; 726 z6 = z7; 727 z7 = 0; 728 ez -= 32; 729 } 730 while (z0 < 0x10000 && ez > 1) { 731 z0 = (z0 << 1) | (z1 >> 31); 732 z1 = (z1 << 1) | (z2 >> 31); 733 z2 = (z2 << 1) | (z3 >> 31); 734 z3 = (z3 << 1) | (z4 >> 31); 735 z4 = (z4 << 1) | (z5 >> 31); 736 z5 = (z5 << 1) | (z6 >> 31); 737 z6 = (z6 << 1) | (z7 >> 31); 738 z7 <<= 1; 739 ez--; 740 } 741 } 742 if (z5 | z6 | z7) 743 z4 |= 1; 744 } 745 746 /* get the rounding mode */ 747 rm = fsr >> 30; 748 749 /* strip off the integer bit, if there is one */ 750 ibit = z0 & 0x10000; 751 if (ibit) 752 z0 -= 0x10000; 753 else { 754 ez = 0; 755 if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */ 756 zz.i[0] = rm == FSR_RM ? 0x80000000 : 0; 757 zz.i[1] = zz.i[2] = zz.i[3] = 0; 758 __fenv_setfsr32(&fsr); 759 return (zz.q); 760 } 761 } 762 763 /* 764 * flip the sense of directed roundings if the result is negative; 765 * the logic below applies to a positive result 766 */ 767 if (sz) 768 rm ^= rm >> 1; 769 770 /* round and raise exceptions */ 771 if (z4) { 772 fsr |= FSR_NXC; 773 774 /* decide whether to round the fraction up */ 775 if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u || 776 (z4 == 0x80000000u && (z3 & 1))))) { 777 /* round up and renormalize if necessary */ 778 if (++z3 == 0) 779 if (++z2 == 0) 780 if (++z1 == 0) 781 if (++z0 == 0x10000) { 782 z0 = 0; 783 ez++; 784 } 785 } 786 } 787 788 /* check for under/overflow */ 789 if (ez >= 0x7fff) { 790 if (rm == FSR_RN || rm == FSR_RP) { 791 zz.i[0] = sz | 0x7fff0000; 792 zz.i[1] = zz.i[2] = zz.i[3] = 0; 793 } else { 794 zz.i[0] = sz | 0x7ffeffff; 795 zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff; 796 } 797 fsr |= FSR_OFC | FSR_NXC; 798 } else { 799 zz.i[0] = sz | (ez << 16) | z0; 800 zz.i[1] = z1; 801 zz.i[2] = z2; 802 zz.i[3] = z3; 803 804 /* 805 * !ibit => exact result was tiny before rounding, 806 * z4 nonzero => result delivered is inexact 807 */ 808 if (!ibit) { 809 if (z4) 810 fsr |= FSR_UFC | FSR_NXC; 811 else if (fsr & FSR_UFM) 812 fsr |= FSR_UFC; 813 } 814 } 815 816 /* restore the fsr and emulate exceptions as needed */ 817 if ((fsr & FSR_CEXC) & (fsr >> 23)) { 818 __fenv_setfsr32(&fsr); 819 if (fsr & FSR_OFC) { 820 dummy = huge; 821 dummy *= huge; 822 } else if (fsr & FSR_UFC) { 823 dummy = tiny; 824 if (fsr & FSR_NXC) 825 dummy *= tiny; 826 else 827 dummy -= tiny2; 828 } else { 829 dummy = huge; 830 dummy += tiny; 831 } 832 } else { 833 fsr |= (fsr & 0x1f) << 5; 834 __fenv_setfsr32(&fsr); 835 } 836 return (zz.q); 837 } 838 839 #elif defined(__x86) 840 841 static const union { 842 unsigned i[2]; 843 double d; 844 } C[] = { 845 { 0, 0x3fe00000u }, 846 { 0, 0x40000000u }, 847 { 0, 0x3df00000u }, 848 { 0, 0x3bf00000u }, 849 { 0, 0x41f00000u }, 850 { 0, 0x43e00000u }, 851 { 0, 0x7fe00000u }, 852 { 0, 0x00100000u }, 853 { 0, 0x00100001u } 854 }; 855 856 #define half C[0].d 857 #define two C[1].d 858 #define twom32 C[2].d 859 #define twom64 C[3].d 860 #define two32 C[4].d 861 #define two63 C[5].d 862 #define huge C[6].d 863 #define tiny C[7].d 864 #define tiny2 C[8].d 865 866 #if defined(__amd64) 867 #define NI 4 868 #else 869 #define NI 3 870 #endif 871 872 /* 873 * fmal for x86: 80-bit extended double precision, little-endian 874 */ 875 long double 876 __fmal(long double x, long double y, long double z) { 877 union { 878 unsigned i[NI]; 879 long double e; 880 } xx, yy, zz; 881 long double xhi, yhi, xlo, ylo, t; 882 unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4; 883 unsigned oldcwsw, cwsw, rm, sticky, carry; 884 int ex, ey, ez, exy, sxy, sz, e, tinyafter; 885 volatile double dummy; 886 887 /* extract the exponents of the arguments */ 888 xx.e = x; 889 yy.e = y; 890 zz.e = z; 891 ex = xx.i[2] & 0x7fff; 892 ey = yy.i[2] & 0x7fff; 893 ez = zz.i[2] & 0x7fff; 894 895 /* dispense with inf, nan, and zero cases */ 896 if (ex == 0x7fff || ey == 0x7fff || (ex | xx.i[1] | xx.i[0]) == 0 || 897 (ey | yy.i[1] | yy.i[0]) == 0) /* x or y is inf, nan, or 0 */ 898 return (x * y + z); 899 900 if (ez == 0x7fff) /* z is inf or nan */ 901 return (x + z); /* avoid spurious under/overflow in x * y */ 902 903 if ((ez | zz.i[1] | zz.i[0]) == 0) /* z is zero */ 904 /* 905 * x * y isn't zero but could underflow to zero, 906 * so don't add z, lest we perturb the sign 907 */ 908 return (x * y); 909 910 /* 911 * now x, y, and z are all finite and nonzero; extract signs and 912 * normalize the significands (this will raise the denormal operand 913 * exception if need be) 914 */ 915 sxy = (xx.i[2] ^ yy.i[2]) & 0x8000; 916 sz = zz.i[2] & 0x8000; 917 if (!ex) { 918 xx.e = x * two63; 919 ex = (xx.i[2] & 0x7fff) - 63; 920 } 921 if (!ey) { 922 yy.e = y * two63; 923 ey = (yy.i[2] & 0x7fff) - 63; 924 } 925 if (!ez) { 926 zz.e = z * two63; 927 ez = (zz.i[2] & 0x7fff) - 63; 928 } 929 930 /* 931 * save the control and status words, mask all exceptions, and 932 * set rounding to 64-bit precision and toward-zero 933 */ 934 __fenv_getcwsw(&oldcwsw); 935 cwsw = (oldcwsw & 0xf0c0ffff) | 0x0f3f0000; 936 __fenv_setcwsw(&cwsw); 937 938 /* multiply x*y to 128 bits */ 939 exy = ex + ey - 0x3fff; 940 xx.i[2] = 0x3fff; 941 yy.i[2] = 0x3fff; 942 x = xx.e; 943 y = yy.e; 944 xhi = ((x + twom32) + two32) - two32; 945 yhi = ((y + twom32) + two32) - two32; 946 xlo = x - xhi; 947 ylo = y - yhi; 948 x *= y; 949 y = ((xhi * yhi - x) + xhi * ylo + xlo * yhi) + xlo * ylo; 950 if (x >= two) { 951 x *= half; 952 y *= half; 953 exy++; 954 } 955 956 /* extract the significands */ 957 xx.e = x; 958 xy0 = xx.i[1]; 959 xy1 = xx.i[0]; 960 yy.e = t = y + twom32; 961 xy2 = yy.i[0]; 962 yy.e = (y - (t - twom32)) + twom64; 963 xy3 = yy.i[0]; 964 xy4 = 0; 965 z0 = zz.i[1]; 966 z1 = zz.i[0]; 967 z2 = z3 = z4 = 0; 968 969 /* 970 * now x*y is represented by sxy, exy, and xy[0-4], and z is 971 * represented likewise; swap if need be so |xy| <= |z| 972 */ 973 if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && 974 (xy1 > z1 || (xy1 == z1 && (xy2 | xy3) != 0)))))) { 975 e = sxy; sxy = sz; sz = e; 976 e = exy; exy = ez; ez = e; 977 e = xy0; xy0 = z0; z0 = e; 978 e = xy1; xy1 = z1; z1 = e; 979 z2 = xy2; xy2 = 0; 980 z3 = xy3; xy3 = 0; 981 } 982 983 /* shift the significand of xy keeping a sticky bit */ 984 e = ez - exy; 985 if (e > 130) { 986 xy0 = xy1 = xy2 = xy3 = 0; 987 xy4 = 1; 988 } else if (e >= 128) { 989 sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e)); 990 xy4 = xy0 >> (e - 128); 991 if (sticky) 992 xy4 |= 1; 993 xy0 = xy1 = xy2 = xy3 = 0; 994 } else if (e >= 96) { 995 sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e)); 996 xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e)); 997 if (sticky) 998 xy4 |= 1; 999 xy3 = xy0 >> (e - 96); 1000 xy0 = xy1 = xy2 = 0; 1001 } else if (e >= 64) { 1002 sticky = xy3 | ((xy2 << 1) << (95 - e)); 1003 xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e)); 1004 if (sticky) 1005 xy4 |= 1; 1006 xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e)); 1007 xy2 = xy0 >> (e - 64); 1008 xy0 = xy1 = 0; 1009 } else if (e >= 32) { 1010 sticky = (xy3 << 1) << (63 - e); 1011 xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e)); 1012 if (sticky) 1013 xy4 |= 1; 1014 xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e)); 1015 xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e)); 1016 xy1 = xy0 >> (e - 32); 1017 xy0 = 0; 1018 } else if (e) { 1019 xy4 = (xy3 << 1) << (31 - e); 1020 xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e)); 1021 xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e)); 1022 xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e)); 1023 xy0 >>= e; 1024 } 1025 1026 /* if this is a magnitude subtract, negate the significand of xy */ 1027 if (sxy ^ sz) { 1028 xy0 = ~xy0; 1029 xy1 = ~xy1; 1030 xy2 = ~xy2; 1031 xy3 = ~xy3; 1032 xy4 = -xy4; 1033 if (xy4 == 0) 1034 if (++xy3 == 0) 1035 if (++xy2 == 0) 1036 if (++xy1 == 0) 1037 xy0++; 1038 } 1039 1040 /* add, propagating carries */ 1041 z4 += xy4; 1042 carry = (z4 < xy4); 1043 z3 += xy3; 1044 if (carry) { 1045 z3++; 1046 carry = (z3 <= xy3); 1047 } else 1048 carry = (z3 < xy3); 1049 z2 += xy2; 1050 if (carry) { 1051 z2++; 1052 carry = (z2 <= xy2); 1053 } else 1054 carry = (z2 < xy2); 1055 z1 += xy1; 1056 if (carry) { 1057 z1++; 1058 carry = (z1 <= xy1); 1059 } else 1060 carry = (z1 < xy1); 1061 z0 += xy0; 1062 if (carry) { 1063 z0++; 1064 carry = (z0 <= xy0); 1065 } else 1066 carry = (z0 < xy0); 1067 1068 /* for a magnitude subtract, ignore the last carry out */ 1069 if (sxy ^ sz) 1070 carry = 0; 1071 1072 /* postnormalize and collect rounding information into z2 */ 1073 if (ez < 1) { 1074 /* result is tiny; shift right until exponent is within range */ 1075 e = 1 - ez; 1076 if (e > 67) { 1077 z2 = 1; /* result can't be exactly zero */ 1078 z0 = z1 = 0; 1079 } else if (e >= 64) { 1080 sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e)); 1081 z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e)); 1082 if (sticky) 1083 z2 |= 1; 1084 z1 = carry >> (e - 64); 1085 z0 = 0; 1086 } else if (e >= 32) { 1087 sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e)); 1088 z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e)); 1089 if (sticky) 1090 z2 |= 1; 1091 z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e)); 1092 z0 = carry >> (e - 32); 1093 } else { 1094 sticky = z4 | z3 | (z2 << 1) << (31 - e); 1095 z2 = (z2 >> e) | ((z1 << 1) << (31 - e)); 1096 if (sticky) 1097 z2 |= 1; 1098 z1 = (z1 >> e) | ((z0 << 1) << (31 - e)); 1099 z0 = (z0 >> e) | ((carry << 1) << (31 - e)); 1100 } 1101 ez = 1; 1102 } else if (carry) { 1103 /* carry out; shift right by one */ 1104 sticky = (z2 & 1) | z3 | z4; 1105 z2 = (z2 >> 1) | (z1 << 31); 1106 if (sticky) 1107 z2 |= 1; 1108 z1 = (z1 >> 1) | (z0 << 31); 1109 z0 = (z0 >> 1) | 0x80000000; 1110 ez++; 1111 } else { 1112 if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) { 1113 /* 1114 * borrow/cancellation; shift left as much as 1115 * exponent allows 1116 */ 1117 while (!z0 && ez >= 33) { 1118 z0 = z1; 1119 z1 = z2; 1120 z2 = z3; 1121 z3 = z4; 1122 z4 = 0; 1123 ez -= 32; 1124 } 1125 while (z0 < 0x80000000u && ez > 1) { 1126 z0 = (z0 << 1) | (z1 >> 31); 1127 z1 = (z1 << 1) | (z2 >> 31); 1128 z2 = (z2 << 1) | (z3 >> 31); 1129 z3 = (z3 << 1) | (z4 >> 31); 1130 z4 <<= 1; 1131 ez--; 1132 } 1133 } 1134 if (z3 | z4) 1135 z2 |= 1; 1136 } 1137 1138 /* get the rounding mode */ 1139 rm = oldcwsw & 0x0c000000; 1140 1141 /* adjust exponent if result is subnormal */ 1142 tinyafter = 0; 1143 if (!(z0 & 0x80000000)) { 1144 ez = 0; 1145 tinyafter = 1; 1146 if (!(z0 | z1 | z2)) { /* exact zero */ 1147 zz.i[2] = rm == FCW_RM ? 0x8000 : 0; 1148 zz.i[1] = zz.i[0] = 0; 1149 __fenv_setcwsw(&oldcwsw); 1150 return (zz.e); 1151 } 1152 } 1153 1154 /* 1155 * flip the sense of directed roundings if the result is negative; 1156 * the logic below applies to a positive result 1157 */ 1158 if (sz && (rm == FCW_RM || rm == FCW_RP)) 1159 rm = (FCW_RM + FCW_RP) - rm; 1160 1161 /* round */ 1162 if (z2) { 1163 if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u || 1164 (z2 == 0x80000000u && (z1 & 1))))) { 1165 /* round up and renormalize if necessary */ 1166 if (++z1 == 0) { 1167 if (++z0 == 0) { 1168 z0 = 0x80000000; 1169 ez++; 1170 } else if (z0 == 0x80000000) { 1171 /* rounded up to smallest normal */ 1172 ez = 1; 1173 if ((rm == FCW_RP && z2 > 1174 0x80000000u) || (rm == FCW_RN && 1175 z2 >= 0xc0000000u)) 1176 /* 1177 * would have rounded up to 1178 * smallest normal even with 1179 * unbounded range 1180 */ 1181 tinyafter = 0; 1182 } 1183 } 1184 } 1185 } 1186 1187 /* restore the control and status words, check for over/underflow */ 1188 __fenv_setcwsw(&oldcwsw); 1189 if (ez >= 0x7fff) { 1190 if (rm == FCW_RN || rm == FCW_RP) { 1191 zz.i[2] = sz | 0x7fff; 1192 zz.i[1] = 0x80000000; 1193 zz.i[0] = 0; 1194 } else { 1195 zz.i[2] = sz | 0x7ffe; 1196 zz.i[1] = 0xffffffff; 1197 zz.i[0] = 0xffffffff; 1198 } 1199 dummy = huge; 1200 dummy *= huge; 1201 } else { 1202 zz.i[2] = sz | ez; 1203 zz.i[1] = z0; 1204 zz.i[0] = z1; 1205 1206 /* 1207 * tinyafter => result rounded w/ unbounded range would be tiny, 1208 * z2 nonzero => result delivered is inexact 1209 */ 1210 if (tinyafter) { 1211 dummy = tiny; 1212 if (z2) 1213 dummy *= tiny; 1214 else 1215 dummy -= tiny2; 1216 } else if (z2) { 1217 dummy = huge; 1218 dummy += tiny; 1219 } 1220 } 1221 1222 return (zz.e); 1223 } 1224 1225 #else 1226 #error Unknown architecture 1227 #endif 1228