1 /* Software floating-point emulation. Common operations. 2 Copyright (C) 1997,1998,1999,2006,2007 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 Contributed by Richard Henderson (rth@cygnus.com), 5 Jakub Jelinek (jj@ultra.linux.cz), 6 David S. Miller (davem@redhat.com) and 7 Peter Maydell (pmaydell@chiark.greenend.org.uk). 8 9 The GNU C Library is free software; you can redistribute it and/or 10 modify it under the terms of the GNU Lesser General Public 11 License as published by the Free Software Foundation; either 12 version 2.1 of the License, or (at your option) any later version. 13 14 In addition to the permissions in the GNU Lesser General Public 15 License, the Free Software Foundation gives you unlimited 16 permission to link the compiled version of this file into 17 combinations with other programs, and to distribute those 18 combinations without any restriction coming from the use of this 19 file. (The Lesser General Public License restrictions do apply in 20 other respects; for example, they cover modification of the file, 21 and distribution when not linked into a combine executable.) 22 23 The GNU C Library is distributed in the hope that it will be useful, 24 but WITHOUT ANY WARRANTY; without even the implied warranty of 25 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 26 Lesser General Public License for more details. 27 28 You should have received a copy of the GNU Lesser General Public 29 License along with the GNU C Library; if not, see 30 <http://www.gnu.org/licenses/>. */ 31 32 #define _FP_DECL(wc, X) \ 33 _FP_I_TYPE X##_c __attribute__((unused)), X##_s, X##_e; \ 34 _FP_FRAC_DECL_##wc(X) 35 36 /* 37 * Finish truely unpacking a native fp value by classifying the kind 38 * of fp value and normalizing both the exponent and the fraction. 39 */ 40 41 #define _FP_UNPACK_CANONICAL(fs, wc, X) \ 42 do { \ 43 switch (X##_e) \ 44 { \ 45 default: \ 46 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \ 47 _FP_FRAC_SLL_##wc(X, _FP_WORKBITS); \ 48 X##_e -= _FP_EXPBIAS_##fs; \ 49 X##_c = FP_CLS_NORMAL; \ 50 break; \ 51 \ 52 case 0: \ 53 if (_FP_FRAC_ZEROP_##wc(X)) \ 54 X##_c = FP_CLS_ZERO; \ 55 else \ 56 { \ 57 /* a denormalized number */ \ 58 _FP_I_TYPE _shift; \ 59 _FP_FRAC_CLZ_##wc(_shift, X); \ 60 _shift -= _FP_FRACXBITS_##fs; \ 61 _FP_FRAC_SLL_##wc(X, (_shift+_FP_WORKBITS)); \ 62 X##_e -= _FP_EXPBIAS_##fs - 1 + _shift; \ 63 X##_c = FP_CLS_NORMAL; \ 64 FP_SET_EXCEPTION(FP_EX_DENORM); \ 65 } \ 66 break; \ 67 \ 68 case _FP_EXPMAX_##fs: \ 69 if (_FP_FRAC_ZEROP_##wc(X)) \ 70 X##_c = FP_CLS_INF; \ 71 else \ 72 { \ 73 X##_c = FP_CLS_NAN; \ 74 /* Check for signaling NaN */ \ 75 if (!(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \ 76 FP_SET_EXCEPTION(FP_EX_INVALID); \ 77 } \ 78 break; \ 79 } \ 80 } while (0) 81 82 /* Finish unpacking an fp value in semi-raw mode: the mantissa is 83 shifted by _FP_WORKBITS but the implicit MSB is not inserted and 84 other classification is not done. */ 85 #define _FP_UNPACK_SEMIRAW(fs, wc, X) _FP_FRAC_SLL_##wc(X, _FP_WORKBITS) 86 87 /* A semi-raw value has overflowed to infinity. Adjust the mantissa 88 and exponent appropriately. */ 89 #define _FP_OVERFLOW_SEMIRAW(fs, wc, X) \ 90 do { \ 91 if (FP_ROUNDMODE == FP_RND_NEAREST \ 92 || (FP_ROUNDMODE == FP_RND_PINF && !X##_s) \ 93 || (FP_ROUNDMODE == FP_RND_MINF && X##_s)) \ 94 { \ 95 X##_e = _FP_EXPMAX_##fs; \ 96 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \ 97 } \ 98 else \ 99 { \ 100 X##_e = _FP_EXPMAX_##fs - 1; \ 101 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \ 102 } \ 103 FP_SET_EXCEPTION(FP_EX_INEXACT); \ 104 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \ 105 } while (0) 106 107 /* Check for a semi-raw value being a signaling NaN and raise the 108 invalid exception if so. */ 109 #define _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X) \ 110 do { \ 111 if (X##_e == _FP_EXPMAX_##fs \ 112 && !_FP_FRAC_ZEROP_##wc(X) \ 113 && !(_FP_FRAC_HIGH_##fs(X) & _FP_QNANBIT_SH_##fs)) \ 114 FP_SET_EXCEPTION(FP_EX_INVALID); \ 115 } while (0) 116 117 /* Choose a NaN result from an operation on two semi-raw NaN 118 values. */ 119 #define _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP) \ 120 do { \ 121 /* _FP_CHOOSENAN expects raw values, so shift as required. */ \ 122 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \ 123 _FP_FRAC_SRL_##wc(Y, _FP_WORKBITS); \ 124 _FP_CHOOSENAN(fs, wc, R, X, Y, OP); \ 125 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \ 126 } while (0) 127 128 /* Test whether a biased exponent is normal (not zero or maximum). */ 129 #define _FP_EXP_NORMAL(fs, wc, X) (((X##_e + 1) & _FP_EXPMAX_##fs) > 1) 130 131 /* Prepare to pack an fp value in semi-raw mode: the mantissa is 132 rounded and shifted right, with the rounding possibly increasing 133 the exponent (including changing a finite value to infinity). */ 134 #define _FP_PACK_SEMIRAW(fs, wc, X) \ 135 do { \ 136 _FP_ROUND(wc, X); \ 137 if (_FP_FRAC_HIGH_##fs(X) \ 138 & (_FP_OVERFLOW_##fs >> 1)) \ 139 { \ 140 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_OVERFLOW_##fs >> 1); \ 141 X##_e++; \ 142 if (X##_e == _FP_EXPMAX_##fs) \ 143 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \ 144 } \ 145 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \ 146 if (!_FP_EXP_NORMAL(fs, wc, X) && !_FP_FRAC_ZEROP_##wc(X)) \ 147 { \ 148 if (X##_e == 0) \ 149 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \ 150 else \ 151 { \ 152 if (!_FP_KEEPNANFRACP) \ 153 { \ 154 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \ 155 X##_s = _FP_NANSIGN_##fs; \ 156 } \ 157 else \ 158 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \ 159 } \ 160 } \ 161 } while (0) 162 163 /* 164 * Before packing the bits back into the native fp result, take care 165 * of such mundane things as rounding and overflow. Also, for some 166 * kinds of fp values, the original parts may not have been fully 167 * extracted -- but that is ok, we can regenerate them now. 168 */ 169 170 #define _FP_PACK_CANONICAL(fs, wc, X) \ 171 do { \ 172 switch (X##_c) \ 173 { \ 174 case FP_CLS_NORMAL: \ 175 X##_e += _FP_EXPBIAS_##fs; \ 176 if (X##_e > 0) \ 177 { \ 178 _FP_ROUND(wc, X); \ 179 if (_FP_FRAC_OVERP_##wc(fs, X)) \ 180 { \ 181 _FP_FRAC_CLEAR_OVERP_##wc(fs, X); \ 182 X##_e++; \ 183 } \ 184 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \ 185 if (X##_e >= _FP_EXPMAX_##fs) \ 186 { \ 187 /* overflow */ \ 188 switch (FP_ROUNDMODE) \ 189 { \ 190 case FP_RND_NEAREST: \ 191 X##_c = FP_CLS_INF; \ 192 break; \ 193 case FP_RND_PINF: \ 194 if (!X##_s) X##_c = FP_CLS_INF; \ 195 break; \ 196 case FP_RND_MINF: \ 197 if (X##_s) X##_c = FP_CLS_INF; \ 198 break; \ 199 } \ 200 if (X##_c == FP_CLS_INF) \ 201 { \ 202 /* Overflow to infinity */ \ 203 X##_e = _FP_EXPMAX_##fs; \ 204 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \ 205 } \ 206 else \ 207 { \ 208 /* Overflow to maximum normal */ \ 209 X##_e = _FP_EXPMAX_##fs - 1; \ 210 _FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc); \ 211 } \ 212 FP_SET_EXCEPTION(FP_EX_OVERFLOW); \ 213 FP_SET_EXCEPTION(FP_EX_INEXACT); \ 214 } \ 215 } \ 216 else \ 217 { \ 218 /* we've got a denormalized number */ \ 219 X##_e = -X##_e + 1; \ 220 if (X##_e <= _FP_WFRACBITS_##fs) \ 221 { \ 222 _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs); \ 223 _FP_ROUND(wc, X); \ 224 if (_FP_FRAC_HIGH_##fs(X) \ 225 & (_FP_OVERFLOW_##fs >> 1)) \ 226 { \ 227 X##_e = 1; \ 228 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \ 229 } \ 230 else \ 231 { \ 232 X##_e = 0; \ 233 _FP_FRAC_SRL_##wc(X, _FP_WORKBITS); \ 234 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \ 235 } \ 236 } \ 237 else \ 238 { \ 239 /* underflow to zero */ \ 240 X##_e = 0; \ 241 if (!_FP_FRAC_ZEROP_##wc(X)) \ 242 { \ 243 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \ 244 _FP_ROUND(wc, X); \ 245 _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS); \ 246 } \ 247 FP_SET_EXCEPTION(FP_EX_UNDERFLOW); \ 248 } \ 249 } \ 250 break; \ 251 \ 252 case FP_CLS_ZERO: \ 253 X##_e = 0; \ 254 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \ 255 break; \ 256 \ 257 case FP_CLS_INF: \ 258 X##_e = _FP_EXPMAX_##fs; \ 259 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \ 260 break; \ 261 \ 262 case FP_CLS_NAN: \ 263 X##_e = _FP_EXPMAX_##fs; \ 264 if (!_FP_KEEPNANFRACP) \ 265 { \ 266 _FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs); \ 267 X##_s = _FP_NANSIGN_##fs; \ 268 } \ 269 else \ 270 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs; \ 271 break; \ 272 } \ 273 } while (0) 274 275 /* This one accepts raw argument and not cooked, returns 276 * 1 if X is a signaling NaN. 277 */ 278 #define _FP_ISSIGNAN(fs, wc, X) \ 279 ({ \ 280 int __ret = 0; \ 281 if (X##_e == _FP_EXPMAX_##fs) \ 282 { \ 283 if (!_FP_FRAC_ZEROP_##wc(X) \ 284 && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs)) \ 285 __ret = 1; \ 286 } \ 287 __ret; \ 288 }) 289 290 291 292 293 294 /* Addition on semi-raw values. */ 295 #define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP) \ 296 do { \ 297 if (X##_s == Y##_s) \ 298 { \ 299 /* Addition. */ \ 300 R##_s = X##_s; \ 301 int ediff = X##_e - Y##_e; \ 302 if (ediff > 0) \ 303 { \ 304 R##_e = X##_e; \ 305 if (Y##_e == 0) \ 306 { \ 307 /* Y is zero or denormalized. */ \ 308 if (_FP_FRAC_ZEROP_##wc(Y)) \ 309 { \ 310 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 311 _FP_FRAC_COPY_##wc(R, X); \ 312 goto add_done; \ 313 } \ 314 else \ 315 { \ 316 FP_SET_EXCEPTION(FP_EX_DENORM); \ 317 ediff--; \ 318 if (ediff == 0) \ 319 { \ 320 _FP_FRAC_ADD_##wc(R, X, Y); \ 321 goto add3; \ 322 } \ 323 if (X##_e == _FP_EXPMAX_##fs) \ 324 { \ 325 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 326 _FP_FRAC_COPY_##wc(R, X); \ 327 goto add_done; \ 328 } \ 329 goto add1; \ 330 } \ 331 } \ 332 else if (X##_e == _FP_EXPMAX_##fs) \ 333 { \ 334 /* X is NaN or Inf, Y is normal. */ \ 335 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 336 _FP_FRAC_COPY_##wc(R, X); \ 337 goto add_done; \ 338 } \ 339 \ 340 /* Insert implicit MSB of Y. */ \ 341 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \ 342 \ 343 add1: \ 344 /* Shift the mantissa of Y to the right EDIFF steps; \ 345 remember to account later for the implicit MSB of X. */ \ 346 if (ediff <= _FP_WFRACBITS_##fs) \ 347 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \ 348 else if (!_FP_FRAC_ZEROP_##wc(Y)) \ 349 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \ 350 _FP_FRAC_ADD_##wc(R, X, Y); \ 351 } \ 352 else if (ediff < 0) \ 353 { \ 354 ediff = -ediff; \ 355 R##_e = Y##_e; \ 356 if (X##_e == 0) \ 357 { \ 358 /* X is zero or denormalized. */ \ 359 if (_FP_FRAC_ZEROP_##wc(X)) \ 360 { \ 361 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 362 _FP_FRAC_COPY_##wc(R, Y); \ 363 goto add_done; \ 364 } \ 365 else \ 366 { \ 367 FP_SET_EXCEPTION(FP_EX_DENORM); \ 368 ediff--; \ 369 if (ediff == 0) \ 370 { \ 371 _FP_FRAC_ADD_##wc(R, Y, X); \ 372 goto add3; \ 373 } \ 374 if (Y##_e == _FP_EXPMAX_##fs) \ 375 { \ 376 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 377 _FP_FRAC_COPY_##wc(R, Y); \ 378 goto add_done; \ 379 } \ 380 goto add2; \ 381 } \ 382 } \ 383 else if (Y##_e == _FP_EXPMAX_##fs) \ 384 { \ 385 /* Y is NaN or Inf, X is normal. */ \ 386 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 387 _FP_FRAC_COPY_##wc(R, Y); \ 388 goto add_done; \ 389 } \ 390 \ 391 /* Insert implicit MSB of X. */ \ 392 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \ 393 \ 394 add2: \ 395 /* Shift the mantissa of X to the right EDIFF steps; \ 396 remember to account later for the implicit MSB of Y. */ \ 397 if (ediff <= _FP_WFRACBITS_##fs) \ 398 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \ 399 else if (!_FP_FRAC_ZEROP_##wc(X)) \ 400 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \ 401 _FP_FRAC_ADD_##wc(R, Y, X); \ 402 } \ 403 else \ 404 { \ 405 /* ediff == 0. */ \ 406 if (!_FP_EXP_NORMAL(fs, wc, X)) \ 407 { \ 408 if (X##_e == 0) \ 409 { \ 410 /* X and Y are zero or denormalized. */ \ 411 R##_e = 0; \ 412 if (_FP_FRAC_ZEROP_##wc(X)) \ 413 { \ 414 if (!_FP_FRAC_ZEROP_##wc(Y)) \ 415 FP_SET_EXCEPTION(FP_EX_DENORM); \ 416 _FP_FRAC_COPY_##wc(R, Y); \ 417 goto add_done; \ 418 } \ 419 else if (_FP_FRAC_ZEROP_##wc(Y)) \ 420 { \ 421 FP_SET_EXCEPTION(FP_EX_DENORM); \ 422 _FP_FRAC_COPY_##wc(R, X); \ 423 goto add_done; \ 424 } \ 425 else \ 426 { \ 427 FP_SET_EXCEPTION(FP_EX_DENORM); \ 428 _FP_FRAC_ADD_##wc(R, X, Y); \ 429 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \ 430 { \ 431 /* Normalized result. */ \ 432 _FP_FRAC_HIGH_##fs(R) \ 433 &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \ 434 R##_e = 1; \ 435 } \ 436 goto add_done; \ 437 } \ 438 } \ 439 else \ 440 { \ 441 /* X and Y are NaN or Inf. */ \ 442 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 443 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 444 R##_e = _FP_EXPMAX_##fs; \ 445 if (_FP_FRAC_ZEROP_##wc(X)) \ 446 _FP_FRAC_COPY_##wc(R, Y); \ 447 else if (_FP_FRAC_ZEROP_##wc(Y)) \ 448 _FP_FRAC_COPY_##wc(R, X); \ 449 else \ 450 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \ 451 goto add_done; \ 452 } \ 453 } \ 454 /* The exponents of X and Y, both normal, are equal. The \ 455 implicit MSBs will always add to increase the \ 456 exponent. */ \ 457 _FP_FRAC_ADD_##wc(R, X, Y); \ 458 R##_e = X##_e + 1; \ 459 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \ 460 if (R##_e == _FP_EXPMAX_##fs) \ 461 /* Overflow to infinity (depending on rounding mode). */ \ 462 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \ 463 goto add_done; \ 464 } \ 465 add3: \ 466 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \ 467 { \ 468 /* Overflow. */ \ 469 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \ 470 R##_e++; \ 471 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \ 472 if (R##_e == _FP_EXPMAX_##fs) \ 473 /* Overflow to infinity (depending on rounding mode). */ \ 474 _FP_OVERFLOW_SEMIRAW(fs, wc, R); \ 475 } \ 476 add_done: ; \ 477 } \ 478 else \ 479 { \ 480 /* Subtraction. */ \ 481 int ediff = X##_e - Y##_e; \ 482 if (ediff > 0) \ 483 { \ 484 R##_e = X##_e; \ 485 R##_s = X##_s; \ 486 if (Y##_e == 0) \ 487 { \ 488 /* Y is zero or denormalized. */ \ 489 if (_FP_FRAC_ZEROP_##wc(Y)) \ 490 { \ 491 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 492 _FP_FRAC_COPY_##wc(R, X); \ 493 goto sub_done; \ 494 } \ 495 else \ 496 { \ 497 FP_SET_EXCEPTION(FP_EX_DENORM); \ 498 ediff--; \ 499 if (ediff == 0) \ 500 { \ 501 _FP_FRAC_SUB_##wc(R, X, Y); \ 502 goto sub3; \ 503 } \ 504 if (X##_e == _FP_EXPMAX_##fs) \ 505 { \ 506 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 507 _FP_FRAC_COPY_##wc(R, X); \ 508 goto sub_done; \ 509 } \ 510 goto sub1; \ 511 } \ 512 } \ 513 else if (X##_e == _FP_EXPMAX_##fs) \ 514 { \ 515 /* X is NaN or Inf, Y is normal. */ \ 516 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 517 _FP_FRAC_COPY_##wc(R, X); \ 518 goto sub_done; \ 519 } \ 520 \ 521 /* Insert implicit MSB of Y. */ \ 522 _FP_FRAC_HIGH_##fs(Y) |= _FP_IMPLBIT_SH_##fs; \ 523 \ 524 sub1: \ 525 /* Shift the mantissa of Y to the right EDIFF steps; \ 526 remember to account later for the implicit MSB of X. */ \ 527 if (ediff <= _FP_WFRACBITS_##fs) \ 528 _FP_FRAC_SRS_##wc(Y, ediff, _FP_WFRACBITS_##fs); \ 529 else if (!_FP_FRAC_ZEROP_##wc(Y)) \ 530 _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc); \ 531 _FP_FRAC_SUB_##wc(R, X, Y); \ 532 } \ 533 else if (ediff < 0) \ 534 { \ 535 ediff = -ediff; \ 536 R##_e = Y##_e; \ 537 R##_s = Y##_s; \ 538 if (X##_e == 0) \ 539 { \ 540 /* X is zero or denormalized. */ \ 541 if (_FP_FRAC_ZEROP_##wc(X)) \ 542 { \ 543 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 544 _FP_FRAC_COPY_##wc(R, Y); \ 545 goto sub_done; \ 546 } \ 547 else \ 548 { \ 549 FP_SET_EXCEPTION(FP_EX_DENORM); \ 550 ediff--; \ 551 if (ediff == 0) \ 552 { \ 553 _FP_FRAC_SUB_##wc(R, Y, X); \ 554 goto sub3; \ 555 } \ 556 if (Y##_e == _FP_EXPMAX_##fs) \ 557 { \ 558 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 559 _FP_FRAC_COPY_##wc(R, Y); \ 560 goto sub_done; \ 561 } \ 562 goto sub2; \ 563 } \ 564 } \ 565 else if (Y##_e == _FP_EXPMAX_##fs) \ 566 { \ 567 /* Y is NaN or Inf, X is normal. */ \ 568 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 569 _FP_FRAC_COPY_##wc(R, Y); \ 570 goto sub_done; \ 571 } \ 572 \ 573 /* Insert implicit MSB of X. */ \ 574 _FP_FRAC_HIGH_##fs(X) |= _FP_IMPLBIT_SH_##fs; \ 575 \ 576 sub2: \ 577 /* Shift the mantissa of X to the right EDIFF steps; \ 578 remember to account later for the implicit MSB of Y. */ \ 579 if (ediff <= _FP_WFRACBITS_##fs) \ 580 _FP_FRAC_SRS_##wc(X, ediff, _FP_WFRACBITS_##fs); \ 581 else if (!_FP_FRAC_ZEROP_##wc(X)) \ 582 _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc); \ 583 _FP_FRAC_SUB_##wc(R, Y, X); \ 584 } \ 585 else \ 586 { \ 587 /* ediff == 0. */ \ 588 if (!_FP_EXP_NORMAL(fs, wc, X)) \ 589 { \ 590 if (X##_e == 0) \ 591 { \ 592 /* X and Y are zero or denormalized. */ \ 593 R##_e = 0; \ 594 if (_FP_FRAC_ZEROP_##wc(X)) \ 595 { \ 596 _FP_FRAC_COPY_##wc(R, Y); \ 597 if (_FP_FRAC_ZEROP_##wc(Y)) \ 598 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \ 599 else \ 600 { \ 601 FP_SET_EXCEPTION(FP_EX_DENORM); \ 602 R##_s = Y##_s; \ 603 } \ 604 goto sub_done; \ 605 } \ 606 else if (_FP_FRAC_ZEROP_##wc(Y)) \ 607 { \ 608 FP_SET_EXCEPTION(FP_EX_DENORM); \ 609 _FP_FRAC_COPY_##wc(R, X); \ 610 R##_s = X##_s; \ 611 goto sub_done; \ 612 } \ 613 else \ 614 { \ 615 FP_SET_EXCEPTION(FP_EX_DENORM); \ 616 _FP_FRAC_SUB_##wc(R, X, Y); \ 617 R##_s = X##_s; \ 618 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \ 619 { \ 620 /* |X| < |Y|, negate result. */ \ 621 _FP_FRAC_SUB_##wc(R, Y, X); \ 622 R##_s = Y##_s; \ 623 } \ 624 else if (_FP_FRAC_ZEROP_##wc(R)) \ 625 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \ 626 goto sub_done; \ 627 } \ 628 } \ 629 else \ 630 { \ 631 /* X and Y are NaN or Inf, of opposite signs. */ \ 632 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, X); \ 633 _FP_CHECK_SIGNAN_SEMIRAW(fs, wc, Y); \ 634 R##_e = _FP_EXPMAX_##fs; \ 635 if (_FP_FRAC_ZEROP_##wc(X)) \ 636 { \ 637 if (_FP_FRAC_ZEROP_##wc(Y)) \ 638 { \ 639 /* Inf - Inf. */ \ 640 R##_s = _FP_NANSIGN_##fs; \ 641 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \ 642 _FP_FRAC_SLL_##wc(R, _FP_WORKBITS); \ 643 FP_SET_EXCEPTION(FP_EX_INVALID); \ 644 } \ 645 else \ 646 { \ 647 /* Inf - NaN. */ \ 648 R##_s = Y##_s; \ 649 _FP_FRAC_COPY_##wc(R, Y); \ 650 } \ 651 } \ 652 else \ 653 { \ 654 if (_FP_FRAC_ZEROP_##wc(Y)) \ 655 { \ 656 /* NaN - Inf. */ \ 657 R##_s = X##_s; \ 658 _FP_FRAC_COPY_##wc(R, X); \ 659 } \ 660 else \ 661 { \ 662 /* NaN - NaN. */ \ 663 _FP_CHOOSENAN_SEMIRAW(fs, wc, R, X, Y, OP); \ 664 } \ 665 } \ 666 goto sub_done; \ 667 } \ 668 } \ 669 /* The exponents of X and Y, both normal, are equal. The \ 670 implicit MSBs cancel. */ \ 671 R##_e = X##_e; \ 672 _FP_FRAC_SUB_##wc(R, X, Y); \ 673 R##_s = X##_s; \ 674 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \ 675 { \ 676 /* |X| < |Y|, negate result. */ \ 677 _FP_FRAC_SUB_##wc(R, Y, X); \ 678 R##_s = Y##_s; \ 679 } \ 680 else if (_FP_FRAC_ZEROP_##wc(R)) \ 681 { \ 682 R##_e = 0; \ 683 R##_s = (FP_ROUNDMODE == FP_RND_MINF); \ 684 goto sub_done; \ 685 } \ 686 goto norm; \ 687 } \ 688 sub3: \ 689 if (_FP_FRAC_HIGH_##fs(R) & _FP_IMPLBIT_SH_##fs) \ 690 { \ 691 int diff; \ 692 /* Carry into most significant bit of larger one of X and Y, \ 693 canceling it; renormalize. */ \ 694 _FP_FRAC_HIGH_##fs(R) &= _FP_IMPLBIT_SH_##fs - 1; \ 695 norm: \ 696 _FP_FRAC_CLZ_##wc(diff, R); \ 697 diff -= _FP_WFRACXBITS_##fs; \ 698 _FP_FRAC_SLL_##wc(R, diff); \ 699 if (R##_e <= diff) \ 700 { \ 701 /* R is denormalized. */ \ 702 diff = diff - R##_e + 1; \ 703 _FP_FRAC_SRS_##wc(R, diff, _FP_WFRACBITS_##fs); \ 704 R##_e = 0; \ 705 } \ 706 else \ 707 { \ 708 R##_e -= diff; \ 709 _FP_FRAC_HIGH_##fs(R) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \ 710 } \ 711 } \ 712 sub_done: ; \ 713 } \ 714 } while (0) 715 716 #define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+') 717 #define _FP_SUB(fs, wc, R, X, Y) \ 718 do { \ 719 if (!(Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) Y##_s ^= 1; \ 720 _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-'); \ 721 } while (0) 722 723 724 /* 725 * Main negation routine. FIXME -- when we care about setting exception 726 * bits reliably, this will not do. We should examine all of the fp classes. 727 */ 728 729 #define _FP_NEG(fs, wc, R, X) \ 730 do { \ 731 _FP_FRAC_COPY_##wc(R, X); \ 732 R##_c = X##_c; \ 733 R##_e = X##_e; \ 734 R##_s = 1 ^ X##_s; \ 735 } while (0) 736 737 738 /* 739 * Main multiplication routine. The input values should be cooked. 740 */ 741 742 #define _FP_MUL(fs, wc, R, X, Y) \ 743 do { \ 744 R##_s = X##_s ^ Y##_s; \ 745 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \ 746 { \ 747 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \ 748 R##_c = FP_CLS_NORMAL; \ 749 R##_e = X##_e + Y##_e + 1; \ 750 \ 751 _FP_MUL_MEAT_##fs(R,X,Y); \ 752 \ 753 if (_FP_FRAC_OVERP_##wc(fs, R)) \ 754 _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs); \ 755 else \ 756 R##_e--; \ 757 break; \ 758 \ 759 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \ 760 _FP_CHOOSENAN(fs, wc, R, X, Y, '*'); \ 761 break; \ 762 \ 763 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \ 764 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \ 765 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \ 766 R##_s = X##_s; \ 767 \ 768 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \ 769 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \ 770 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \ 771 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \ 772 _FP_FRAC_COPY_##wc(R, X); \ 773 R##_c = X##_c; \ 774 break; \ 775 \ 776 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \ 777 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \ 778 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \ 779 R##_s = Y##_s; \ 780 \ 781 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \ 782 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \ 783 _FP_FRAC_COPY_##wc(R, Y); \ 784 R##_c = Y##_c; \ 785 break; \ 786 \ 787 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \ 788 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \ 789 R##_s = _FP_NANSIGN_##fs; \ 790 R##_c = FP_CLS_NAN; \ 791 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \ 792 FP_SET_EXCEPTION(FP_EX_INVALID); \ 793 break; \ 794 \ 795 default: \ 796 abort(); \ 797 } \ 798 } while (0) 799 800 801 /* 802 * Main division routine. The input values should be cooked. 803 */ 804 805 #define _FP_DIV(fs, wc, R, X, Y) \ 806 do { \ 807 R##_s = X##_s ^ Y##_s; \ 808 switch (_FP_CLS_COMBINE(X##_c, Y##_c)) \ 809 { \ 810 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL): \ 811 R##_c = FP_CLS_NORMAL; \ 812 R##_e = X##_e - Y##_e; \ 813 \ 814 _FP_DIV_MEAT_##fs(R,X,Y); \ 815 break; \ 816 \ 817 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN): \ 818 _FP_CHOOSENAN(fs, wc, R, X, Y, '/'); \ 819 break; \ 820 \ 821 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL): \ 822 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF): \ 823 case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO): \ 824 R##_s = X##_s; \ 825 _FP_FRAC_COPY_##wc(R, X); \ 826 R##_c = X##_c; \ 827 break; \ 828 \ 829 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN): \ 830 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN): \ 831 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN): \ 832 R##_s = Y##_s; \ 833 _FP_FRAC_COPY_##wc(R, Y); \ 834 R##_c = Y##_c; \ 835 break; \ 836 \ 837 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF): \ 838 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF): \ 839 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL): \ 840 R##_c = FP_CLS_ZERO; \ 841 break; \ 842 \ 843 case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO): \ 844 FP_SET_EXCEPTION(FP_EX_DIVZERO); \ 845 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO): \ 846 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL): \ 847 R##_c = FP_CLS_INF; \ 848 break; \ 849 \ 850 case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF): \ 851 case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO): \ 852 R##_s = _FP_NANSIGN_##fs; \ 853 R##_c = FP_CLS_NAN; \ 854 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \ 855 FP_SET_EXCEPTION(FP_EX_INVALID); \ 856 break; \ 857 \ 858 default: \ 859 abort(); \ 860 } \ 861 } while (0) 862 863 864 /* 865 * Main differential comparison routine. The inputs should be raw not 866 * cooked. The return is -1,0,1 for normal values, 2 otherwise. 867 */ 868 869 #define _FP_CMP(fs, wc, ret, X, Y, un) \ 870 do { \ 871 /* NANs are unordered */ \ 872 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \ 873 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \ 874 { \ 875 ret = un; \ 876 } \ 877 else \ 878 { \ 879 int __is_zero_x; \ 880 int __is_zero_y; \ 881 \ 882 __is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0; \ 883 __is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0; \ 884 \ 885 if (__is_zero_x && __is_zero_y) \ 886 ret = 0; \ 887 else if (__is_zero_x) \ 888 ret = Y##_s ? 1 : -1; \ 889 else if (__is_zero_y) \ 890 ret = X##_s ? -1 : 1; \ 891 else if (X##_s != Y##_s) \ 892 ret = X##_s ? -1 : 1; \ 893 else if (X##_e > Y##_e) \ 894 ret = X##_s ? -1 : 1; \ 895 else if (X##_e < Y##_e) \ 896 ret = X##_s ? 1 : -1; \ 897 else if (_FP_FRAC_GT_##wc(X, Y)) \ 898 ret = X##_s ? -1 : 1; \ 899 else if (_FP_FRAC_GT_##wc(Y, X)) \ 900 ret = X##_s ? 1 : -1; \ 901 else \ 902 ret = 0; \ 903 } \ 904 } while (0) 905 906 907 /* Simplification for strict equality. */ 908 909 #define _FP_CMP_EQ(fs, wc, ret, X, Y) \ 910 do { \ 911 /* NANs are unordered */ \ 912 if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \ 913 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))) \ 914 { \ 915 ret = 1; \ 916 } \ 917 else \ 918 { \ 919 ret = !(X##_e == Y##_e \ 920 && _FP_FRAC_EQ_##wc(X, Y) \ 921 && (X##_s == Y##_s || (!X##_e && _FP_FRAC_ZEROP_##wc(X)))); \ 922 } \ 923 } while (0) 924 925 /* Version to test unordered. */ 926 927 #define _FP_CMP_UNORD(fs, wc, ret, X, Y) \ 928 do { \ 929 ret = ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X)) \ 930 || (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y))); \ 931 } while (0) 932 933 /* 934 * Main square root routine. The input value should be cooked. 935 */ 936 937 #define _FP_SQRT(fs, wc, R, X) \ 938 do { \ 939 _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S); \ 940 _FP_W_TYPE q; \ 941 switch (X##_c) \ 942 { \ 943 case FP_CLS_NAN: \ 944 _FP_FRAC_COPY_##wc(R, X); \ 945 R##_s = X##_s; \ 946 R##_c = FP_CLS_NAN; \ 947 break; \ 948 case FP_CLS_INF: \ 949 if (X##_s) \ 950 { \ 951 R##_s = _FP_NANSIGN_##fs; \ 952 R##_c = FP_CLS_NAN; /* NAN */ \ 953 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \ 954 FP_SET_EXCEPTION(FP_EX_INVALID); \ 955 } \ 956 else \ 957 { \ 958 R##_s = 0; \ 959 R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */ \ 960 } \ 961 break; \ 962 case FP_CLS_ZERO: \ 963 R##_s = X##_s; \ 964 R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */ \ 965 break; \ 966 case FP_CLS_NORMAL: \ 967 R##_s = 0; \ 968 if (X##_s) \ 969 { \ 970 R##_c = FP_CLS_NAN; /* sNAN */ \ 971 R##_s = _FP_NANSIGN_##fs; \ 972 _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs); \ 973 FP_SET_EXCEPTION(FP_EX_INVALID); \ 974 break; \ 975 } \ 976 R##_c = FP_CLS_NORMAL; \ 977 if (X##_e & 1) \ 978 _FP_FRAC_SLL_##wc(X, 1); \ 979 R##_e = X##_e >> 1; \ 980 _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc); \ 981 _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc); \ 982 q = _FP_OVERFLOW_##fs >> 1; \ 983 _FP_SQRT_MEAT_##wc(R, S, T, X, q); \ 984 } \ 985 } while (0) 986 987 /* 988 * Convert from FP to integer. Input is raw. 989 */ 990 991 /* RSIGNED can have following values: 992 * 0: the number is required to be 0..(2^rsize)-1, if not, NV is set plus 993 * the result is either 0 or (2^rsize)-1 depending on the sign in such 994 * case. 995 * 1: the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not, 996 * NV is set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1 997 * depending on the sign in such case. 998 * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is 999 * set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1 1000 * depending on the sign in such case. 1001 */ 1002 #define _FP_TO_INT(fs, wc, r, X, rsize, rsigned) \ 1003 do { \ 1004 if (X##_e < _FP_EXPBIAS_##fs) \ 1005 { \ 1006 r = 0; \ 1007 if (X##_e == 0) \ 1008 { \ 1009 if (!_FP_FRAC_ZEROP_##wc(X)) \ 1010 { \ 1011 FP_SET_EXCEPTION(FP_EX_INEXACT); \ 1012 FP_SET_EXCEPTION(FP_EX_DENORM); \ 1013 } \ 1014 } \ 1015 else \ 1016 FP_SET_EXCEPTION(FP_EX_INEXACT); \ 1017 } \ 1018 else if (X##_e >= _FP_EXPBIAS_##fs + rsize - (rsigned > 0 || X##_s) \ 1019 || (!rsigned && X##_s)) \ 1020 { \ 1021 /* Overflow or converting to the most negative integer. */ \ 1022 if (rsigned) \ 1023 { \ 1024 r = 1; \ 1025 r <<= rsize - 1; \ 1026 r -= 1 - X##_s; \ 1027 } else { \ 1028 r = 0; \ 1029 if (X##_s) \ 1030 r = ~r; \ 1031 } \ 1032 \ 1033 if (rsigned && X##_s && X##_e == _FP_EXPBIAS_##fs + rsize - 1) \ 1034 { \ 1035 /* Possibly converting to most negative integer; check the \ 1036 mantissa. */ \ 1037 int inexact = 0; \ 1038 (void)((_FP_FRACBITS_##fs > rsize) \ 1039 ? ({ _FP_FRAC_SRST_##wc(X, inexact, \ 1040 _FP_FRACBITS_##fs - rsize, \ 1041 _FP_FRACBITS_##fs); 0; }) \ 1042 : 0); \ 1043 if (!_FP_FRAC_ZEROP_##wc(X)) \ 1044 FP_SET_EXCEPTION(FP_EX_INVALID); \ 1045 else if (inexact) \ 1046 FP_SET_EXCEPTION(FP_EX_INEXACT); \ 1047 } \ 1048 else \ 1049 FP_SET_EXCEPTION(FP_EX_INVALID); \ 1050 } \ 1051 else \ 1052 { \ 1053 _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs; \ 1054 if (X##_e >= _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs - 1) \ 1055 { \ 1056 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \ 1057 r <<= X##_e - _FP_EXPBIAS_##fs - _FP_FRACBITS_##fs + 1; \ 1058 } \ 1059 else \ 1060 { \ 1061 int inexact; \ 1062 _FP_FRAC_SRST_##wc(X, inexact, \ 1063 (_FP_FRACBITS_##fs + _FP_EXPBIAS_##fs - 1 \ 1064 - X##_e), \ 1065 _FP_FRACBITS_##fs); \ 1066 if (inexact) \ 1067 FP_SET_EXCEPTION(FP_EX_INEXACT); \ 1068 _FP_FRAC_ASSEMBLE_##wc(r, X, rsize); \ 1069 } \ 1070 if (rsigned && X##_s) \ 1071 r = -r; \ 1072 } \ 1073 } while (0) 1074 1075 /* Convert integer to fp. Output is raw. RTYPE is unsigned even if 1076 input is signed. */ 1077 #define _FP_FROM_INT(fs, wc, X, r, rsize, rtype) \ 1078 do { \ 1079 if (r) \ 1080 { \ 1081 rtype ur_; \ 1082 \ 1083 if ((X##_s = (r < 0))) \ 1084 r = -(rtype)r; \ 1085 \ 1086 ur_ = (rtype) r; \ 1087 (void)((rsize <= _FP_W_TYPE_SIZE) \ 1088 ? ({ \ 1089 int lz_; \ 1090 __FP_CLZ(lz_, (_FP_W_TYPE)ur_); \ 1091 X##_e = _FP_EXPBIAS_##fs + _FP_W_TYPE_SIZE - 1 - lz_; \ 1092 }) \ 1093 : ((rsize <= 2 * _FP_W_TYPE_SIZE) \ 1094 ? ({ \ 1095 int lz_; \ 1096 __FP_CLZ_2(lz_, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE), \ 1097 (_FP_W_TYPE)ur_); \ 1098 X##_e = (_FP_EXPBIAS_##fs + 2 * _FP_W_TYPE_SIZE - 1 \ 1099 - lz_); \ 1100 }) \ 1101 : (abort(), 0))); \ 1102 \ 1103 if (rsize - 1 + _FP_EXPBIAS_##fs >= _FP_EXPMAX_##fs \ 1104 && X##_e >= _FP_EXPMAX_##fs) \ 1105 { \ 1106 /* Exponent too big; overflow to infinity. (May also \ 1107 happen after rounding below.) */ \ 1108 _FP_OVERFLOW_SEMIRAW(fs, wc, X); \ 1109 goto pack_semiraw; \ 1110 } \ 1111 \ 1112 if (rsize <= _FP_FRACBITS_##fs \ 1113 || X##_e < _FP_EXPBIAS_##fs + _FP_FRACBITS_##fs) \ 1114 { \ 1115 /* Exactly representable; shift left. */ \ 1116 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \ 1117 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \ 1118 + _FP_FRACBITS_##fs - 1 - X##_e)); \ 1119 } \ 1120 else \ 1121 { \ 1122 /* More bits in integer than in floating type; need to \ 1123 round. */ \ 1124 if (_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 < X##_e) \ 1125 ur_ = ((ur_ >> (X##_e - _FP_EXPBIAS_##fs \ 1126 - _FP_WFRACBITS_##fs + 1)) \ 1127 | ((ur_ << (rsize - (X##_e - _FP_EXPBIAS_##fs \ 1128 - _FP_WFRACBITS_##fs + 1))) \ 1129 != 0)); \ 1130 _FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize); \ 1131 if ((_FP_EXPBIAS_##fs + _FP_WFRACBITS_##fs - 1 - X##_e) > 0) \ 1132 _FP_FRAC_SLL_##wc(X, (_FP_EXPBIAS_##fs \ 1133 + _FP_WFRACBITS_##fs - 1 - X##_e)); \ 1134 _FP_FRAC_HIGH_##fs(X) &= ~(_FP_W_TYPE)_FP_IMPLBIT_SH_##fs; \ 1135 pack_semiraw: \ 1136 _FP_PACK_SEMIRAW(fs, wc, X); \ 1137 } \ 1138 } \ 1139 else \ 1140 { \ 1141 X##_s = 0; \ 1142 X##_e = 0; \ 1143 _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc); \ 1144 } \ 1145 } while (0) 1146 1147 1148 /* Extend from a narrower floating-point format to a wider one. Input 1149 and output are raw. */ 1150 #define FP_EXTEND(dfs,sfs,dwc,swc,D,S) \ 1151 do { \ 1152 if (_FP_FRACBITS_##dfs < _FP_FRACBITS_##sfs \ 1153 || (_FP_EXPMAX_##dfs - _FP_EXPBIAS_##dfs \ 1154 < _FP_EXPMAX_##sfs - _FP_EXPBIAS_##sfs) \ 1155 || (_FP_EXPBIAS_##dfs < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1 \ 1156 && _FP_EXPBIAS_##dfs != _FP_EXPBIAS_##sfs)) \ 1157 abort(); \ 1158 D##_s = S##_s; \ 1159 _FP_FRAC_COPY_##dwc##_##swc(D, S); \ 1160 if (_FP_EXP_NORMAL(sfs, swc, S)) \ 1161 { \ 1162 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \ 1163 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs - _FP_FRACBITS_##sfs)); \ 1164 } \ 1165 else \ 1166 { \ 1167 if (S##_e == 0) \ 1168 { \ 1169 if (_FP_FRAC_ZEROP_##swc(S)) \ 1170 D##_e = 0; \ 1171 else if (_FP_EXPBIAS_##dfs \ 1172 < _FP_EXPBIAS_##sfs + _FP_FRACBITS_##sfs - 1) \ 1173 { \ 1174 FP_SET_EXCEPTION(FP_EX_DENORM); \ 1175 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \ 1176 - _FP_FRACBITS_##sfs)); \ 1177 D##_e = 0; \ 1178 } \ 1179 else \ 1180 { \ 1181 int _lz; \ 1182 FP_SET_EXCEPTION(FP_EX_DENORM); \ 1183 _FP_FRAC_CLZ_##swc(_lz, S); \ 1184 _FP_FRAC_SLL_##dwc(D, \ 1185 _lz + _FP_FRACBITS_##dfs \ 1186 - _FP_FRACTBITS_##sfs); \ 1187 D##_e = (_FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs + 1 \ 1188 + _FP_FRACXBITS_##sfs - _lz); \ 1189 } \ 1190 } \ 1191 else \ 1192 { \ 1193 D##_e = _FP_EXPMAX_##dfs; \ 1194 if (!_FP_FRAC_ZEROP_##swc(S)) \ 1195 { \ 1196 if (!(_FP_FRAC_HIGH_RAW_##sfs(S) & _FP_QNANBIT_##sfs)) \ 1197 FP_SET_EXCEPTION(FP_EX_INVALID); \ 1198 _FP_FRAC_SLL_##dwc(D, (_FP_FRACBITS_##dfs \ 1199 - _FP_FRACBITS_##sfs)); \ 1200 } \ 1201 } \ 1202 } \ 1203 } while (0) 1204 1205 /* Truncate from a wider floating-point format to a narrower one. 1206 Input and output are semi-raw. */ 1207 #define FP_TRUNC(dfs,sfs,dwc,swc,D,S) \ 1208 do { \ 1209 if (_FP_FRACBITS_##sfs < _FP_FRACBITS_##dfs \ 1210 || (_FP_EXPBIAS_##sfs < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1 \ 1211 && _FP_EXPBIAS_##sfs != _FP_EXPBIAS_##dfs)) \ 1212 abort(); \ 1213 D##_s = S##_s; \ 1214 if (_FP_EXP_NORMAL(sfs, swc, S)) \ 1215 { \ 1216 D##_e = S##_e + _FP_EXPBIAS_##dfs - _FP_EXPBIAS_##sfs; \ 1217 if (D##_e >= _FP_EXPMAX_##dfs) \ 1218 _FP_OVERFLOW_SEMIRAW(dfs, dwc, D); \ 1219 else \ 1220 { \ 1221 if (D##_e <= 0) \ 1222 { \ 1223 if (D##_e < 1 - _FP_FRACBITS_##dfs) \ 1224 { \ 1225 _FP_FRAC_SET_##swc(S, _FP_ZEROFRAC_##swc); \ 1226 _FP_FRAC_LOW_##swc(S) |= 1; \ 1227 } \ 1228 else \ 1229 { \ 1230 _FP_FRAC_HIGH_##sfs(S) |= _FP_IMPLBIT_SH_##sfs; \ 1231 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \ 1232 - _FP_WFRACBITS_##dfs + 1 - D##_e), \ 1233 _FP_WFRACBITS_##sfs); \ 1234 } \ 1235 D##_e = 0; \ 1236 } \ 1237 else \ 1238 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \ 1239 - _FP_WFRACBITS_##dfs), \ 1240 _FP_WFRACBITS_##sfs); \ 1241 _FP_FRAC_COPY_##dwc##_##swc(D, S); \ 1242 } \ 1243 } \ 1244 else \ 1245 { \ 1246 if (S##_e == 0) \ 1247 { \ 1248 D##_e = 0; \ 1249 if (_FP_FRAC_ZEROP_##swc(S)) \ 1250 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \ 1251 else \ 1252 { \ 1253 FP_SET_EXCEPTION(FP_EX_DENORM); \ 1254 if (_FP_EXPBIAS_##sfs \ 1255 < _FP_EXPBIAS_##dfs + _FP_FRACBITS_##dfs - 1) \ 1256 { \ 1257 _FP_FRAC_SRS_##swc(S, (_FP_WFRACBITS_##sfs \ 1258 - _FP_WFRACBITS_##dfs), \ 1259 _FP_WFRACBITS_##sfs); \ 1260 _FP_FRAC_COPY_##dwc##_##swc(D, S); \ 1261 } \ 1262 else \ 1263 { \ 1264 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \ 1265 _FP_FRAC_LOW_##dwc(D) |= 1; \ 1266 } \ 1267 } \ 1268 } \ 1269 else \ 1270 { \ 1271 D##_e = _FP_EXPMAX_##dfs; \ 1272 if (_FP_FRAC_ZEROP_##swc(S)) \ 1273 _FP_FRAC_SET_##dwc(D, _FP_ZEROFRAC_##dwc); \ 1274 else \ 1275 { \ 1276 _FP_CHECK_SIGNAN_SEMIRAW(sfs, swc, S); \ 1277 _FP_FRAC_SRL_##swc(S, (_FP_WFRACBITS_##sfs \ 1278 - _FP_WFRACBITS_##dfs)); \ 1279 _FP_FRAC_COPY_##dwc##_##swc(D, S); \ 1280 /* Semi-raw NaN must have all workbits cleared. */ \ 1281 _FP_FRAC_LOW_##dwc(D) \ 1282 &= ~(_FP_W_TYPE) ((1 << _FP_WORKBITS) - 1); \ 1283 _FP_FRAC_HIGH_##dfs(D) |= _FP_QNANBIT_SH_##dfs; \ 1284 } \ 1285 } \ 1286 } \ 1287 } while (0) 1288 1289 /* 1290 * Helper primitives. 1291 */ 1292 1293 /* Count leading zeros in a word. */ 1294 1295 #ifndef __FP_CLZ 1296 /* GCC 3.4 and later provide the builtins for us. */ 1297 #define __FP_CLZ(r, x) \ 1298 do { \ 1299 if (sizeof (_FP_W_TYPE) == sizeof (unsigned int)) \ 1300 r = __builtin_clz (x); \ 1301 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long)) \ 1302 r = __builtin_clzl (x); \ 1303 else if (sizeof (_FP_W_TYPE) == sizeof (unsigned long long)) \ 1304 r = __builtin_clzll (x); \ 1305 else \ 1306 abort (); \ 1307 } while (0) 1308 #endif /* ndef __FP_CLZ */ 1309 1310 #define _FP_DIV_HELP_imm(q, r, n, d) \ 1311 do { \ 1312 q = n / d, r = n % d; \ 1313 } while (0) 1314 1315 1316 /* A restoring bit-by-bit division primitive. */ 1317 1318 #define _FP_DIV_MEAT_N_loop(fs, wc, R, X, Y) \ 1319 do { \ 1320 int count = _FP_WFRACBITS_##fs; \ 1321 _FP_FRAC_DECL_##wc (u); \ 1322 _FP_FRAC_DECL_##wc (v); \ 1323 _FP_FRAC_COPY_##wc (u, X); \ 1324 _FP_FRAC_COPY_##wc (v, Y); \ 1325 _FP_FRAC_SET_##wc (R, _FP_ZEROFRAC_##wc); \ 1326 /* Normalize U and V. */ \ 1327 _FP_FRAC_SLL_##wc (u, _FP_WFRACXBITS_##fs); \ 1328 _FP_FRAC_SLL_##wc (v, _FP_WFRACXBITS_##fs); \ 1329 /* First round. Since the operands are normalized, either the \ 1330 first or second bit will be set in the fraction. Produce a \ 1331 normalized result by checking which and adjusting the loop \ 1332 count and exponent accordingly. */ \ 1333 if (_FP_FRAC_GE_1 (u, v)) \ 1334 { \ 1335 _FP_FRAC_SUB_##wc (u, u, v); \ 1336 _FP_FRAC_LOW_##wc (R) |= 1; \ 1337 count--; \ 1338 } \ 1339 else \ 1340 R##_e--; \ 1341 /* Subsequent rounds. */ \ 1342 do { \ 1343 int msb = (_FP_WS_TYPE) _FP_FRAC_HIGH_##wc (u) < 0; \ 1344 _FP_FRAC_SLL_##wc (u, 1); \ 1345 _FP_FRAC_SLL_##wc (R, 1); \ 1346 if (msb || _FP_FRAC_GE_1 (u, v)) \ 1347 { \ 1348 _FP_FRAC_SUB_##wc (u, u, v); \ 1349 _FP_FRAC_LOW_##wc (R) |= 1; \ 1350 } \ 1351 } while (--count > 0); \ 1352 /* If there's anything left in U, the result is inexact. */ \ 1353 _FP_FRAC_LOW_##wc (R) |= !_FP_FRAC_ZEROP_##wc (u); \ 1354 } while (0) 1355 1356 #define _FP_DIV_MEAT_1_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 1, R, X, Y) 1357 #define _FP_DIV_MEAT_2_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 2, R, X, Y) 1358 #define _FP_DIV_MEAT_4_loop(fs, R, X, Y) _FP_DIV_MEAT_N_loop (fs, 4, R, X, Y) 1359