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