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, Version 1.0 only 6 * (the "License"). You may not use this file except in compliance 7 * with the License. 8 * 9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10 * or http://www.opensolaris.org/os/licensing. 11 * See the License for the specific language governing permissions 12 * and limitations under the License. 13 * 14 * When distributing Covered Code, include this CDDL HEADER in each 15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16 * If applicable, add the following below this CDDL HEADER, with the 17 * fields enclosed by brackets "[]" replaced with your own identifying 18 * information: Portions Copyright [yyyy] [name of copyright owner] 19 * 20 * CDDL HEADER END 21 */ 22 /* 23 * Copyright 2003 Sun Microsystems, Inc. All rights reserved. 24 * Use is subject to license terms. 25 */ 26 27 #include "quad.h" 28 29 static const double C[] = { 30 0.0, 31 1.0, 32 68719476736.0, 33 402653184.0, 34 24.0, 35 16.0, 36 3.66210937500000000000e-04, 37 1.52587890625000000000e-05, 38 1.43051147460937500000e-06, 39 5.96046447753906250000e-08, 40 3.72529029846191406250e-09, 41 2.18278728425502777100e-11, 42 8.52651282912120223045e-14, 43 3.55271367880050092936e-15, 44 1.30104260698260532081e-18, 45 8.67361737988403547206e-19, 46 2.16840434497100886801e-19, 47 3.17637355220362627151e-22, 48 7.75481824268463445192e-26, 49 4.62223186652936604733e-33, 50 9.62964972193617926528e-35, 51 4.70197740328915003187e-38, 52 2.75506488473973634680e-40, 53 }; 54 55 #define zero C[0] 56 #define one C[1] 57 #define two36 C[2] 58 #define three2p27 C[3] 59 #define three2p3 C[4] 60 #define two4 C[5] 61 #define three2m13 C[6] 62 #define twom16 C[7] 63 #define three2m21 C[8] 64 #define twom24 C[9] 65 #define twom28 C[10] 66 #define three2m37 C[11] 67 #define three2m45 C[12] 68 #define twom48 C[13] 69 #define three2m61 C[14] 70 #define twom60 C[15] 71 #define twom62 C[16] 72 #define three2m73 C[17] 73 #define three2m85 C[18] 74 #define three2m109 C[19] 75 #define twom113 C[20] 76 #define twom124 C[21] 77 #define three2m133 C[22] 78 79 static const unsigned int 80 fsr_re = 0x00000000u, 81 fsr_rn = 0xc0000000u; 82 83 #ifdef __sparcv9 84 85 /* 86 * _Qp_div(pz, x, y) sets *pz = *x / *y. 87 */ 88 void 89 _Qp_div(union longdouble *pz, const union longdouble *x, 90 const union longdouble *y) 91 92 #else 93 94 /* 95 * _Q_div(x, y) returns *x / *y. 96 */ 97 union longdouble 98 _Q_div(const union longdouble *x, const union longdouble *y) 99 100 #endif /* __sparcv9 */ 101 102 { 103 union longdouble z; 104 union xdouble u; 105 double c, d, ry, xx[4], yy[5], zz[5]; 106 unsigned int xm, ym, fsr, lx, ly, wx[3], wy[3]; 107 unsigned int msw, frac2, frac3, frac4, rm; 108 int ibit, ex, ey, ez, sign; 109 110 xm = x->l.msw & 0x7fffffff; 111 ym = y->l.msw & 0x7fffffff; 112 sign = (x->l.msw ^ y->l.msw) & ~0x7fffffff; 113 114 __quad_getfsrp(&fsr); 115 116 /* handle nan and inf cases */ 117 if (xm >= 0x7fff0000 || ym >= 0x7fff0000) { 118 /* handle nan cases according to V9 app. B */ 119 if (QUAD_ISNAN(*y)) { 120 if (!(y->l.msw & 0x8000)) { 121 /* snan, signal invalid */ 122 if (fsr & FSR_NVM) { 123 __quad_fdivq(x, y, &Z); 124 } else { 125 Z = *y; 126 Z.l.msw |= 0x8000; 127 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | 128 FSR_NVC; 129 __quad_setfsrp(&fsr); 130 } 131 QUAD_RETURN(Z); 132 } else if (QUAD_ISNAN(*x) && !(x->l.msw & 0x8000)) { 133 /* snan, signal invalid */ 134 if (fsr & FSR_NVM) { 135 __quad_fdivq(x, y, &Z); 136 } else { 137 Z = *x; 138 Z.l.msw |= 0x8000; 139 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | 140 FSR_NVC; 141 __quad_setfsrp(&fsr); 142 } 143 QUAD_RETURN(Z); 144 } 145 Z = *y; 146 QUAD_RETURN(Z); 147 } 148 if (QUAD_ISNAN(*x)) { 149 if (!(x->l.msw & 0x8000)) { 150 /* snan, signal invalid */ 151 if (fsr & FSR_NVM) { 152 __quad_fdivq(x, y, &Z); 153 } else { 154 Z = *x; 155 Z.l.msw |= 0x8000; 156 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | 157 FSR_NVC; 158 __quad_setfsrp(&fsr); 159 } 160 QUAD_RETURN(Z); 161 } 162 Z = *x; 163 QUAD_RETURN(Z); 164 } 165 if (xm == 0x7fff0000) { 166 /* x is inf */ 167 if (ym == 0x7fff0000) { 168 /* inf / inf, signal invalid */ 169 if (fsr & FSR_NVM) { 170 __quad_fdivq(x, y, &Z); 171 } else { 172 Z.l.msw = 0x7fffffff; 173 Z.l.frac2 = Z.l.frac3 = 174 Z.l.frac4 = 0xffffffff; 175 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | 176 FSR_NVC; 177 __quad_setfsrp(&fsr); 178 } 179 QUAD_RETURN(Z); 180 } 181 /* inf / finite, return inf */ 182 Z.l.msw = sign | 0x7fff0000; 183 Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0; 184 QUAD_RETURN(Z); 185 } 186 /* y is inf */ 187 Z.l.msw = sign; 188 Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0; 189 QUAD_RETURN(Z); 190 } 191 192 /* handle zero cases */ 193 if (xm == 0 || ym == 0) { 194 if (QUAD_ISZERO(*x)) { 195 if (QUAD_ISZERO(*y)) { 196 /* zero / zero, signal invalid */ 197 if (fsr & FSR_NVM) { 198 __quad_fdivq(x, y, &Z); 199 } else { 200 Z.l.msw = 0x7fffffff; 201 Z.l.frac2 = Z.l.frac3 = 202 Z.l.frac4 = 0xffffffff; 203 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | 204 FSR_NVC; 205 __quad_setfsrp(&fsr); 206 } 207 QUAD_RETURN(Z); 208 } 209 /* zero / nonzero, return zero */ 210 Z.l.msw = sign; 211 Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0; 212 QUAD_RETURN(Z); 213 } 214 if (QUAD_ISZERO(*y)) { 215 /* nonzero / zero, signal zero divide */ 216 if (fsr & FSR_DZM) { 217 __quad_fdivq(x, y, &Z); 218 } else { 219 Z.l.msw = sign | 0x7fff0000; 220 Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0; 221 fsr = (fsr & ~FSR_CEXC) | FSR_DZA | FSR_DZC; 222 __quad_setfsrp(&fsr); 223 } 224 QUAD_RETURN(Z); 225 } 226 } 227 228 /* now x and y are finite, nonzero */ 229 __quad_setfsrp(&fsr_re); 230 231 /* get their normalized significands and exponents */ 232 ex = (int)(xm >> 16); 233 lx = xm & 0xffff; 234 if (ex) { 235 lx |= 0x10000; 236 wx[0] = x->l.frac2; 237 wx[1] = x->l.frac3; 238 wx[2] = x->l.frac4; 239 } else { 240 if (lx | (x->l.frac2 & 0xfffe0000)) { 241 wx[0] = x->l.frac2; 242 wx[1] = x->l.frac3; 243 wx[2] = x->l.frac4; 244 ex = 1; 245 } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) { 246 lx = x->l.frac2; 247 wx[0] = x->l.frac3; 248 wx[1] = x->l.frac4; 249 wx[2] = 0; 250 ex = -31; 251 } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) { 252 lx = x->l.frac3; 253 wx[0] = x->l.frac4; 254 wx[1] = wx[2] = 0; 255 ex = -63; 256 } else { 257 lx = x->l.frac4; 258 wx[0] = wx[1] = wx[2] = 0; 259 ex = -95; 260 } 261 while ((lx & 0x10000) == 0) { 262 lx = (lx << 1) | (wx[0] >> 31); 263 wx[0] = (wx[0] << 1) | (wx[1] >> 31); 264 wx[1] = (wx[1] << 1) | (wx[2] >> 31); 265 wx[2] <<= 1; 266 ex--; 267 } 268 } 269 ez = ex; 270 271 ey = (int)(ym >> 16); 272 ly = ym & 0xffff; 273 if (ey) { 274 ly |= 0x10000; 275 wy[0] = y->l.frac2; 276 wy[1] = y->l.frac3; 277 wy[2] = y->l.frac4; 278 } else { 279 if (ly | (y->l.frac2 & 0xfffe0000)) { 280 wy[0] = y->l.frac2; 281 wy[1] = y->l.frac3; 282 wy[2] = y->l.frac4; 283 ey = 1; 284 } else if (y->l.frac2 | (y->l.frac3 & 0xfffe0000)) { 285 ly = y->l.frac2; 286 wy[0] = y->l.frac3; 287 wy[1] = y->l.frac4; 288 wy[2] = 0; 289 ey = -31; 290 } else if (y->l.frac3 | (y->l.frac4 & 0xfffe0000)) { 291 ly = y->l.frac3; 292 wy[0] = y->l.frac4; 293 wy[1] = wy[2] = 0; 294 ey = -63; 295 } else { 296 ly = y->l.frac4; 297 wy[0] = wy[1] = wy[2] = 0; 298 ey = -95; 299 } 300 while ((ly & 0x10000) == 0) { 301 ly = (ly << 1) | (wy[0] >> 31); 302 wy[0] = (wy[0] << 1) | (wy[1] >> 31); 303 wy[1] = (wy[1] << 1) | (wy[2] >> 31); 304 wy[2] <<= 1; 305 ey--; 306 } 307 } 308 ez -= ey - 0x3fff; 309 310 /* extract the significands into doubles */ 311 c = twom16; 312 xx[0] = (double)((int)lx) * c; 313 yy[0] = (double)((int)ly) * c; 314 315 c *= twom24; 316 xx[0] += (double)((int)(wx[0] >> 8)) * c; 317 yy[1] = (double)((int)(wy[0] >> 8)) * c; 318 319 c *= twom24; 320 xx[1] = (double)((int)(((wx[0] << 16) | 321 (wx[1] >> 16)) & 0xffffff)) * c; 322 yy[2] = (double)((int)(((wy[0] << 16) | 323 (wy[1] >> 16)) & 0xffffff)) * c; 324 325 c *= twom24; 326 xx[2] = (double)((int)(((wx[1] << 8) | 327 (wx[2] >> 24)) & 0xffffff)) * c; 328 yy[3] = (double)((int)(((wy[1] << 8) | 329 (wy[2] >> 24)) & 0xffffff)) * c; 330 331 c *= twom24; 332 xx[3] = (double)((int)(wx[2] & 0xffffff)) * c; 333 yy[4] = (double)((int)(wy[2] & 0xffffff)) * c; 334 335 /* approximate the reciprocal of y */ 336 ry = one / ((yy[0] + yy[1]) + yy[2]); 337 338 /* compute the first five "digits" of the quotient */ 339 zz[0] = (ry * (xx[0] + xx[1]) + three2p27) - three2p27; 340 xx[0] = ((xx[0] - zz[0] * yy[0]) - zz[0] * yy[1]) + xx[1]; 341 d = zz[0] * yy[2]; 342 c = (d + three2m13) - three2m13; 343 xx[0] -= c; 344 xx[1] = xx[2] - (d - c); 345 d = zz[0] * yy[3]; 346 c = (d + three2m37) - three2m37; 347 xx[1] -= c; 348 xx[2] = xx[3] - (d - c); 349 d = zz[0] * yy[4]; 350 c = (d + three2m61) - three2m61; 351 xx[2] -= c; 352 xx[3] = c - d; 353 354 zz[1] = (ry * (xx[0] + xx[1]) + three2p3) - three2p3; 355 xx[0] = ((xx[0] - zz[1] * yy[0]) - zz[1] * yy[1]) + xx[1]; 356 d = zz[1] * yy[2]; 357 c = (d + three2m37) - three2m37; 358 xx[0] -= c; 359 xx[1] = xx[2] - (d - c); 360 d = zz[1] * yy[3]; 361 c = (d + three2m61) - three2m61; 362 xx[1] -= c; 363 xx[2] = xx[3] - (d - c); 364 d = zz[1] * yy[4]; 365 c = (d + three2m85) - three2m85; 366 xx[2] -= c; 367 xx[3] = c - d; 368 369 zz[2] = (ry * (xx[0] + xx[1]) + three2m21) - three2m21; 370 xx[0] = ((xx[0] - zz[2] * yy[0]) - zz[2] * yy[1]) + xx[1]; 371 d = zz[2] * yy[2]; 372 c = (d + three2m61) - three2m61; 373 xx[0] -= c; 374 xx[1] = xx[2] - (d - c); 375 d = zz[2] * yy[3]; 376 c = (d + three2m85) - three2m85; 377 xx[1] -= c; 378 xx[2] = xx[3] - (d - c); 379 d = zz[2] * yy[4]; 380 c = (d + three2m109) - three2m109; 381 xx[2] -= c; 382 xx[3] = c - d; 383 384 zz[3] = (ry * (xx[0] + xx[1]) + three2m45) - three2m45; 385 xx[0] = ((xx[0] - zz[3] * yy[0]) - zz[3] * yy[1]) + xx[1]; 386 d = zz[3] * yy[2]; 387 c = (d + three2m85) - three2m85; 388 xx[0] -= c; 389 xx[1] = xx[2] - (d - c); 390 d = zz[3] * yy[3]; 391 c = (d + three2m109) - three2m109; 392 xx[1] -= c; 393 xx[2] = xx[3] - (d - c); 394 d = zz[3] * yy[4]; 395 c = (d + three2m133) - three2m133; 396 xx[2] -= c; 397 xx[3] = c - d; 398 399 zz[4] = (ry * (xx[0] + xx[1]) + three2m73) - three2m73; 400 401 /* reduce to three doubles, making sure zz[1] is positive */ 402 zz[0] += zz[1] - twom48; 403 zz[1] = twom48 + zz[2] + zz[3]; 404 zz[2] = zz[4]; 405 406 /* if the third term might lie on a rounding boundary, perturb it */ 407 if (zz[2] == (twom62 + zz[2]) - twom62) { 408 /* here we just need to get the sign of the remainder */ 409 c = (((((xx[0] - zz[4] * yy[0]) - zz[4] * yy[1]) + xx[1]) + 410 (xx[2] - zz[4] * yy[2])) + (xx[3] - zz[4] * yy[3])) 411 - zz[4] * yy[4]; 412 if (c < zero) 413 zz[2] -= twom124; 414 else if (c > zero) 415 zz[2] += twom124; 416 } 417 418 /* 419 * propagate carries/borrows, using round-to-negative-infinity mode 420 * to make all terms nonnegative (note that we can't encounter a 421 * borrow so large that the roundoff is unrepresentable because 422 * we took care to make zz[1] positive above) 423 */ 424 __quad_setfsrp(&fsr_rn); 425 c = zz[1] + zz[2]; 426 zz[2] += (zz[1] - c); 427 zz[1] = c; 428 c = zz[0] + zz[1]; 429 zz[1] += (zz[0] - c); 430 zz[0] = c; 431 432 /* check for borrow */ 433 if (c < one) { 434 /* postnormalize */ 435 zz[0] += zz[0]; 436 zz[1] += zz[1]; 437 zz[2] += zz[2]; 438 ez--; 439 } 440 441 /* if exponent > 0 strip off integer bit, else denormalize */ 442 if (ez > 0) { 443 ibit = 1; 444 zz[0] -= one; 445 } else { 446 ibit = 0; 447 if (ez > -128) 448 u.l.hi = (unsigned int)(0x3fe + ez) << 20; 449 else 450 u.l.hi = 0x37e00000; 451 u.l.lo = 0; 452 zz[0] *= u.d; 453 zz[1] *= u.d; 454 zz[2] *= u.d; 455 ez = 0; 456 } 457 458 /* the first 48 bits of fraction come from zz[0] */ 459 u.d = d = two36 + zz[0]; 460 msw = u.l.lo; 461 zz[0] -= (d - two36); 462 463 u.d = d = two4 + zz[0]; 464 frac2 = u.l.lo; 465 zz[0] -= (d - two4); 466 467 /* the next 32 come from zz[0] and zz[1] */ 468 u.d = d = twom28 + (zz[0] + zz[1]); 469 frac3 = u.l.lo; 470 zz[0] -= (d - twom28); 471 472 /* condense the remaining fraction; errors here won't matter */ 473 c = zz[0] + zz[1]; 474 zz[1] = ((zz[0] - c) + zz[1]) + zz[2]; 475 zz[0] = c; 476 477 /* get the last word of fraction */ 478 u.d = d = twom60 + (zz[0] + zz[1]); 479 frac4 = u.l.lo; 480 zz[0] -= (d - twom60); 481 482 /* keep track of what's left for rounding; note that the error */ 483 /* in computing c will be non-negative due to rounding mode */ 484 c = zz[0] + zz[1]; 485 486 /* get the rounding mode, fudging directed rounding modes */ 487 /* as though the result were positive */ 488 rm = fsr >> 30; 489 if (sign) 490 rm ^= (rm >> 1); 491 492 /* round and raise exceptions */ 493 fsr &= ~FSR_CEXC; 494 if (c != zero) { 495 fsr |= FSR_NXC; 496 497 /* decide whether to round the fraction up */ 498 if (rm == FSR_RP || (rm == FSR_RN && (c > twom113 || 499 (c == twom113 && ((frac4 & 1) || (c - zz[0] != 500 zz[1])))))) { 501 /* round up and renormalize if necessary */ 502 if (++frac4 == 0) 503 if (++frac3 == 0) 504 if (++frac2 == 0) 505 if (++msw == 0x10000) { 506 msw = 0; 507 ez++; 508 } 509 } 510 } 511 512 /* check for under/overflow */ 513 if (ez >= 0x7fff) { 514 if (rm == FSR_RN || rm == FSR_RP) { 515 z.l.msw = sign | 0x7fff0000; 516 z.l.frac2 = z.l.frac3 = z.l.frac4 = 0; 517 } else { 518 z.l.msw = sign | 0x7ffeffff; 519 z.l.frac2 = z.l.frac3 = z.l.frac4 = 0xffffffff; 520 } 521 fsr |= FSR_OFC | FSR_NXC; 522 } else { 523 z.l.msw = sign | (ez << 16) | msw; 524 z.l.frac2 = frac2; 525 z.l.frac3 = frac3; 526 z.l.frac4 = frac4; 527 528 /* !ibit => exact result was tiny before rounding, */ 529 /* t nonzero => result delivered is inexact */ 530 if (!ibit) { 531 if (c != zero) 532 fsr |= FSR_UFC | FSR_NXC; 533 else if (fsr & FSR_UFM) 534 fsr |= FSR_UFC; 535 } 536 } 537 538 if ((fsr & FSR_CEXC) & (fsr >> 23)) { 539 __quad_setfsrp(&fsr); 540 __quad_fdivq(x, y, &Z); 541 } else { 542 Z = z; 543 fsr |= (fsr & 0x1f) << 5; 544 __quad_setfsrp(&fsr); 545 } 546 QUAD_RETURN(Z); 547 } 548