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