1 /* 2 =============================================================================== 3 4 This C source file is part of the SoftFloat IEC/IEEE Floating-point 5 Arithmetic Package, Release 2. 6 7 Written by John R. Hauser. This work was made possible in part by the 8 International Computer Science Institute, located at Suite 600, 1947 Center 9 Street, Berkeley, California 94704. Funding was partially provided by the 10 National Science Foundation under grant MIP-9311980. The original version 11 of this code was written as part of a project to build a fixed-point vector 12 processor in collaboration with the University of California at Berkeley, 13 overseen by Profs. Nelson Morgan and John Wawrzynek. More information 14 is available through the web page 15 http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt 16 17 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort 18 has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT 19 TIMES RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO 20 PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY 21 AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE. 22 23 Derivative works are acceptable, even for commercial purposes, so long as 24 (1) they include prominent notice that the work is derivative, and (2) they 25 include prominent notice akin to these three paragraphs for those parts of 26 this code that are retained. 27 28 =============================================================================== 29 */ 30 31 #include <asm/div64.h> 32 33 #include "fpa11.h" 34 //#include "milieu.h" 35 //#include "softfloat.h" 36 37 /* 38 ------------------------------------------------------------------------------- 39 Primitive arithmetic functions, including multi-word arithmetic, and 40 division and square root approximations. (Can be specialized to target if 41 desired.) 42 ------------------------------------------------------------------------------- 43 */ 44 #include "softfloat-macros" 45 46 /* 47 ------------------------------------------------------------------------------- 48 Functions and definitions to determine: (1) whether tininess for underflow 49 is detected before or after rounding by default, (2) what (if anything) 50 happens when exceptions are raised, (3) how signaling NaNs are distinguished 51 from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs 52 are propagated from function inputs to output. These details are target- 53 specific. 54 ------------------------------------------------------------------------------- 55 */ 56 #include "softfloat-specialize" 57 58 /* 59 ------------------------------------------------------------------------------- 60 Takes a 64-bit fixed-point value `absZ' with binary point between bits 6 61 and 7, and returns the properly rounded 32-bit integer corresponding to the 62 input. If `zSign' is nonzero, the input is negated before being converted 63 to an integer. Bit 63 of `absZ' must be zero. Ordinarily, the fixed-point 64 input is simply rounded to an integer, with the inexact exception raised if 65 the input cannot be represented exactly as an integer. If the fixed-point 66 input is too large, however, the invalid exception is raised and the largest 67 positive or negative integer is returned. 68 ------------------------------------------------------------------------------- 69 */ 70 static int32 roundAndPackInt32( struct roundingData *roundData, flag zSign, bits64 absZ ) 71 { 72 int8 roundingMode; 73 flag roundNearestEven; 74 int8 roundIncrement, roundBits; 75 int32 z; 76 77 roundingMode = roundData->mode; 78 roundNearestEven = ( roundingMode == float_round_nearest_even ); 79 roundIncrement = 0x40; 80 if ( ! roundNearestEven ) { 81 if ( roundingMode == float_round_to_zero ) { 82 roundIncrement = 0; 83 } 84 else { 85 roundIncrement = 0x7F; 86 if ( zSign ) { 87 if ( roundingMode == float_round_up ) roundIncrement = 0; 88 } 89 else { 90 if ( roundingMode == float_round_down ) roundIncrement = 0; 91 } 92 } 93 } 94 roundBits = absZ & 0x7F; 95 absZ = ( absZ + roundIncrement )>>7; 96 absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven ); 97 z = absZ; 98 if ( zSign ) z = - z; 99 if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) { 100 roundData->exception |= float_flag_invalid; 101 return zSign ? 0x80000000 : 0x7FFFFFFF; 102 } 103 if ( roundBits ) roundData->exception |= float_flag_inexact; 104 return z; 105 106 } 107 108 /* 109 ------------------------------------------------------------------------------- 110 Returns the fraction bits of the single-precision floating-point value `a'. 111 ------------------------------------------------------------------------------- 112 */ 113 INLINE bits32 extractFloat32Frac( float32 a ) 114 { 115 116 return a & 0x007FFFFF; 117 118 } 119 120 /* 121 ------------------------------------------------------------------------------- 122 Returns the exponent bits of the single-precision floating-point value `a'. 123 ------------------------------------------------------------------------------- 124 */ 125 INLINE int16 extractFloat32Exp( float32 a ) 126 { 127 128 return ( a>>23 ) & 0xFF; 129 130 } 131 132 /* 133 ------------------------------------------------------------------------------- 134 Returns the sign bit of the single-precision floating-point value `a'. 135 ------------------------------------------------------------------------------- 136 */ 137 #if 0 /* in softfloat.h */ 138 INLINE flag extractFloat32Sign( float32 a ) 139 { 140 141 return a>>31; 142 143 } 144 #endif 145 146 /* 147 ------------------------------------------------------------------------------- 148 Normalizes the subnormal single-precision floating-point value represented 149 by the denormalized significand `aSig'. The normalized exponent and 150 significand are stored at the locations pointed to by `zExpPtr' and 151 `zSigPtr', respectively. 152 ------------------------------------------------------------------------------- 153 */ 154 static void 155 normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr ) 156 { 157 int8 shiftCount; 158 159 shiftCount = countLeadingZeros32( aSig ) - 8; 160 *zSigPtr = aSig<<shiftCount; 161 *zExpPtr = 1 - shiftCount; 162 163 } 164 165 /* 166 ------------------------------------------------------------------------------- 167 Packs the sign `zSign', exponent `zExp', and significand `zSig' into a 168 single-precision floating-point value, returning the result. After being 169 shifted into the proper positions, the three fields are simply added 170 together to form the result. This means that any integer portion of `zSig' 171 will be added into the exponent. Since a properly normalized significand 172 will have an integer portion equal to 1, the `zExp' input should be 1 less 173 than the desired result exponent whenever `zSig' is a complete, normalized 174 significand. 175 ------------------------------------------------------------------------------- 176 */ 177 INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig ) 178 { 179 #if 0 180 float32 f; 181 __asm__("@ packFloat32 \n\ 182 mov %0, %1, asl #31 \n\ 183 orr %0, %2, asl #23 \n\ 184 orr %0, %3" 185 : /* no outputs */ 186 : "g" (f), "g" (zSign), "g" (zExp), "g" (zSig) 187 : "cc"); 188 return f; 189 #else 190 return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig; 191 #endif 192 } 193 194 /* 195 ------------------------------------------------------------------------------- 196 Takes an abstract floating-point value having sign `zSign', exponent `zExp', 197 and significand `zSig', and returns the proper single-precision floating- 198 point value corresponding to the abstract input. Ordinarily, the abstract 199 value is simply rounded and packed into the single-precision format, with 200 the inexact exception raised if the abstract input cannot be represented 201 exactly. If the abstract value is too large, however, the overflow and 202 inexact exceptions are raised and an infinity or maximal finite value is 203 returned. If the abstract value is too small, the input value is rounded to 204 a subnormal number, and the underflow and inexact exceptions are raised if 205 the abstract input cannot be represented exactly as a subnormal single- 206 precision floating-point number. 207 The input significand `zSig' has its binary point between bits 30 208 and 29, which is 7 bits to the left of the usual location. This shifted 209 significand must be normalized or smaller. If `zSig' is not normalized, 210 `zExp' must be 0; in that case, the result returned is a subnormal number, 211 and it must not require rounding. In the usual case that `zSig' is 212 normalized, `zExp' must be 1 less than the ``true'' floating-point exponent. 213 The handling of underflow and overflow follows the IEC/IEEE Standard for 214 Binary Floating-point Arithmetic. 215 ------------------------------------------------------------------------------- 216 */ 217 static float32 roundAndPackFloat32( struct roundingData *roundData, flag zSign, int16 zExp, bits32 zSig ) 218 { 219 int8 roundingMode; 220 flag roundNearestEven; 221 int8 roundIncrement, roundBits; 222 flag isTiny; 223 224 roundingMode = roundData->mode; 225 roundNearestEven = ( roundingMode == float_round_nearest_even ); 226 roundIncrement = 0x40; 227 if ( ! roundNearestEven ) { 228 if ( roundingMode == float_round_to_zero ) { 229 roundIncrement = 0; 230 } 231 else { 232 roundIncrement = 0x7F; 233 if ( zSign ) { 234 if ( roundingMode == float_round_up ) roundIncrement = 0; 235 } 236 else { 237 if ( roundingMode == float_round_down ) roundIncrement = 0; 238 } 239 } 240 } 241 roundBits = zSig & 0x7F; 242 if ( 0xFD <= (bits16) zExp ) { 243 if ( ( 0xFD < zExp ) 244 || ( ( zExp == 0xFD ) 245 && ( (sbits32) ( zSig + roundIncrement ) < 0 ) ) 246 ) { 247 roundData->exception |= float_flag_overflow | float_flag_inexact; 248 return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 ); 249 } 250 if ( zExp < 0 ) { 251 isTiny = 252 ( float_detect_tininess == float_tininess_before_rounding ) 253 || ( zExp < -1 ) 254 || ( zSig + roundIncrement < 0x80000000 ); 255 shift32RightJamming( zSig, - zExp, &zSig ); 256 zExp = 0; 257 roundBits = zSig & 0x7F; 258 if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow; 259 } 260 } 261 if ( roundBits ) roundData->exception |= float_flag_inexact; 262 zSig = ( zSig + roundIncrement )>>7; 263 zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven ); 264 if ( zSig == 0 ) zExp = 0; 265 return packFloat32( zSign, zExp, zSig ); 266 267 } 268 269 /* 270 ------------------------------------------------------------------------------- 271 Takes an abstract floating-point value having sign `zSign', exponent `zExp', 272 and significand `zSig', and returns the proper single-precision floating- 273 point value corresponding to the abstract input. This routine is just like 274 `roundAndPackFloat32' except that `zSig' does not have to be normalized in 275 any way. In all cases, `zExp' must be 1 less than the ``true'' floating- 276 point exponent. 277 ------------------------------------------------------------------------------- 278 */ 279 static float32 280 normalizeRoundAndPackFloat32( struct roundingData *roundData, flag zSign, int16 zExp, bits32 zSig ) 281 { 282 int8 shiftCount; 283 284 shiftCount = countLeadingZeros32( zSig ) - 1; 285 return roundAndPackFloat32( roundData, zSign, zExp - shiftCount, zSig<<shiftCount ); 286 287 } 288 289 /* 290 ------------------------------------------------------------------------------- 291 Returns the fraction bits of the double-precision floating-point value `a'. 292 ------------------------------------------------------------------------------- 293 */ 294 INLINE bits64 extractFloat64Frac( float64 a ) 295 { 296 297 return a & LIT64( 0x000FFFFFFFFFFFFF ); 298 299 } 300 301 /* 302 ------------------------------------------------------------------------------- 303 Returns the exponent bits of the double-precision floating-point value `a'. 304 ------------------------------------------------------------------------------- 305 */ 306 INLINE int16 extractFloat64Exp( float64 a ) 307 { 308 309 return ( a>>52 ) & 0x7FF; 310 311 } 312 313 /* 314 ------------------------------------------------------------------------------- 315 Returns the sign bit of the double-precision floating-point value `a'. 316 ------------------------------------------------------------------------------- 317 */ 318 #if 0 /* in softfloat.h */ 319 INLINE flag extractFloat64Sign( float64 a ) 320 { 321 322 return a>>63; 323 324 } 325 #endif 326 327 /* 328 ------------------------------------------------------------------------------- 329 Normalizes the subnormal double-precision floating-point value represented 330 by the denormalized significand `aSig'. The normalized exponent and 331 significand are stored at the locations pointed to by `zExpPtr' and 332 `zSigPtr', respectively. 333 ------------------------------------------------------------------------------- 334 */ 335 static void 336 normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr ) 337 { 338 int8 shiftCount; 339 340 shiftCount = countLeadingZeros64( aSig ) - 11; 341 *zSigPtr = aSig<<shiftCount; 342 *zExpPtr = 1 - shiftCount; 343 344 } 345 346 /* 347 ------------------------------------------------------------------------------- 348 Packs the sign `zSign', exponent `zExp', and significand `zSig' into a 349 double-precision floating-point value, returning the result. After being 350 shifted into the proper positions, the three fields are simply added 351 together to form the result. This means that any integer portion of `zSig' 352 will be added into the exponent. Since a properly normalized significand 353 will have an integer portion equal to 1, the `zExp' input should be 1 less 354 than the desired result exponent whenever `zSig' is a complete, normalized 355 significand. 356 ------------------------------------------------------------------------------- 357 */ 358 INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig ) 359 { 360 361 return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig; 362 363 } 364 365 /* 366 ------------------------------------------------------------------------------- 367 Takes an abstract floating-point value having sign `zSign', exponent `zExp', 368 and significand `zSig', and returns the proper double-precision floating- 369 point value corresponding to the abstract input. Ordinarily, the abstract 370 value is simply rounded and packed into the double-precision format, with 371 the inexact exception raised if the abstract input cannot be represented 372 exactly. If the abstract value is too large, however, the overflow and 373 inexact exceptions are raised and an infinity or maximal finite value is 374 returned. If the abstract value is too small, the input value is rounded to 375 a subnormal number, and the underflow and inexact exceptions are raised if 376 the abstract input cannot be represented exactly as a subnormal double- 377 precision floating-point number. 378 The input significand `zSig' has its binary point between bits 62 379 and 61, which is 10 bits to the left of the usual location. This shifted 380 significand must be normalized or smaller. If `zSig' is not normalized, 381 `zExp' must be 0; in that case, the result returned is a subnormal number, 382 and it must not require rounding. In the usual case that `zSig' is 383 normalized, `zExp' must be 1 less than the ``true'' floating-point exponent. 384 The handling of underflow and overflow follows the IEC/IEEE Standard for 385 Binary Floating-point Arithmetic. 386 ------------------------------------------------------------------------------- 387 */ 388 static float64 roundAndPackFloat64( struct roundingData *roundData, flag zSign, int16 zExp, bits64 zSig ) 389 { 390 int8 roundingMode; 391 flag roundNearestEven; 392 int16 roundIncrement, roundBits; 393 flag isTiny; 394 395 roundingMode = roundData->mode; 396 roundNearestEven = ( roundingMode == float_round_nearest_even ); 397 roundIncrement = 0x200; 398 if ( ! roundNearestEven ) { 399 if ( roundingMode == float_round_to_zero ) { 400 roundIncrement = 0; 401 } 402 else { 403 roundIncrement = 0x3FF; 404 if ( zSign ) { 405 if ( roundingMode == float_round_up ) roundIncrement = 0; 406 } 407 else { 408 if ( roundingMode == float_round_down ) roundIncrement = 0; 409 } 410 } 411 } 412 roundBits = zSig & 0x3FF; 413 if ( 0x7FD <= (bits16) zExp ) { 414 if ( ( 0x7FD < zExp ) 415 || ( ( zExp == 0x7FD ) 416 && ( (sbits64) ( zSig + roundIncrement ) < 0 ) ) 417 ) { 418 //register int lr = __builtin_return_address(0); 419 //printk("roundAndPackFloat64 called from 0x%08x\n",lr); 420 roundData->exception |= float_flag_overflow | float_flag_inexact; 421 return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 ); 422 } 423 if ( zExp < 0 ) { 424 isTiny = 425 ( float_detect_tininess == float_tininess_before_rounding ) 426 || ( zExp < -1 ) 427 || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) ); 428 shift64RightJamming( zSig, - zExp, &zSig ); 429 zExp = 0; 430 roundBits = zSig & 0x3FF; 431 if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow; 432 } 433 } 434 if ( roundBits ) roundData->exception |= float_flag_inexact; 435 zSig = ( zSig + roundIncrement )>>10; 436 zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven ); 437 if ( zSig == 0 ) zExp = 0; 438 return packFloat64( zSign, zExp, zSig ); 439 440 } 441 442 /* 443 ------------------------------------------------------------------------------- 444 Takes an abstract floating-point value having sign `zSign', exponent `zExp', 445 and significand `zSig', and returns the proper double-precision floating- 446 point value corresponding to the abstract input. This routine is just like 447 `roundAndPackFloat64' except that `zSig' does not have to be normalized in 448 any way. In all cases, `zExp' must be 1 less than the ``true'' floating- 449 point exponent. 450 ------------------------------------------------------------------------------- 451 */ 452 static float64 453 normalizeRoundAndPackFloat64( struct roundingData *roundData, flag zSign, int16 zExp, bits64 zSig ) 454 { 455 int8 shiftCount; 456 457 shiftCount = countLeadingZeros64( zSig ) - 1; 458 return roundAndPackFloat64( roundData, zSign, zExp - shiftCount, zSig<<shiftCount ); 459 460 } 461 462 #ifdef FLOATX80 463 464 /* 465 ------------------------------------------------------------------------------- 466 Returns the fraction bits of the extended double-precision floating-point 467 value `a'. 468 ------------------------------------------------------------------------------- 469 */ 470 INLINE bits64 extractFloatx80Frac( floatx80 a ) 471 { 472 473 return a.low; 474 475 } 476 477 /* 478 ------------------------------------------------------------------------------- 479 Returns the exponent bits of the extended double-precision floating-point 480 value `a'. 481 ------------------------------------------------------------------------------- 482 */ 483 INLINE int32 extractFloatx80Exp( floatx80 a ) 484 { 485 486 return a.high & 0x7FFF; 487 488 } 489 490 /* 491 ------------------------------------------------------------------------------- 492 Returns the sign bit of the extended double-precision floating-point value 493 `a'. 494 ------------------------------------------------------------------------------- 495 */ 496 INLINE flag extractFloatx80Sign( floatx80 a ) 497 { 498 499 return a.high>>15; 500 501 } 502 503 /* 504 ------------------------------------------------------------------------------- 505 Normalizes the subnormal extended double-precision floating-point value 506 represented by the denormalized significand `aSig'. The normalized exponent 507 and significand are stored at the locations pointed to by `zExpPtr' and 508 `zSigPtr', respectively. 509 ------------------------------------------------------------------------------- 510 */ 511 static void 512 normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr ) 513 { 514 int8 shiftCount; 515 516 shiftCount = countLeadingZeros64( aSig ); 517 *zSigPtr = aSig<<shiftCount; 518 *zExpPtr = 1 - shiftCount; 519 520 } 521 522 /* 523 ------------------------------------------------------------------------------- 524 Packs the sign `zSign', exponent `zExp', and significand `zSig' into an 525 extended double-precision floating-point value, returning the result. 526 ------------------------------------------------------------------------------- 527 */ 528 INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig ) 529 { 530 floatx80 z; 531 532 z.low = zSig; 533 z.high = ( ( (bits16) zSign )<<15 ) + zExp; 534 z.__padding = 0; 535 return z; 536 537 } 538 539 /* 540 ------------------------------------------------------------------------------- 541 Takes an abstract floating-point value having sign `zSign', exponent `zExp', 542 and extended significand formed by the concatenation of `zSig0' and `zSig1', 543 and returns the proper extended double-precision floating-point value 544 corresponding to the abstract input. Ordinarily, the abstract value is 545 rounded and packed into the extended double-precision format, with the 546 inexact exception raised if the abstract input cannot be represented 547 exactly. If the abstract value is too large, however, the overflow and 548 inexact exceptions are raised and an infinity or maximal finite value is 549 returned. If the abstract value is too small, the input value is rounded to 550 a subnormal number, and the underflow and inexact exceptions are raised if 551 the abstract input cannot be represented exactly as a subnormal extended 552 double-precision floating-point number. 553 If `roundingPrecision' is 32 or 64, the result is rounded to the same 554 number of bits as single or double precision, respectively. Otherwise, the 555 result is rounded to the full precision of the extended double-precision 556 format. 557 The input significand must be normalized or smaller. If the input 558 significand is not normalized, `zExp' must be 0; in that case, the result 559 returned is a subnormal number, and it must not require rounding. The 560 handling of underflow and overflow follows the IEC/IEEE Standard for Binary 561 Floating-point Arithmetic. 562 ------------------------------------------------------------------------------- 563 */ 564 static floatx80 565 roundAndPackFloatx80( 566 struct roundingData *roundData, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 567 ) 568 { 569 int8 roundingMode, roundingPrecision; 570 flag roundNearestEven, increment, isTiny; 571 int64 roundIncrement, roundMask, roundBits; 572 573 roundingMode = roundData->mode; 574 roundingPrecision = roundData->precision; 575 roundNearestEven = ( roundingMode == float_round_nearest_even ); 576 if ( roundingPrecision == 80 ) goto precision80; 577 if ( roundingPrecision == 64 ) { 578 roundIncrement = LIT64( 0x0000000000000400 ); 579 roundMask = LIT64( 0x00000000000007FF ); 580 } 581 else if ( roundingPrecision == 32 ) { 582 roundIncrement = LIT64( 0x0000008000000000 ); 583 roundMask = LIT64( 0x000000FFFFFFFFFF ); 584 } 585 else { 586 goto precision80; 587 } 588 zSig0 |= ( zSig1 != 0 ); 589 if ( ! roundNearestEven ) { 590 if ( roundingMode == float_round_to_zero ) { 591 roundIncrement = 0; 592 } 593 else { 594 roundIncrement = roundMask; 595 if ( zSign ) { 596 if ( roundingMode == float_round_up ) roundIncrement = 0; 597 } 598 else { 599 if ( roundingMode == float_round_down ) roundIncrement = 0; 600 } 601 } 602 } 603 roundBits = zSig0 & roundMask; 604 if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) { 605 if ( ( 0x7FFE < zExp ) 606 || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) ) 607 ) { 608 goto overflow; 609 } 610 if ( zExp <= 0 ) { 611 isTiny = 612 ( float_detect_tininess == float_tininess_before_rounding ) 613 || ( zExp < 0 ) 614 || ( zSig0 <= zSig0 + roundIncrement ); 615 shift64RightJamming( zSig0, 1 - zExp, &zSig0 ); 616 zExp = 0; 617 roundBits = zSig0 & roundMask; 618 if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow; 619 if ( roundBits ) roundData->exception |= float_flag_inexact; 620 zSig0 += roundIncrement; 621 if ( (sbits64) zSig0 < 0 ) zExp = 1; 622 roundIncrement = roundMask + 1; 623 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) { 624 roundMask |= roundIncrement; 625 } 626 zSig0 &= ~ roundMask; 627 return packFloatx80( zSign, zExp, zSig0 ); 628 } 629 } 630 if ( roundBits ) roundData->exception |= float_flag_inexact; 631 zSig0 += roundIncrement; 632 if ( zSig0 < roundIncrement ) { 633 ++zExp; 634 zSig0 = LIT64( 0x8000000000000000 ); 635 } 636 roundIncrement = roundMask + 1; 637 if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) { 638 roundMask |= roundIncrement; 639 } 640 zSig0 &= ~ roundMask; 641 if ( zSig0 == 0 ) zExp = 0; 642 return packFloatx80( zSign, zExp, zSig0 ); 643 precision80: 644 increment = ( (sbits64) zSig1 < 0 ); 645 if ( ! roundNearestEven ) { 646 if ( roundingMode == float_round_to_zero ) { 647 increment = 0; 648 } 649 else { 650 if ( zSign ) { 651 increment = ( roundingMode == float_round_down ) && zSig1; 652 } 653 else { 654 increment = ( roundingMode == float_round_up ) && zSig1; 655 } 656 } 657 } 658 if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) { 659 if ( ( 0x7FFE < zExp ) 660 || ( ( zExp == 0x7FFE ) 661 && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) ) 662 && increment 663 ) 664 ) { 665 roundMask = 0; 666 overflow: 667 roundData->exception |= float_flag_overflow | float_flag_inexact; 668 if ( ( roundingMode == float_round_to_zero ) 669 || ( zSign && ( roundingMode == float_round_up ) ) 670 || ( ! zSign && ( roundingMode == float_round_down ) ) 671 ) { 672 return packFloatx80( zSign, 0x7FFE, ~ roundMask ); 673 } 674 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 675 } 676 if ( zExp <= 0 ) { 677 isTiny = 678 ( float_detect_tininess == float_tininess_before_rounding ) 679 || ( zExp < 0 ) 680 || ! increment 681 || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) ); 682 shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 ); 683 zExp = 0; 684 if ( isTiny && zSig1 ) roundData->exception |= float_flag_underflow; 685 if ( zSig1 ) roundData->exception |= float_flag_inexact; 686 if ( roundNearestEven ) { 687 increment = ( (sbits64) zSig1 < 0 ); 688 } 689 else { 690 if ( zSign ) { 691 increment = ( roundingMode == float_round_down ) && zSig1; 692 } 693 else { 694 increment = ( roundingMode == float_round_up ) && zSig1; 695 } 696 } 697 if ( increment ) { 698 ++zSig0; 699 zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven ); 700 if ( (sbits64) zSig0 < 0 ) zExp = 1; 701 } 702 return packFloatx80( zSign, zExp, zSig0 ); 703 } 704 } 705 if ( zSig1 ) roundData->exception |= float_flag_inexact; 706 if ( increment ) { 707 ++zSig0; 708 if ( zSig0 == 0 ) { 709 ++zExp; 710 zSig0 = LIT64( 0x8000000000000000 ); 711 } 712 else { 713 zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven ); 714 } 715 } 716 else { 717 if ( zSig0 == 0 ) zExp = 0; 718 } 719 720 return packFloatx80( zSign, zExp, zSig0 ); 721 } 722 723 /* 724 ------------------------------------------------------------------------------- 725 Takes an abstract floating-point value having sign `zSign', exponent 726 `zExp', and significand formed by the concatenation of `zSig0' and `zSig1', 727 and returns the proper extended double-precision floating-point value 728 corresponding to the abstract input. This routine is just like 729 `roundAndPackFloatx80' except that the input significand does not have to be 730 normalized. 731 ------------------------------------------------------------------------------- 732 */ 733 static floatx80 734 normalizeRoundAndPackFloatx80( 735 struct roundingData *roundData, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 736 ) 737 { 738 int8 shiftCount; 739 740 if ( zSig0 == 0 ) { 741 zSig0 = zSig1; 742 zSig1 = 0; 743 zExp -= 64; 744 } 745 shiftCount = countLeadingZeros64( zSig0 ); 746 shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 ); 747 zExp -= shiftCount; 748 return 749 roundAndPackFloatx80( roundData, zSign, zExp, zSig0, zSig1 ); 750 751 } 752 753 #endif 754 755 /* 756 ------------------------------------------------------------------------------- 757 Returns the result of converting the 32-bit two's complement integer `a' to 758 the single-precision floating-point format. The conversion is performed 759 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 760 ------------------------------------------------------------------------------- 761 */ 762 float32 int32_to_float32(struct roundingData *roundData, int32 a) 763 { 764 flag zSign; 765 766 if ( a == 0 ) return 0; 767 if ( a == 0x80000000 ) return packFloat32( 1, 0x9E, 0 ); 768 zSign = ( a < 0 ); 769 return normalizeRoundAndPackFloat32( roundData, zSign, 0x9C, zSign ? - a : a ); 770 771 } 772 773 /* 774 ------------------------------------------------------------------------------- 775 Returns the result of converting the 32-bit two's complement integer `a' to 776 the double-precision floating-point format. The conversion is performed 777 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 778 ------------------------------------------------------------------------------- 779 */ 780 float64 int32_to_float64( int32 a ) 781 { 782 flag aSign; 783 uint32 absA; 784 int8 shiftCount; 785 bits64 zSig; 786 787 if ( a == 0 ) return 0; 788 aSign = ( a < 0 ); 789 absA = aSign ? - a : a; 790 shiftCount = countLeadingZeros32( absA ) + 21; 791 zSig = absA; 792 return packFloat64( aSign, 0x432 - shiftCount, zSig<<shiftCount ); 793 794 } 795 796 #ifdef FLOATX80 797 798 /* 799 ------------------------------------------------------------------------------- 800 Returns the result of converting the 32-bit two's complement integer `a' 801 to the extended double-precision floating-point format. The conversion 802 is performed according to the IEC/IEEE Standard for Binary Floating-point 803 Arithmetic. 804 ------------------------------------------------------------------------------- 805 */ 806 floatx80 int32_to_floatx80( int32 a ) 807 { 808 flag zSign; 809 uint32 absA; 810 int8 shiftCount; 811 bits64 zSig; 812 813 if ( a == 0 ) return packFloatx80( 0, 0, 0 ); 814 zSign = ( a < 0 ); 815 absA = zSign ? - a : a; 816 shiftCount = countLeadingZeros32( absA ) + 32; 817 zSig = absA; 818 return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount ); 819 820 } 821 822 #endif 823 824 /* 825 ------------------------------------------------------------------------------- 826 Returns the result of converting the single-precision floating-point value 827 `a' to the 32-bit two's complement integer format. The conversion is 828 performed according to the IEC/IEEE Standard for Binary Floating-point 829 Arithmetic---which means in particular that the conversion is rounded 830 according to the current rounding mode. If `a' is a NaN, the largest 831 positive integer is returned. Otherwise, if the conversion overflows, the 832 largest integer with the same sign as `a' is returned. 833 ------------------------------------------------------------------------------- 834 */ 835 int32 float32_to_int32( struct roundingData *roundData, float32 a ) 836 { 837 flag aSign; 838 int16 aExp, shiftCount; 839 bits32 aSig; 840 bits64 zSig; 841 842 aSig = extractFloat32Frac( a ); 843 aExp = extractFloat32Exp( a ); 844 aSign = extractFloat32Sign( a ); 845 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0; 846 if ( aExp ) aSig |= 0x00800000; 847 shiftCount = 0xAF - aExp; 848 zSig = aSig; 849 zSig <<= 32; 850 if ( 0 < shiftCount ) shift64RightJamming( zSig, shiftCount, &zSig ); 851 return roundAndPackInt32( roundData, aSign, zSig ); 852 853 } 854 855 /* 856 ------------------------------------------------------------------------------- 857 Returns the result of converting the single-precision floating-point value 858 `a' to the 32-bit two's complement integer format. The conversion is 859 performed according to the IEC/IEEE Standard for Binary Floating-point 860 Arithmetic, except that the conversion is always rounded toward zero. If 861 `a' is a NaN, the largest positive integer is returned. Otherwise, if the 862 conversion overflows, the largest integer with the same sign as `a' is 863 returned. 864 ------------------------------------------------------------------------------- 865 */ 866 int32 float32_to_int32_round_to_zero( float32 a ) 867 { 868 flag aSign; 869 int16 aExp, shiftCount; 870 bits32 aSig; 871 int32 z; 872 873 aSig = extractFloat32Frac( a ); 874 aExp = extractFloat32Exp( a ); 875 aSign = extractFloat32Sign( a ); 876 shiftCount = aExp - 0x9E; 877 if ( 0 <= shiftCount ) { 878 if ( a == 0xCF000000 ) return 0x80000000; 879 float_raise( float_flag_invalid ); 880 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF; 881 return 0x80000000; 882 } 883 else if ( aExp <= 0x7E ) { 884 if ( aExp | aSig ) float_raise( float_flag_inexact ); 885 return 0; 886 } 887 aSig = ( aSig | 0x00800000 )<<8; 888 z = aSig>>( - shiftCount ); 889 if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) { 890 float_raise( float_flag_inexact ); 891 } 892 return aSign ? - z : z; 893 894 } 895 896 /* 897 ------------------------------------------------------------------------------- 898 Returns the result of converting the single-precision floating-point value 899 `a' to the double-precision floating-point format. The conversion is 900 performed according to the IEC/IEEE Standard for Binary Floating-point 901 Arithmetic. 902 ------------------------------------------------------------------------------- 903 */ 904 float64 float32_to_float64( float32 a ) 905 { 906 flag aSign; 907 int16 aExp; 908 bits32 aSig; 909 910 aSig = extractFloat32Frac( a ); 911 aExp = extractFloat32Exp( a ); 912 aSign = extractFloat32Sign( a ); 913 if ( aExp == 0xFF ) { 914 if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) ); 915 return packFloat64( aSign, 0x7FF, 0 ); 916 } 917 if ( aExp == 0 ) { 918 if ( aSig == 0 ) return packFloat64( aSign, 0, 0 ); 919 normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 920 --aExp; 921 } 922 return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 ); 923 924 } 925 926 #ifdef FLOATX80 927 928 /* 929 ------------------------------------------------------------------------------- 930 Returns the result of converting the single-precision floating-point value 931 `a' to the extended double-precision floating-point format. The conversion 932 is performed according to the IEC/IEEE Standard for Binary Floating-point 933 Arithmetic. 934 ------------------------------------------------------------------------------- 935 */ 936 floatx80 float32_to_floatx80( float32 a ) 937 { 938 flag aSign; 939 int16 aExp; 940 bits32 aSig; 941 942 aSig = extractFloat32Frac( a ); 943 aExp = extractFloat32Exp( a ); 944 aSign = extractFloat32Sign( a ); 945 if ( aExp == 0xFF ) { 946 if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a ) ); 947 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 948 } 949 if ( aExp == 0 ) { 950 if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 ); 951 normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 952 } 953 aSig |= 0x00800000; 954 return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 ); 955 956 } 957 958 #endif 959 960 /* 961 ------------------------------------------------------------------------------- 962 Rounds the single-precision floating-point value `a' to an integer, and 963 returns the result as a single-precision floating-point value. The 964 operation is performed according to the IEC/IEEE Standard for Binary 965 Floating-point Arithmetic. 966 ------------------------------------------------------------------------------- 967 */ 968 float32 float32_round_to_int( struct roundingData *roundData, float32 a ) 969 { 970 flag aSign; 971 int16 aExp; 972 bits32 lastBitMask, roundBitsMask; 973 int8 roundingMode; 974 float32 z; 975 976 aExp = extractFloat32Exp( a ); 977 if ( 0x96 <= aExp ) { 978 if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) { 979 return propagateFloat32NaN( a, a ); 980 } 981 return a; 982 } 983 roundingMode = roundData->mode; 984 if ( aExp <= 0x7E ) { 985 if ( (bits32) ( a<<1 ) == 0 ) return a; 986 roundData->exception |= float_flag_inexact; 987 aSign = extractFloat32Sign( a ); 988 switch ( roundingMode ) { 989 case float_round_nearest_even: 990 if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) { 991 return packFloat32( aSign, 0x7F, 0 ); 992 } 993 break; 994 case float_round_down: 995 return aSign ? 0xBF800000 : 0; 996 case float_round_up: 997 return aSign ? 0x80000000 : 0x3F800000; 998 } 999 return packFloat32( aSign, 0, 0 ); 1000 } 1001 lastBitMask = 1; 1002 lastBitMask <<= 0x96 - aExp; 1003 roundBitsMask = lastBitMask - 1; 1004 z = a; 1005 if ( roundingMode == float_round_nearest_even ) { 1006 z += lastBitMask>>1; 1007 if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask; 1008 } 1009 else if ( roundingMode != float_round_to_zero ) { 1010 if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) { 1011 z += roundBitsMask; 1012 } 1013 } 1014 z &= ~ roundBitsMask; 1015 if ( z != a ) roundData->exception |= float_flag_inexact; 1016 return z; 1017 1018 } 1019 1020 /* 1021 ------------------------------------------------------------------------------- 1022 Returns the result of adding the absolute values of the single-precision 1023 floating-point values `a' and `b'. If `zSign' is true, the sum is negated 1024 before being returned. `zSign' is ignored if the result is a NaN. The 1025 addition is performed according to the IEC/IEEE Standard for Binary 1026 Floating-point Arithmetic. 1027 ------------------------------------------------------------------------------- 1028 */ 1029 static float32 addFloat32Sigs( struct roundingData *roundData, float32 a, float32 b, flag zSign ) 1030 { 1031 int16 aExp, bExp, zExp; 1032 bits32 aSig, bSig, zSig; 1033 int16 expDiff; 1034 1035 aSig = extractFloat32Frac( a ); 1036 aExp = extractFloat32Exp( a ); 1037 bSig = extractFloat32Frac( b ); 1038 bExp = extractFloat32Exp( b ); 1039 expDiff = aExp - bExp; 1040 aSig <<= 6; 1041 bSig <<= 6; 1042 if ( 0 < expDiff ) { 1043 if ( aExp == 0xFF ) { 1044 if ( aSig ) return propagateFloat32NaN( a, b ); 1045 return a; 1046 } 1047 if ( bExp == 0 ) { 1048 --expDiff; 1049 } 1050 else { 1051 bSig |= 0x20000000; 1052 } 1053 shift32RightJamming( bSig, expDiff, &bSig ); 1054 zExp = aExp; 1055 } 1056 else if ( expDiff < 0 ) { 1057 if ( bExp == 0xFF ) { 1058 if ( bSig ) return propagateFloat32NaN( a, b ); 1059 return packFloat32( zSign, 0xFF, 0 ); 1060 } 1061 if ( aExp == 0 ) { 1062 ++expDiff; 1063 } 1064 else { 1065 aSig |= 0x20000000; 1066 } 1067 shift32RightJamming( aSig, - expDiff, &aSig ); 1068 zExp = bExp; 1069 } 1070 else { 1071 if ( aExp == 0xFF ) { 1072 if ( aSig | bSig ) return propagateFloat32NaN( a, b ); 1073 return a; 1074 } 1075 if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 ); 1076 zSig = 0x40000000 + aSig + bSig; 1077 zExp = aExp; 1078 goto roundAndPack; 1079 } 1080 aSig |= 0x20000000; 1081 zSig = ( aSig + bSig )<<1; 1082 --zExp; 1083 if ( (sbits32) zSig < 0 ) { 1084 zSig = aSig + bSig; 1085 ++zExp; 1086 } 1087 roundAndPack: 1088 return roundAndPackFloat32( roundData, zSign, zExp, zSig ); 1089 1090 } 1091 1092 /* 1093 ------------------------------------------------------------------------------- 1094 Returns the result of subtracting the absolute values of the single- 1095 precision floating-point values `a' and `b'. If `zSign' is true, the 1096 difference is negated before being returned. `zSign' is ignored if the 1097 result is a NaN. The subtraction is performed according to the IEC/IEEE 1098 Standard for Binary Floating-point Arithmetic. 1099 ------------------------------------------------------------------------------- 1100 */ 1101 static float32 subFloat32Sigs( struct roundingData *roundData, float32 a, float32 b, flag zSign ) 1102 { 1103 int16 aExp, bExp, zExp; 1104 bits32 aSig, bSig, zSig; 1105 int16 expDiff; 1106 1107 aSig = extractFloat32Frac( a ); 1108 aExp = extractFloat32Exp( a ); 1109 bSig = extractFloat32Frac( b ); 1110 bExp = extractFloat32Exp( b ); 1111 expDiff = aExp - bExp; 1112 aSig <<= 7; 1113 bSig <<= 7; 1114 if ( 0 < expDiff ) goto aExpBigger; 1115 if ( expDiff < 0 ) goto bExpBigger; 1116 if ( aExp == 0xFF ) { 1117 if ( aSig | bSig ) return propagateFloat32NaN( a, b ); 1118 roundData->exception |= float_flag_invalid; 1119 return float32_default_nan; 1120 } 1121 if ( aExp == 0 ) { 1122 aExp = 1; 1123 bExp = 1; 1124 } 1125 if ( bSig < aSig ) goto aBigger; 1126 if ( aSig < bSig ) goto bBigger; 1127 return packFloat32( roundData->mode == float_round_down, 0, 0 ); 1128 bExpBigger: 1129 if ( bExp == 0xFF ) { 1130 if ( bSig ) return propagateFloat32NaN( a, b ); 1131 return packFloat32( zSign ^ 1, 0xFF, 0 ); 1132 } 1133 if ( aExp == 0 ) { 1134 ++expDiff; 1135 } 1136 else { 1137 aSig |= 0x40000000; 1138 } 1139 shift32RightJamming( aSig, - expDiff, &aSig ); 1140 bSig |= 0x40000000; 1141 bBigger: 1142 zSig = bSig - aSig; 1143 zExp = bExp; 1144 zSign ^= 1; 1145 goto normalizeRoundAndPack; 1146 aExpBigger: 1147 if ( aExp == 0xFF ) { 1148 if ( aSig ) return propagateFloat32NaN( a, b ); 1149 return a; 1150 } 1151 if ( bExp == 0 ) { 1152 --expDiff; 1153 } 1154 else { 1155 bSig |= 0x40000000; 1156 } 1157 shift32RightJamming( bSig, expDiff, &bSig ); 1158 aSig |= 0x40000000; 1159 aBigger: 1160 zSig = aSig - bSig; 1161 zExp = aExp; 1162 normalizeRoundAndPack: 1163 --zExp; 1164 return normalizeRoundAndPackFloat32( roundData, zSign, zExp, zSig ); 1165 1166 } 1167 1168 /* 1169 ------------------------------------------------------------------------------- 1170 Returns the result of adding the single-precision floating-point values `a' 1171 and `b'. The operation is performed according to the IEC/IEEE Standard for 1172 Binary Floating-point Arithmetic. 1173 ------------------------------------------------------------------------------- 1174 */ 1175 float32 float32_add( struct roundingData *roundData, float32 a, float32 b ) 1176 { 1177 flag aSign, bSign; 1178 1179 aSign = extractFloat32Sign( a ); 1180 bSign = extractFloat32Sign( b ); 1181 if ( aSign == bSign ) { 1182 return addFloat32Sigs( roundData, a, b, aSign ); 1183 } 1184 else { 1185 return subFloat32Sigs( roundData, a, b, aSign ); 1186 } 1187 1188 } 1189 1190 /* 1191 ------------------------------------------------------------------------------- 1192 Returns the result of subtracting the single-precision floating-point values 1193 `a' and `b'. The operation is performed according to the IEC/IEEE Standard 1194 for Binary Floating-point Arithmetic. 1195 ------------------------------------------------------------------------------- 1196 */ 1197 float32 float32_sub( struct roundingData *roundData, float32 a, float32 b ) 1198 { 1199 flag aSign, bSign; 1200 1201 aSign = extractFloat32Sign( a ); 1202 bSign = extractFloat32Sign( b ); 1203 if ( aSign == bSign ) { 1204 return subFloat32Sigs( roundData, a, b, aSign ); 1205 } 1206 else { 1207 return addFloat32Sigs( roundData, a, b, aSign ); 1208 } 1209 1210 } 1211 1212 /* 1213 ------------------------------------------------------------------------------- 1214 Returns the result of multiplying the single-precision floating-point values 1215 `a' and `b'. The operation is performed according to the IEC/IEEE Standard 1216 for Binary Floating-point Arithmetic. 1217 ------------------------------------------------------------------------------- 1218 */ 1219 float32 float32_mul( struct roundingData *roundData, float32 a, float32 b ) 1220 { 1221 flag aSign, bSign, zSign; 1222 int16 aExp, bExp, zExp; 1223 bits32 aSig, bSig; 1224 bits64 zSig64; 1225 bits32 zSig; 1226 1227 aSig = extractFloat32Frac( a ); 1228 aExp = extractFloat32Exp( a ); 1229 aSign = extractFloat32Sign( a ); 1230 bSig = extractFloat32Frac( b ); 1231 bExp = extractFloat32Exp( b ); 1232 bSign = extractFloat32Sign( b ); 1233 zSign = aSign ^ bSign; 1234 if ( aExp == 0xFF ) { 1235 if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) { 1236 return propagateFloat32NaN( a, b ); 1237 } 1238 if ( ( bExp | bSig ) == 0 ) { 1239 roundData->exception |= float_flag_invalid; 1240 return float32_default_nan; 1241 } 1242 return packFloat32( zSign, 0xFF, 0 ); 1243 } 1244 if ( bExp == 0xFF ) { 1245 if ( bSig ) return propagateFloat32NaN( a, b ); 1246 if ( ( aExp | aSig ) == 0 ) { 1247 roundData->exception |= float_flag_invalid; 1248 return float32_default_nan; 1249 } 1250 return packFloat32( zSign, 0xFF, 0 ); 1251 } 1252 if ( aExp == 0 ) { 1253 if ( aSig == 0 ) return packFloat32( zSign, 0, 0 ); 1254 normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1255 } 1256 if ( bExp == 0 ) { 1257 if ( bSig == 0 ) return packFloat32( zSign, 0, 0 ); 1258 normalizeFloat32Subnormal( bSig, &bExp, &bSig ); 1259 } 1260 zExp = aExp + bExp - 0x7F; 1261 aSig = ( aSig | 0x00800000 )<<7; 1262 bSig = ( bSig | 0x00800000 )<<8; 1263 shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 ); 1264 zSig = zSig64; 1265 if ( 0 <= (sbits32) ( zSig<<1 ) ) { 1266 zSig <<= 1; 1267 --zExp; 1268 } 1269 return roundAndPackFloat32( roundData, zSign, zExp, zSig ); 1270 1271 } 1272 1273 /* 1274 ------------------------------------------------------------------------------- 1275 Returns the result of dividing the single-precision floating-point value `a' 1276 by the corresponding value `b'. The operation is performed according to the 1277 IEC/IEEE Standard for Binary Floating-point Arithmetic. 1278 ------------------------------------------------------------------------------- 1279 */ 1280 float32 float32_div( struct roundingData *roundData, float32 a, float32 b ) 1281 { 1282 flag aSign, bSign, zSign; 1283 int16 aExp, bExp, zExp; 1284 bits32 aSig, bSig, zSig; 1285 1286 aSig = extractFloat32Frac( a ); 1287 aExp = extractFloat32Exp( a ); 1288 aSign = extractFloat32Sign( a ); 1289 bSig = extractFloat32Frac( b ); 1290 bExp = extractFloat32Exp( b ); 1291 bSign = extractFloat32Sign( b ); 1292 zSign = aSign ^ bSign; 1293 if ( aExp == 0xFF ) { 1294 if ( aSig ) return propagateFloat32NaN( a, b ); 1295 if ( bExp == 0xFF ) { 1296 if ( bSig ) return propagateFloat32NaN( a, b ); 1297 roundData->exception |= float_flag_invalid; 1298 return float32_default_nan; 1299 } 1300 return packFloat32( zSign, 0xFF, 0 ); 1301 } 1302 if ( bExp == 0xFF ) { 1303 if ( bSig ) return propagateFloat32NaN( a, b ); 1304 return packFloat32( zSign, 0, 0 ); 1305 } 1306 if ( bExp == 0 ) { 1307 if ( bSig == 0 ) { 1308 if ( ( aExp | aSig ) == 0 ) { 1309 roundData->exception |= float_flag_invalid; 1310 return float32_default_nan; 1311 } 1312 roundData->exception |= float_flag_divbyzero; 1313 return packFloat32( zSign, 0xFF, 0 ); 1314 } 1315 normalizeFloat32Subnormal( bSig, &bExp, &bSig ); 1316 } 1317 if ( aExp == 0 ) { 1318 if ( aSig == 0 ) return packFloat32( zSign, 0, 0 ); 1319 normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1320 } 1321 zExp = aExp - bExp + 0x7D; 1322 aSig = ( aSig | 0x00800000 )<<7; 1323 bSig = ( bSig | 0x00800000 )<<8; 1324 if ( bSig <= ( aSig + aSig ) ) { 1325 aSig >>= 1; 1326 ++zExp; 1327 } 1328 { 1329 bits64 tmp = ( (bits64) aSig )<<32; 1330 do_div( tmp, bSig ); 1331 zSig = tmp; 1332 } 1333 if ( ( zSig & 0x3F ) == 0 ) { 1334 zSig |= ( ( (bits64) bSig ) * zSig != ( (bits64) aSig )<<32 ); 1335 } 1336 return roundAndPackFloat32( roundData, zSign, zExp, zSig ); 1337 1338 } 1339 1340 /* 1341 ------------------------------------------------------------------------------- 1342 Returns the remainder of the single-precision floating-point value `a' 1343 with respect to the corresponding value `b'. The operation is performed 1344 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 1345 ------------------------------------------------------------------------------- 1346 */ 1347 float32 float32_rem( struct roundingData *roundData, float32 a, float32 b ) 1348 { 1349 flag aSign, bSign, zSign; 1350 int16 aExp, bExp, expDiff; 1351 bits32 aSig, bSig; 1352 bits32 q; 1353 bits64 aSig64, bSig64, q64; 1354 bits32 alternateASig; 1355 sbits32 sigMean; 1356 1357 aSig = extractFloat32Frac( a ); 1358 aExp = extractFloat32Exp( a ); 1359 aSign = extractFloat32Sign( a ); 1360 bSig = extractFloat32Frac( b ); 1361 bExp = extractFloat32Exp( b ); 1362 bSign = extractFloat32Sign( b ); 1363 if ( aExp == 0xFF ) { 1364 if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) { 1365 return propagateFloat32NaN( a, b ); 1366 } 1367 roundData->exception |= float_flag_invalid; 1368 return float32_default_nan; 1369 } 1370 if ( bExp == 0xFF ) { 1371 if ( bSig ) return propagateFloat32NaN( a, b ); 1372 return a; 1373 } 1374 if ( bExp == 0 ) { 1375 if ( bSig == 0 ) { 1376 roundData->exception |= float_flag_invalid; 1377 return float32_default_nan; 1378 } 1379 normalizeFloat32Subnormal( bSig, &bExp, &bSig ); 1380 } 1381 if ( aExp == 0 ) { 1382 if ( aSig == 0 ) return a; 1383 normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1384 } 1385 expDiff = aExp - bExp; 1386 aSig |= 0x00800000; 1387 bSig |= 0x00800000; 1388 if ( expDiff < 32 ) { 1389 aSig <<= 8; 1390 bSig <<= 8; 1391 if ( expDiff < 0 ) { 1392 if ( expDiff < -1 ) return a; 1393 aSig >>= 1; 1394 } 1395 q = ( bSig <= aSig ); 1396 if ( q ) aSig -= bSig; 1397 if ( 0 < expDiff ) { 1398 bits64 tmp = ( (bits64) aSig )<<32; 1399 do_div( tmp, bSig ); 1400 q = tmp; 1401 q >>= 32 - expDiff; 1402 bSig >>= 2; 1403 aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q; 1404 } 1405 else { 1406 aSig >>= 2; 1407 bSig >>= 2; 1408 } 1409 } 1410 else { 1411 if ( bSig <= aSig ) aSig -= bSig; 1412 aSig64 = ( (bits64) aSig )<<40; 1413 bSig64 = ( (bits64) bSig )<<40; 1414 expDiff -= 64; 1415 while ( 0 < expDiff ) { 1416 q64 = estimateDiv128To64( aSig64, 0, bSig64 ); 1417 q64 = ( 2 < q64 ) ? q64 - 2 : 0; 1418 aSig64 = - ( ( bSig * q64 )<<38 ); 1419 expDiff -= 62; 1420 } 1421 expDiff += 64; 1422 q64 = estimateDiv128To64( aSig64, 0, bSig64 ); 1423 q64 = ( 2 < q64 ) ? q64 - 2 : 0; 1424 q = q64>>( 64 - expDiff ); 1425 bSig <<= 6; 1426 aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q; 1427 } 1428 do { 1429 alternateASig = aSig; 1430 ++q; 1431 aSig -= bSig; 1432 } while ( 0 <= (sbits32) aSig ); 1433 sigMean = aSig + alternateASig; 1434 if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) { 1435 aSig = alternateASig; 1436 } 1437 zSign = ( (sbits32) aSig < 0 ); 1438 if ( zSign ) aSig = - aSig; 1439 return normalizeRoundAndPackFloat32( roundData, aSign ^ zSign, bExp, aSig ); 1440 1441 } 1442 1443 /* 1444 ------------------------------------------------------------------------------- 1445 Returns the square root of the single-precision floating-point value `a'. 1446 The operation is performed according to the IEC/IEEE Standard for Binary 1447 Floating-point Arithmetic. 1448 ------------------------------------------------------------------------------- 1449 */ 1450 float32 float32_sqrt( struct roundingData *roundData, float32 a ) 1451 { 1452 flag aSign; 1453 int16 aExp, zExp; 1454 bits32 aSig, zSig; 1455 bits64 rem, term; 1456 1457 aSig = extractFloat32Frac( a ); 1458 aExp = extractFloat32Exp( a ); 1459 aSign = extractFloat32Sign( a ); 1460 if ( aExp == 0xFF ) { 1461 if ( aSig ) return propagateFloat32NaN( a, 0 ); 1462 if ( ! aSign ) return a; 1463 roundData->exception |= float_flag_invalid; 1464 return float32_default_nan; 1465 } 1466 if ( aSign ) { 1467 if ( ( aExp | aSig ) == 0 ) return a; 1468 roundData->exception |= float_flag_invalid; 1469 return float32_default_nan; 1470 } 1471 if ( aExp == 0 ) { 1472 if ( aSig == 0 ) return 0; 1473 normalizeFloat32Subnormal( aSig, &aExp, &aSig ); 1474 } 1475 zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E; 1476 aSig = ( aSig | 0x00800000 )<<8; 1477 zSig = estimateSqrt32( aExp, aSig ) + 2; 1478 if ( ( zSig & 0x7F ) <= 5 ) { 1479 if ( zSig < 2 ) { 1480 zSig = 0xFFFFFFFF; 1481 } 1482 else { 1483 aSig >>= aExp & 1; 1484 term = ( (bits64) zSig ) * zSig; 1485 rem = ( ( (bits64) aSig )<<32 ) - term; 1486 while ( (sbits64) rem < 0 ) { 1487 --zSig; 1488 rem += ( ( (bits64) zSig )<<1 ) | 1; 1489 } 1490 zSig |= ( rem != 0 ); 1491 } 1492 } 1493 shift32RightJamming( zSig, 1, &zSig ); 1494 return roundAndPackFloat32( roundData, 0, zExp, zSig ); 1495 1496 } 1497 1498 /* 1499 ------------------------------------------------------------------------------- 1500 Returns 1 if the single-precision floating-point value `a' is equal to the 1501 corresponding value `b', and 0 otherwise. The comparison is performed 1502 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 1503 ------------------------------------------------------------------------------- 1504 */ 1505 flag float32_eq( float32 a, float32 b ) 1506 { 1507 1508 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1509 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1510 ) { 1511 if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) { 1512 float_raise( float_flag_invalid ); 1513 } 1514 return 0; 1515 } 1516 return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1517 1518 } 1519 1520 /* 1521 ------------------------------------------------------------------------------- 1522 Returns 1 if the single-precision floating-point value `a' is less than or 1523 equal to the corresponding value `b', and 0 otherwise. The comparison is 1524 performed according to the IEC/IEEE Standard for Binary Floating-point 1525 Arithmetic. 1526 ------------------------------------------------------------------------------- 1527 */ 1528 flag float32_le( float32 a, float32 b ) 1529 { 1530 flag aSign, bSign; 1531 1532 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1533 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1534 ) { 1535 float_raise( float_flag_invalid ); 1536 return 0; 1537 } 1538 aSign = extractFloat32Sign( a ); 1539 bSign = extractFloat32Sign( b ); 1540 if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1541 return ( a == b ) || ( aSign ^ ( a < b ) ); 1542 1543 } 1544 1545 /* 1546 ------------------------------------------------------------------------------- 1547 Returns 1 if the single-precision floating-point value `a' is less than 1548 the corresponding value `b', and 0 otherwise. The comparison is performed 1549 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 1550 ------------------------------------------------------------------------------- 1551 */ 1552 flag float32_lt( float32 a, float32 b ) 1553 { 1554 flag aSign, bSign; 1555 1556 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1557 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1558 ) { 1559 float_raise( float_flag_invalid ); 1560 return 0; 1561 } 1562 aSign = extractFloat32Sign( a ); 1563 bSign = extractFloat32Sign( b ); 1564 if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 ); 1565 return ( a != b ) && ( aSign ^ ( a < b ) ); 1566 1567 } 1568 1569 /* 1570 ------------------------------------------------------------------------------- 1571 Returns 1 if the single-precision floating-point value `a' is equal to the 1572 corresponding value `b', and 0 otherwise. The invalid exception is raised 1573 if either operand is a NaN. Otherwise, the comparison is performed 1574 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 1575 ------------------------------------------------------------------------------- 1576 */ 1577 flag float32_eq_signaling( float32 a, float32 b ) 1578 { 1579 1580 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1581 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1582 ) { 1583 float_raise( float_flag_invalid ); 1584 return 0; 1585 } 1586 return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1587 1588 } 1589 1590 /* 1591 ------------------------------------------------------------------------------- 1592 Returns 1 if the single-precision floating-point value `a' is less than or 1593 equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not 1594 cause an exception. Otherwise, the comparison is performed according to the 1595 IEC/IEEE Standard for Binary Floating-point Arithmetic. 1596 ------------------------------------------------------------------------------- 1597 */ 1598 flag float32_le_quiet( float32 a, float32 b ) 1599 { 1600 flag aSign, bSign; 1601 //int16 aExp, bExp; 1602 1603 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1604 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1605 ) { 1606 /* Do nothing, even if NaN as we're quiet */ 1607 return 0; 1608 } 1609 aSign = extractFloat32Sign( a ); 1610 bSign = extractFloat32Sign( b ); 1611 if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 ); 1612 return ( a == b ) || ( aSign ^ ( a < b ) ); 1613 1614 } 1615 1616 /* 1617 ------------------------------------------------------------------------------- 1618 Returns 1 if the single-precision floating-point value `a' is less than 1619 the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an 1620 exception. Otherwise, the comparison is performed according to the IEC/IEEE 1621 Standard for Binary Floating-point Arithmetic. 1622 ------------------------------------------------------------------------------- 1623 */ 1624 flag float32_lt_quiet( float32 a, float32 b ) 1625 { 1626 flag aSign, bSign; 1627 1628 if ( ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) ) 1629 || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) ) 1630 ) { 1631 /* Do nothing, even if NaN as we're quiet */ 1632 return 0; 1633 } 1634 aSign = extractFloat32Sign( a ); 1635 bSign = extractFloat32Sign( b ); 1636 if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 ); 1637 return ( a != b ) && ( aSign ^ ( a < b ) ); 1638 1639 } 1640 1641 /* 1642 ------------------------------------------------------------------------------- 1643 Returns the result of converting the double-precision floating-point value 1644 `a' to the 32-bit two's complement integer format. The conversion is 1645 performed according to the IEC/IEEE Standard for Binary Floating-point 1646 Arithmetic---which means in particular that the conversion is rounded 1647 according to the current rounding mode. If `a' is a NaN, the largest 1648 positive integer is returned. Otherwise, if the conversion overflows, the 1649 largest integer with the same sign as `a' is returned. 1650 ------------------------------------------------------------------------------- 1651 */ 1652 int32 float64_to_int32( struct roundingData *roundData, float64 a ) 1653 { 1654 flag aSign; 1655 int16 aExp, shiftCount; 1656 bits64 aSig; 1657 1658 aSig = extractFloat64Frac( a ); 1659 aExp = extractFloat64Exp( a ); 1660 aSign = extractFloat64Sign( a ); 1661 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0; 1662 if ( aExp ) aSig |= LIT64( 0x0010000000000000 ); 1663 shiftCount = 0x42C - aExp; 1664 if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig ); 1665 return roundAndPackInt32( roundData, aSign, aSig ); 1666 1667 } 1668 1669 /* 1670 ------------------------------------------------------------------------------- 1671 Returns the result of converting the double-precision floating-point value 1672 `a' to the 32-bit two's complement integer format. The conversion is 1673 performed according to the IEC/IEEE Standard for Binary Floating-point 1674 Arithmetic, except that the conversion is always rounded toward zero. If 1675 `a' is a NaN, the largest positive integer is returned. Otherwise, if the 1676 conversion overflows, the largest integer with the same sign as `a' is 1677 returned. 1678 ------------------------------------------------------------------------------- 1679 */ 1680 int32 float64_to_int32_round_to_zero( float64 a ) 1681 { 1682 flag aSign; 1683 int16 aExp, shiftCount; 1684 bits64 aSig, savedASig; 1685 int32 z; 1686 1687 aSig = extractFloat64Frac( a ); 1688 aExp = extractFloat64Exp( a ); 1689 aSign = extractFloat64Sign( a ); 1690 shiftCount = 0x433 - aExp; 1691 if ( shiftCount < 21 ) { 1692 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0; 1693 goto invalid; 1694 } 1695 else if ( 52 < shiftCount ) { 1696 if ( aExp || aSig ) float_raise( float_flag_inexact ); 1697 return 0; 1698 } 1699 aSig |= LIT64( 0x0010000000000000 ); 1700 savedASig = aSig; 1701 aSig >>= shiftCount; 1702 z = aSig; 1703 if ( aSign ) z = - z; 1704 if ( ( z < 0 ) ^ aSign ) { 1705 invalid: 1706 float_raise( float_flag_invalid ); 1707 return aSign ? 0x80000000 : 0x7FFFFFFF; 1708 } 1709 if ( ( aSig<<shiftCount ) != savedASig ) { 1710 float_raise( float_flag_inexact ); 1711 } 1712 return z; 1713 1714 } 1715 1716 /* 1717 ------------------------------------------------------------------------------- 1718 Returns the result of converting the double-precision floating-point value 1719 `a' to the 32-bit two's complement unsigned integer format. The conversion 1720 is performed according to the IEC/IEEE Standard for Binary Floating-point 1721 Arithmetic---which means in particular that the conversion is rounded 1722 according to the current rounding mode. If `a' is a NaN, the largest 1723 positive integer is returned. Otherwise, if the conversion overflows, the 1724 largest positive integer is returned. 1725 ------------------------------------------------------------------------------- 1726 */ 1727 int32 float64_to_uint32( struct roundingData *roundData, float64 a ) 1728 { 1729 flag aSign; 1730 int16 aExp, shiftCount; 1731 bits64 aSig; 1732 1733 aSig = extractFloat64Frac( a ); 1734 aExp = extractFloat64Exp( a ); 1735 aSign = 0; //extractFloat64Sign( a ); 1736 //if ( ( aExp == 0x7FF ) && aSig ) aSign = 0; 1737 if ( aExp ) aSig |= LIT64( 0x0010000000000000 ); 1738 shiftCount = 0x42C - aExp; 1739 if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig ); 1740 return roundAndPackInt32( roundData, aSign, aSig ); 1741 } 1742 1743 /* 1744 ------------------------------------------------------------------------------- 1745 Returns the result of converting the double-precision floating-point value 1746 `a' to the 32-bit two's complement integer format. The conversion is 1747 performed according to the IEC/IEEE Standard for Binary Floating-point 1748 Arithmetic, except that the conversion is always rounded toward zero. If 1749 `a' is a NaN, the largest positive integer is returned. Otherwise, if the 1750 conversion overflows, the largest positive integer is returned. 1751 ------------------------------------------------------------------------------- 1752 */ 1753 int32 float64_to_uint32_round_to_zero( float64 a ) 1754 { 1755 flag aSign; 1756 int16 aExp, shiftCount; 1757 bits64 aSig, savedASig; 1758 int32 z; 1759 1760 aSig = extractFloat64Frac( a ); 1761 aExp = extractFloat64Exp( a ); 1762 aSign = extractFloat64Sign( a ); 1763 shiftCount = 0x433 - aExp; 1764 if ( shiftCount < 21 ) { 1765 if ( ( aExp == 0x7FF ) && aSig ) aSign = 0; 1766 goto invalid; 1767 } 1768 else if ( 52 < shiftCount ) { 1769 if ( aExp || aSig ) float_raise( float_flag_inexact ); 1770 return 0; 1771 } 1772 aSig |= LIT64( 0x0010000000000000 ); 1773 savedASig = aSig; 1774 aSig >>= shiftCount; 1775 z = aSig; 1776 if ( aSign ) z = - z; 1777 if ( ( z < 0 ) ^ aSign ) { 1778 invalid: 1779 float_raise( float_flag_invalid ); 1780 return aSign ? 0x80000000 : 0x7FFFFFFF; 1781 } 1782 if ( ( aSig<<shiftCount ) != savedASig ) { 1783 float_raise( float_flag_inexact ); 1784 } 1785 return z; 1786 } 1787 1788 /* 1789 ------------------------------------------------------------------------------- 1790 Returns the result of converting the double-precision floating-point value 1791 `a' to the single-precision floating-point format. The conversion is 1792 performed according to the IEC/IEEE Standard for Binary Floating-point 1793 Arithmetic. 1794 ------------------------------------------------------------------------------- 1795 */ 1796 float32 float64_to_float32( struct roundingData *roundData, float64 a ) 1797 { 1798 flag aSign; 1799 int16 aExp; 1800 bits64 aSig; 1801 bits32 zSig; 1802 1803 aSig = extractFloat64Frac( a ); 1804 aExp = extractFloat64Exp( a ); 1805 aSign = extractFloat64Sign( a ); 1806 if ( aExp == 0x7FF ) { 1807 if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a ) ); 1808 return packFloat32( aSign, 0xFF, 0 ); 1809 } 1810 shift64RightJamming( aSig, 22, &aSig ); 1811 zSig = aSig; 1812 if ( aExp || zSig ) { 1813 zSig |= 0x40000000; 1814 aExp -= 0x381; 1815 } 1816 return roundAndPackFloat32( roundData, aSign, aExp, zSig ); 1817 1818 } 1819 1820 #ifdef FLOATX80 1821 1822 /* 1823 ------------------------------------------------------------------------------- 1824 Returns the result of converting the double-precision floating-point value 1825 `a' to the extended double-precision floating-point format. The conversion 1826 is performed according to the IEC/IEEE Standard for Binary Floating-point 1827 Arithmetic. 1828 ------------------------------------------------------------------------------- 1829 */ 1830 floatx80 float64_to_floatx80( float64 a ) 1831 { 1832 flag aSign; 1833 int16 aExp; 1834 bits64 aSig; 1835 1836 aSig = extractFloat64Frac( a ); 1837 aExp = extractFloat64Exp( a ); 1838 aSign = extractFloat64Sign( a ); 1839 if ( aExp == 0x7FF ) { 1840 if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a ) ); 1841 return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 1842 } 1843 if ( aExp == 0 ) { 1844 if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 ); 1845 normalizeFloat64Subnormal( aSig, &aExp, &aSig ); 1846 } 1847 return 1848 packFloatx80( 1849 aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 ); 1850 1851 } 1852 1853 #endif 1854 1855 /* 1856 ------------------------------------------------------------------------------- 1857 Rounds the double-precision floating-point value `a' to an integer, and 1858 returns the result as a double-precision floating-point value. The 1859 operation is performed according to the IEC/IEEE Standard for Binary 1860 Floating-point Arithmetic. 1861 ------------------------------------------------------------------------------- 1862 */ 1863 float64 float64_round_to_int( struct roundingData *roundData, float64 a ) 1864 { 1865 flag aSign; 1866 int16 aExp; 1867 bits64 lastBitMask, roundBitsMask; 1868 int8 roundingMode; 1869 float64 z; 1870 1871 aExp = extractFloat64Exp( a ); 1872 if ( 0x433 <= aExp ) { 1873 if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) { 1874 return propagateFloat64NaN( a, a ); 1875 } 1876 return a; 1877 } 1878 if ( aExp <= 0x3FE ) { 1879 if ( (bits64) ( a<<1 ) == 0 ) return a; 1880 roundData->exception |= float_flag_inexact; 1881 aSign = extractFloat64Sign( a ); 1882 switch ( roundData->mode ) { 1883 case float_round_nearest_even: 1884 if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) { 1885 return packFloat64( aSign, 0x3FF, 0 ); 1886 } 1887 break; 1888 case float_round_down: 1889 return aSign ? LIT64( 0xBFF0000000000000 ) : 0; 1890 case float_round_up: 1891 return 1892 aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 ); 1893 } 1894 return packFloat64( aSign, 0, 0 ); 1895 } 1896 lastBitMask = 1; 1897 lastBitMask <<= 0x433 - aExp; 1898 roundBitsMask = lastBitMask - 1; 1899 z = a; 1900 roundingMode = roundData->mode; 1901 if ( roundingMode == float_round_nearest_even ) { 1902 z += lastBitMask>>1; 1903 if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask; 1904 } 1905 else if ( roundingMode != float_round_to_zero ) { 1906 if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) { 1907 z += roundBitsMask; 1908 } 1909 } 1910 z &= ~ roundBitsMask; 1911 if ( z != a ) roundData->exception |= float_flag_inexact; 1912 return z; 1913 1914 } 1915 1916 /* 1917 ------------------------------------------------------------------------------- 1918 Returns the result of adding the absolute values of the double-precision 1919 floating-point values `a' and `b'. If `zSign' is true, the sum is negated 1920 before being returned. `zSign' is ignored if the result is a NaN. The 1921 addition is performed according to the IEC/IEEE Standard for Binary 1922 Floating-point Arithmetic. 1923 ------------------------------------------------------------------------------- 1924 */ 1925 static float64 addFloat64Sigs( struct roundingData *roundData, float64 a, float64 b, flag zSign ) 1926 { 1927 int16 aExp, bExp, zExp; 1928 bits64 aSig, bSig, zSig; 1929 int16 expDiff; 1930 1931 aSig = extractFloat64Frac( a ); 1932 aExp = extractFloat64Exp( a ); 1933 bSig = extractFloat64Frac( b ); 1934 bExp = extractFloat64Exp( b ); 1935 expDiff = aExp - bExp; 1936 aSig <<= 9; 1937 bSig <<= 9; 1938 if ( 0 < expDiff ) { 1939 if ( aExp == 0x7FF ) { 1940 if ( aSig ) return propagateFloat64NaN( a, b ); 1941 return a; 1942 } 1943 if ( bExp == 0 ) { 1944 --expDiff; 1945 } 1946 else { 1947 bSig |= LIT64( 0x2000000000000000 ); 1948 } 1949 shift64RightJamming( bSig, expDiff, &bSig ); 1950 zExp = aExp; 1951 } 1952 else if ( expDiff < 0 ) { 1953 if ( bExp == 0x7FF ) { 1954 if ( bSig ) return propagateFloat64NaN( a, b ); 1955 return packFloat64( zSign, 0x7FF, 0 ); 1956 } 1957 if ( aExp == 0 ) { 1958 ++expDiff; 1959 } 1960 else { 1961 aSig |= LIT64( 0x2000000000000000 ); 1962 } 1963 shift64RightJamming( aSig, - expDiff, &aSig ); 1964 zExp = bExp; 1965 } 1966 else { 1967 if ( aExp == 0x7FF ) { 1968 if ( aSig | bSig ) return propagateFloat64NaN( a, b ); 1969 return a; 1970 } 1971 if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 ); 1972 zSig = LIT64( 0x4000000000000000 ) + aSig + bSig; 1973 zExp = aExp; 1974 goto roundAndPack; 1975 } 1976 aSig |= LIT64( 0x2000000000000000 ); 1977 zSig = ( aSig + bSig )<<1; 1978 --zExp; 1979 if ( (sbits64) zSig < 0 ) { 1980 zSig = aSig + bSig; 1981 ++zExp; 1982 } 1983 roundAndPack: 1984 return roundAndPackFloat64( roundData, zSign, zExp, zSig ); 1985 1986 } 1987 1988 /* 1989 ------------------------------------------------------------------------------- 1990 Returns the result of subtracting the absolute values of the double- 1991 precision floating-point values `a' and `b'. If `zSign' is true, the 1992 difference is negated before being returned. `zSign' is ignored if the 1993 result is a NaN. The subtraction is performed according to the IEC/IEEE 1994 Standard for Binary Floating-point Arithmetic. 1995 ------------------------------------------------------------------------------- 1996 */ 1997 static float64 subFloat64Sigs( struct roundingData *roundData, float64 a, float64 b, flag zSign ) 1998 { 1999 int16 aExp, bExp, zExp; 2000 bits64 aSig, bSig, zSig; 2001 int16 expDiff; 2002 2003 aSig = extractFloat64Frac( a ); 2004 aExp = extractFloat64Exp( a ); 2005 bSig = extractFloat64Frac( b ); 2006 bExp = extractFloat64Exp( b ); 2007 expDiff = aExp - bExp; 2008 aSig <<= 10; 2009 bSig <<= 10; 2010 if ( 0 < expDiff ) goto aExpBigger; 2011 if ( expDiff < 0 ) goto bExpBigger; 2012 if ( aExp == 0x7FF ) { 2013 if ( aSig | bSig ) return propagateFloat64NaN( a, b ); 2014 roundData->exception |= float_flag_invalid; 2015 return float64_default_nan; 2016 } 2017 if ( aExp == 0 ) { 2018 aExp = 1; 2019 bExp = 1; 2020 } 2021 if ( bSig < aSig ) goto aBigger; 2022 if ( aSig < bSig ) goto bBigger; 2023 return packFloat64( roundData->mode == float_round_down, 0, 0 ); 2024 bExpBigger: 2025 if ( bExp == 0x7FF ) { 2026 if ( bSig ) return propagateFloat64NaN( a, b ); 2027 return packFloat64( zSign ^ 1, 0x7FF, 0 ); 2028 } 2029 if ( aExp == 0 ) { 2030 ++expDiff; 2031 } 2032 else { 2033 aSig |= LIT64( 0x4000000000000000 ); 2034 } 2035 shift64RightJamming( aSig, - expDiff, &aSig ); 2036 bSig |= LIT64( 0x4000000000000000 ); 2037 bBigger: 2038 zSig = bSig - aSig; 2039 zExp = bExp; 2040 zSign ^= 1; 2041 goto normalizeRoundAndPack; 2042 aExpBigger: 2043 if ( aExp == 0x7FF ) { 2044 if ( aSig ) return propagateFloat64NaN( a, b ); 2045 return a; 2046 } 2047 if ( bExp == 0 ) { 2048 --expDiff; 2049 } 2050 else { 2051 bSig |= LIT64( 0x4000000000000000 ); 2052 } 2053 shift64RightJamming( bSig, expDiff, &bSig ); 2054 aSig |= LIT64( 0x4000000000000000 ); 2055 aBigger: 2056 zSig = aSig - bSig; 2057 zExp = aExp; 2058 normalizeRoundAndPack: 2059 --zExp; 2060 return normalizeRoundAndPackFloat64( roundData, zSign, zExp, zSig ); 2061 2062 } 2063 2064 /* 2065 ------------------------------------------------------------------------------- 2066 Returns the result of adding the double-precision floating-point values `a' 2067 and `b'. The operation is performed according to the IEC/IEEE Standard for 2068 Binary Floating-point Arithmetic. 2069 ------------------------------------------------------------------------------- 2070 */ 2071 float64 float64_add( struct roundingData *roundData, float64 a, float64 b ) 2072 { 2073 flag aSign, bSign; 2074 2075 aSign = extractFloat64Sign( a ); 2076 bSign = extractFloat64Sign( b ); 2077 if ( aSign == bSign ) { 2078 return addFloat64Sigs( roundData, a, b, aSign ); 2079 } 2080 else { 2081 return subFloat64Sigs( roundData, a, b, aSign ); 2082 } 2083 2084 } 2085 2086 /* 2087 ------------------------------------------------------------------------------- 2088 Returns the result of subtracting the double-precision floating-point values 2089 `a' and `b'. The operation is performed according to the IEC/IEEE Standard 2090 for Binary Floating-point Arithmetic. 2091 ------------------------------------------------------------------------------- 2092 */ 2093 float64 float64_sub( struct roundingData *roundData, float64 a, float64 b ) 2094 { 2095 flag aSign, bSign; 2096 2097 aSign = extractFloat64Sign( a ); 2098 bSign = extractFloat64Sign( b ); 2099 if ( aSign == bSign ) { 2100 return subFloat64Sigs( roundData, a, b, aSign ); 2101 } 2102 else { 2103 return addFloat64Sigs( roundData, a, b, aSign ); 2104 } 2105 2106 } 2107 2108 /* 2109 ------------------------------------------------------------------------------- 2110 Returns the result of multiplying the double-precision floating-point values 2111 `a' and `b'. The operation is performed according to the IEC/IEEE Standard 2112 for Binary Floating-point Arithmetic. 2113 ------------------------------------------------------------------------------- 2114 */ 2115 float64 float64_mul( struct roundingData *roundData, float64 a, float64 b ) 2116 { 2117 flag aSign, bSign, zSign; 2118 int16 aExp, bExp, zExp; 2119 bits64 aSig, bSig, zSig0, zSig1; 2120 2121 aSig = extractFloat64Frac( a ); 2122 aExp = extractFloat64Exp( a ); 2123 aSign = extractFloat64Sign( a ); 2124 bSig = extractFloat64Frac( b ); 2125 bExp = extractFloat64Exp( b ); 2126 bSign = extractFloat64Sign( b ); 2127 zSign = aSign ^ bSign; 2128 if ( aExp == 0x7FF ) { 2129 if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) { 2130 return propagateFloat64NaN( a, b ); 2131 } 2132 if ( ( bExp | bSig ) == 0 ) { 2133 roundData->exception |= float_flag_invalid; 2134 return float64_default_nan; 2135 } 2136 return packFloat64( zSign, 0x7FF, 0 ); 2137 } 2138 if ( bExp == 0x7FF ) { 2139 if ( bSig ) return propagateFloat64NaN( a, b ); 2140 if ( ( aExp | aSig ) == 0 ) { 2141 roundData->exception |= float_flag_invalid; 2142 return float64_default_nan; 2143 } 2144 return packFloat64( zSign, 0x7FF, 0 ); 2145 } 2146 if ( aExp == 0 ) { 2147 if ( aSig == 0 ) return packFloat64( zSign, 0, 0 ); 2148 normalizeFloat64Subnormal( aSig, &aExp, &aSig ); 2149 } 2150 if ( bExp == 0 ) { 2151 if ( bSig == 0 ) return packFloat64( zSign, 0, 0 ); 2152 normalizeFloat64Subnormal( bSig, &bExp, &bSig ); 2153 } 2154 zExp = aExp + bExp - 0x3FF; 2155 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10; 2156 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11; 2157 mul64To128( aSig, bSig, &zSig0, &zSig1 ); 2158 zSig0 |= ( zSig1 != 0 ); 2159 if ( 0 <= (sbits64) ( zSig0<<1 ) ) { 2160 zSig0 <<= 1; 2161 --zExp; 2162 } 2163 return roundAndPackFloat64( roundData, zSign, zExp, zSig0 ); 2164 2165 } 2166 2167 /* 2168 ------------------------------------------------------------------------------- 2169 Returns the result of dividing the double-precision floating-point value `a' 2170 by the corresponding value `b'. The operation is performed according to 2171 the IEC/IEEE Standard for Binary Floating-point Arithmetic. 2172 ------------------------------------------------------------------------------- 2173 */ 2174 float64 float64_div( struct roundingData *roundData, float64 a, float64 b ) 2175 { 2176 flag aSign, bSign, zSign; 2177 int16 aExp, bExp, zExp; 2178 bits64 aSig, bSig, zSig; 2179 bits64 rem0, rem1; 2180 bits64 term0, term1; 2181 2182 aSig = extractFloat64Frac( a ); 2183 aExp = extractFloat64Exp( a ); 2184 aSign = extractFloat64Sign( a ); 2185 bSig = extractFloat64Frac( b ); 2186 bExp = extractFloat64Exp( b ); 2187 bSign = extractFloat64Sign( b ); 2188 zSign = aSign ^ bSign; 2189 if ( aExp == 0x7FF ) { 2190 if ( aSig ) return propagateFloat64NaN( a, b ); 2191 if ( bExp == 0x7FF ) { 2192 if ( bSig ) return propagateFloat64NaN( a, b ); 2193 roundData->exception |= float_flag_invalid; 2194 return float64_default_nan; 2195 } 2196 return packFloat64( zSign, 0x7FF, 0 ); 2197 } 2198 if ( bExp == 0x7FF ) { 2199 if ( bSig ) return propagateFloat64NaN( a, b ); 2200 return packFloat64( zSign, 0, 0 ); 2201 } 2202 if ( bExp == 0 ) { 2203 if ( bSig == 0 ) { 2204 if ( ( aExp | aSig ) == 0 ) { 2205 roundData->exception |= float_flag_invalid; 2206 return float64_default_nan; 2207 } 2208 roundData->exception |= float_flag_divbyzero; 2209 return packFloat64( zSign, 0x7FF, 0 ); 2210 } 2211 normalizeFloat64Subnormal( bSig, &bExp, &bSig ); 2212 } 2213 if ( aExp == 0 ) { 2214 if ( aSig == 0 ) return packFloat64( zSign, 0, 0 ); 2215 normalizeFloat64Subnormal( aSig, &aExp, &aSig ); 2216 } 2217 zExp = aExp - bExp + 0x3FD; 2218 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10; 2219 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11; 2220 if ( bSig <= ( aSig + aSig ) ) { 2221 aSig >>= 1; 2222 ++zExp; 2223 } 2224 zSig = estimateDiv128To64( aSig, 0, bSig ); 2225 if ( ( zSig & 0x1FF ) <= 2 ) { 2226 mul64To128( bSig, zSig, &term0, &term1 ); 2227 sub128( aSig, 0, term0, term1, &rem0, &rem1 ); 2228 while ( (sbits64) rem0 < 0 ) { 2229 --zSig; 2230 add128( rem0, rem1, 0, bSig, &rem0, &rem1 ); 2231 } 2232 zSig |= ( rem1 != 0 ); 2233 } 2234 return roundAndPackFloat64( roundData, zSign, zExp, zSig ); 2235 2236 } 2237 2238 /* 2239 ------------------------------------------------------------------------------- 2240 Returns the remainder of the double-precision floating-point value `a' 2241 with respect to the corresponding value `b'. The operation is performed 2242 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 2243 ------------------------------------------------------------------------------- 2244 */ 2245 float64 float64_rem( struct roundingData *roundData, float64 a, float64 b ) 2246 { 2247 flag aSign, bSign, zSign; 2248 int16 aExp, bExp, expDiff; 2249 bits64 aSig, bSig; 2250 bits64 q, alternateASig; 2251 sbits64 sigMean; 2252 2253 aSig = extractFloat64Frac( a ); 2254 aExp = extractFloat64Exp( a ); 2255 aSign = extractFloat64Sign( a ); 2256 bSig = extractFloat64Frac( b ); 2257 bExp = extractFloat64Exp( b ); 2258 bSign = extractFloat64Sign( b ); 2259 if ( aExp == 0x7FF ) { 2260 if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) { 2261 return propagateFloat64NaN( a, b ); 2262 } 2263 roundData->exception |= float_flag_invalid; 2264 return float64_default_nan; 2265 } 2266 if ( bExp == 0x7FF ) { 2267 if ( bSig ) return propagateFloat64NaN( a, b ); 2268 return a; 2269 } 2270 if ( bExp == 0 ) { 2271 if ( bSig == 0 ) { 2272 roundData->exception |= float_flag_invalid; 2273 return float64_default_nan; 2274 } 2275 normalizeFloat64Subnormal( bSig, &bExp, &bSig ); 2276 } 2277 if ( aExp == 0 ) { 2278 if ( aSig == 0 ) return a; 2279 normalizeFloat64Subnormal( aSig, &aExp, &aSig ); 2280 } 2281 expDiff = aExp - bExp; 2282 aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11; 2283 bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11; 2284 if ( expDiff < 0 ) { 2285 if ( expDiff < -1 ) return a; 2286 aSig >>= 1; 2287 } 2288 q = ( bSig <= aSig ); 2289 if ( q ) aSig -= bSig; 2290 expDiff -= 64; 2291 while ( 0 < expDiff ) { 2292 q = estimateDiv128To64( aSig, 0, bSig ); 2293 q = ( 2 < q ) ? q - 2 : 0; 2294 aSig = - ( ( bSig>>2 ) * q ); 2295 expDiff -= 62; 2296 } 2297 expDiff += 64; 2298 if ( 0 < expDiff ) { 2299 q = estimateDiv128To64( aSig, 0, bSig ); 2300 q = ( 2 < q ) ? q - 2 : 0; 2301 q >>= 64 - expDiff; 2302 bSig >>= 2; 2303 aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q; 2304 } 2305 else { 2306 aSig >>= 2; 2307 bSig >>= 2; 2308 } 2309 do { 2310 alternateASig = aSig; 2311 ++q; 2312 aSig -= bSig; 2313 } while ( 0 <= (sbits64) aSig ); 2314 sigMean = aSig + alternateASig; 2315 if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) { 2316 aSig = alternateASig; 2317 } 2318 zSign = ( (sbits64) aSig < 0 ); 2319 if ( zSign ) aSig = - aSig; 2320 return normalizeRoundAndPackFloat64( roundData, aSign ^ zSign, bExp, aSig ); 2321 2322 } 2323 2324 /* 2325 ------------------------------------------------------------------------------- 2326 Returns the square root of the double-precision floating-point value `a'. 2327 The operation is performed according to the IEC/IEEE Standard for Binary 2328 Floating-point Arithmetic. 2329 ------------------------------------------------------------------------------- 2330 */ 2331 float64 float64_sqrt( struct roundingData *roundData, float64 a ) 2332 { 2333 flag aSign; 2334 int16 aExp, zExp; 2335 bits64 aSig, zSig; 2336 bits64 rem0, rem1, term0, term1; //, shiftedRem; 2337 //float64 z; 2338 2339 aSig = extractFloat64Frac( a ); 2340 aExp = extractFloat64Exp( a ); 2341 aSign = extractFloat64Sign( a ); 2342 if ( aExp == 0x7FF ) { 2343 if ( aSig ) return propagateFloat64NaN( a, a ); 2344 if ( ! aSign ) return a; 2345 roundData->exception |= float_flag_invalid; 2346 return float64_default_nan; 2347 } 2348 if ( aSign ) { 2349 if ( ( aExp | aSig ) == 0 ) return a; 2350 roundData->exception |= float_flag_invalid; 2351 return float64_default_nan; 2352 } 2353 if ( aExp == 0 ) { 2354 if ( aSig == 0 ) return 0; 2355 normalizeFloat64Subnormal( aSig, &aExp, &aSig ); 2356 } 2357 zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE; 2358 aSig |= LIT64( 0x0010000000000000 ); 2359 zSig = estimateSqrt32( aExp, aSig>>21 ); 2360 zSig <<= 31; 2361 aSig <<= 9 - ( aExp & 1 ); 2362 zSig = estimateDiv128To64( aSig, 0, zSig ) + zSig + 2; 2363 if ( ( zSig & 0x3FF ) <= 5 ) { 2364 if ( zSig < 2 ) { 2365 zSig = LIT64( 0xFFFFFFFFFFFFFFFF ); 2366 } 2367 else { 2368 aSig <<= 2; 2369 mul64To128( zSig, zSig, &term0, &term1 ); 2370 sub128( aSig, 0, term0, term1, &rem0, &rem1 ); 2371 while ( (sbits64) rem0 < 0 ) { 2372 --zSig; 2373 shortShift128Left( 0, zSig, 1, &term0, &term1 ); 2374 term1 |= 1; 2375 add128( rem0, rem1, term0, term1, &rem0, &rem1 ); 2376 } 2377 zSig |= ( ( rem0 | rem1 ) != 0 ); 2378 } 2379 } 2380 shift64RightJamming( zSig, 1, &zSig ); 2381 return roundAndPackFloat64( roundData, 0, zExp, zSig ); 2382 2383 } 2384 2385 /* 2386 ------------------------------------------------------------------------------- 2387 Returns 1 if the double-precision floating-point value `a' is equal to the 2388 corresponding value `b', and 0 otherwise. The comparison is performed 2389 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 2390 ------------------------------------------------------------------------------- 2391 */ 2392 flag float64_eq( float64 a, float64 b ) 2393 { 2394 2395 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) ) 2396 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) ) 2397 ) { 2398 if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) { 2399 float_raise( float_flag_invalid ); 2400 } 2401 return 0; 2402 } 2403 return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 ); 2404 2405 } 2406 2407 /* 2408 ------------------------------------------------------------------------------- 2409 Returns 1 if the double-precision floating-point value `a' is less than or 2410 equal to the corresponding value `b', and 0 otherwise. The comparison is 2411 performed according to the IEC/IEEE Standard for Binary Floating-point 2412 Arithmetic. 2413 ------------------------------------------------------------------------------- 2414 */ 2415 flag float64_le( float64 a, float64 b ) 2416 { 2417 flag aSign, bSign; 2418 2419 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) ) 2420 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) ) 2421 ) { 2422 float_raise( float_flag_invalid ); 2423 return 0; 2424 } 2425 aSign = extractFloat64Sign( a ); 2426 bSign = extractFloat64Sign( b ); 2427 if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 ); 2428 return ( a == b ) || ( aSign ^ ( a < b ) ); 2429 2430 } 2431 2432 /* 2433 ------------------------------------------------------------------------------- 2434 Returns 1 if the double-precision floating-point value `a' is less than 2435 the corresponding value `b', and 0 otherwise. The comparison is performed 2436 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 2437 ------------------------------------------------------------------------------- 2438 */ 2439 flag float64_lt( float64 a, float64 b ) 2440 { 2441 flag aSign, bSign; 2442 2443 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) ) 2444 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) ) 2445 ) { 2446 float_raise( float_flag_invalid ); 2447 return 0; 2448 } 2449 aSign = extractFloat64Sign( a ); 2450 bSign = extractFloat64Sign( b ); 2451 if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 ); 2452 return ( a != b ) && ( aSign ^ ( a < b ) ); 2453 2454 } 2455 2456 /* 2457 ------------------------------------------------------------------------------- 2458 Returns 1 if the double-precision floating-point value `a' is equal to the 2459 corresponding value `b', and 0 otherwise. The invalid exception is raised 2460 if either operand is a NaN. Otherwise, the comparison is performed 2461 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 2462 ------------------------------------------------------------------------------- 2463 */ 2464 flag float64_eq_signaling( float64 a, float64 b ) 2465 { 2466 2467 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) ) 2468 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) ) 2469 ) { 2470 float_raise( float_flag_invalid ); 2471 return 0; 2472 } 2473 return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 ); 2474 2475 } 2476 2477 /* 2478 ------------------------------------------------------------------------------- 2479 Returns 1 if the double-precision floating-point value `a' is less than or 2480 equal to the corresponding value `b', and 0 otherwise. Quiet NaNs do not 2481 cause an exception. Otherwise, the comparison is performed according to the 2482 IEC/IEEE Standard for Binary Floating-point Arithmetic. 2483 ------------------------------------------------------------------------------- 2484 */ 2485 flag float64_le_quiet( float64 a, float64 b ) 2486 { 2487 flag aSign, bSign; 2488 //int16 aExp, bExp; 2489 2490 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) ) 2491 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) ) 2492 ) { 2493 /* Do nothing, even if NaN as we're quiet */ 2494 return 0; 2495 } 2496 aSign = extractFloat64Sign( a ); 2497 bSign = extractFloat64Sign( b ); 2498 if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 ); 2499 return ( a == b ) || ( aSign ^ ( a < b ) ); 2500 2501 } 2502 2503 /* 2504 ------------------------------------------------------------------------------- 2505 Returns 1 if the double-precision floating-point value `a' is less than 2506 the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause an 2507 exception. Otherwise, the comparison is performed according to the IEC/IEEE 2508 Standard for Binary Floating-point Arithmetic. 2509 ------------------------------------------------------------------------------- 2510 */ 2511 flag float64_lt_quiet( float64 a, float64 b ) 2512 { 2513 flag aSign, bSign; 2514 2515 if ( ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) ) 2516 || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) ) 2517 ) { 2518 /* Do nothing, even if NaN as we're quiet */ 2519 return 0; 2520 } 2521 aSign = extractFloat64Sign( a ); 2522 bSign = extractFloat64Sign( b ); 2523 if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 ); 2524 return ( a != b ) && ( aSign ^ ( a < b ) ); 2525 2526 } 2527 2528 #ifdef FLOATX80 2529 2530 /* 2531 ------------------------------------------------------------------------------- 2532 Returns the result of converting the extended double-precision floating- 2533 point value `a' to the 32-bit two's complement integer format. The 2534 conversion is performed according to the IEC/IEEE Standard for Binary 2535 Floating-point Arithmetic---which means in particular that the conversion 2536 is rounded according to the current rounding mode. If `a' is a NaN, the 2537 largest positive integer is returned. Otherwise, if the conversion 2538 overflows, the largest integer with the same sign as `a' is returned. 2539 ------------------------------------------------------------------------------- 2540 */ 2541 int32 floatx80_to_int32( struct roundingData *roundData, floatx80 a ) 2542 { 2543 flag aSign; 2544 int32 aExp, shiftCount; 2545 bits64 aSig; 2546 2547 aSig = extractFloatx80Frac( a ); 2548 aExp = extractFloatx80Exp( a ); 2549 aSign = extractFloatx80Sign( a ); 2550 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0; 2551 shiftCount = 0x4037 - aExp; 2552 if ( shiftCount <= 0 ) shiftCount = 1; 2553 shift64RightJamming( aSig, shiftCount, &aSig ); 2554 return roundAndPackInt32( roundData, aSign, aSig ); 2555 2556 } 2557 2558 /* 2559 ------------------------------------------------------------------------------- 2560 Returns the result of converting the extended double-precision floating- 2561 point value `a' to the 32-bit two's complement integer format. The 2562 conversion is performed according to the IEC/IEEE Standard for Binary 2563 Floating-point Arithmetic, except that the conversion is always rounded 2564 toward zero. If `a' is a NaN, the largest positive integer is returned. 2565 Otherwise, if the conversion overflows, the largest integer with the same 2566 sign as `a' is returned. 2567 ------------------------------------------------------------------------------- 2568 */ 2569 int32 floatx80_to_int32_round_to_zero( floatx80 a ) 2570 { 2571 flag aSign; 2572 int32 aExp, shiftCount; 2573 bits64 aSig, savedASig; 2574 int32 z; 2575 2576 aSig = extractFloatx80Frac( a ); 2577 aExp = extractFloatx80Exp( a ); 2578 aSign = extractFloatx80Sign( a ); 2579 shiftCount = 0x403E - aExp; 2580 if ( shiftCount < 32 ) { 2581 if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0; 2582 goto invalid; 2583 } 2584 else if ( 63 < shiftCount ) { 2585 if ( aExp || aSig ) float_raise( float_flag_inexact ); 2586 return 0; 2587 } 2588 savedASig = aSig; 2589 aSig >>= shiftCount; 2590 z = aSig; 2591 if ( aSign ) z = - z; 2592 if ( ( z < 0 ) ^ aSign ) { 2593 invalid: 2594 float_raise( float_flag_invalid ); 2595 return aSign ? 0x80000000 : 0x7FFFFFFF; 2596 } 2597 if ( ( aSig<<shiftCount ) != savedASig ) { 2598 float_raise( float_flag_inexact ); 2599 } 2600 return z; 2601 2602 } 2603 2604 /* 2605 ------------------------------------------------------------------------------- 2606 Returns the result of converting the extended double-precision floating- 2607 point value `a' to the single-precision floating-point format. The 2608 conversion is performed according to the IEC/IEEE Standard for Binary 2609 Floating-point Arithmetic. 2610 ------------------------------------------------------------------------------- 2611 */ 2612 float32 floatx80_to_float32( struct roundingData *roundData, floatx80 a ) 2613 { 2614 flag aSign; 2615 int32 aExp; 2616 bits64 aSig; 2617 2618 aSig = extractFloatx80Frac( a ); 2619 aExp = extractFloatx80Exp( a ); 2620 aSign = extractFloatx80Sign( a ); 2621 if ( aExp == 0x7FFF ) { 2622 if ( (bits64) ( aSig<<1 ) ) { 2623 return commonNaNToFloat32( floatx80ToCommonNaN( a ) ); 2624 } 2625 return packFloat32( aSign, 0xFF, 0 ); 2626 } 2627 shift64RightJamming( aSig, 33, &aSig ); 2628 if ( aExp || aSig ) aExp -= 0x3F81; 2629 return roundAndPackFloat32( roundData, aSign, aExp, aSig ); 2630 2631 } 2632 2633 /* 2634 ------------------------------------------------------------------------------- 2635 Returns the result of converting the extended double-precision floating- 2636 point value `a' to the double-precision floating-point format. The 2637 conversion is performed according to the IEC/IEEE Standard for Binary 2638 Floating-point Arithmetic. 2639 ------------------------------------------------------------------------------- 2640 */ 2641 float64 floatx80_to_float64( struct roundingData *roundData, floatx80 a ) 2642 { 2643 flag aSign; 2644 int32 aExp; 2645 bits64 aSig, zSig; 2646 2647 aSig = extractFloatx80Frac( a ); 2648 aExp = extractFloatx80Exp( a ); 2649 aSign = extractFloatx80Sign( a ); 2650 if ( aExp == 0x7FFF ) { 2651 if ( (bits64) ( aSig<<1 ) ) { 2652 return commonNaNToFloat64( floatx80ToCommonNaN( a ) ); 2653 } 2654 return packFloat64( aSign, 0x7FF, 0 ); 2655 } 2656 shift64RightJamming( aSig, 1, &zSig ); 2657 if ( aExp || aSig ) aExp -= 0x3C01; 2658 return roundAndPackFloat64( roundData, aSign, aExp, zSig ); 2659 2660 } 2661 2662 /* 2663 ------------------------------------------------------------------------------- 2664 Rounds the extended double-precision floating-point value `a' to an integer, 2665 and returns the result as an extended quadruple-precision floating-point 2666 value. The operation is performed according to the IEC/IEEE Standard for 2667 Binary Floating-point Arithmetic. 2668 ------------------------------------------------------------------------------- 2669 */ 2670 floatx80 floatx80_round_to_int( struct roundingData *roundData, floatx80 a ) 2671 { 2672 flag aSign; 2673 int32 aExp; 2674 bits64 lastBitMask, roundBitsMask; 2675 int8 roundingMode; 2676 floatx80 z; 2677 2678 aExp = extractFloatx80Exp( a ); 2679 if ( 0x403E <= aExp ) { 2680 if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) { 2681 return propagateFloatx80NaN( a, a ); 2682 } 2683 return a; 2684 } 2685 if ( aExp <= 0x3FFE ) { 2686 if ( ( aExp == 0 ) 2687 && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) { 2688 return a; 2689 } 2690 roundData->exception |= float_flag_inexact; 2691 aSign = extractFloatx80Sign( a ); 2692 switch ( roundData->mode ) { 2693 case float_round_nearest_even: 2694 if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 ) 2695 ) { 2696 return 2697 packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) ); 2698 } 2699 break; 2700 case float_round_down: 2701 return 2702 aSign ? 2703 packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) ) 2704 : packFloatx80( 0, 0, 0 ); 2705 case float_round_up: 2706 return 2707 aSign ? packFloatx80( 1, 0, 0 ) 2708 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) ); 2709 } 2710 return packFloatx80( aSign, 0, 0 ); 2711 } 2712 lastBitMask = 1; 2713 lastBitMask <<= 0x403E - aExp; 2714 roundBitsMask = lastBitMask - 1; 2715 z = a; 2716 roundingMode = roundData->mode; 2717 if ( roundingMode == float_round_nearest_even ) { 2718 z.low += lastBitMask>>1; 2719 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask; 2720 } 2721 else if ( roundingMode != float_round_to_zero ) { 2722 if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) { 2723 z.low += roundBitsMask; 2724 } 2725 } 2726 z.low &= ~ roundBitsMask; 2727 if ( z.low == 0 ) { 2728 ++z.high; 2729 z.low = LIT64( 0x8000000000000000 ); 2730 } 2731 if ( z.low != a.low ) roundData->exception |= float_flag_inexact; 2732 return z; 2733 2734 } 2735 2736 /* 2737 ------------------------------------------------------------------------------- 2738 Returns the result of adding the absolute values of the extended double- 2739 precision floating-point values `a' and `b'. If `zSign' is true, the sum is 2740 negated before being returned. `zSign' is ignored if the result is a NaN. 2741 The addition is performed according to the IEC/IEEE Standard for Binary 2742 Floating-point Arithmetic. 2743 ------------------------------------------------------------------------------- 2744 */ 2745 static floatx80 addFloatx80Sigs( struct roundingData *roundData, floatx80 a, floatx80 b, flag zSign ) 2746 { 2747 int32 aExp, bExp, zExp; 2748 bits64 aSig, bSig, zSig0, zSig1; 2749 int32 expDiff; 2750 2751 aSig = extractFloatx80Frac( a ); 2752 aExp = extractFloatx80Exp( a ); 2753 bSig = extractFloatx80Frac( b ); 2754 bExp = extractFloatx80Exp( b ); 2755 expDiff = aExp - bExp; 2756 if ( 0 < expDiff ) { 2757 if ( aExp == 0x7FFF ) { 2758 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b ); 2759 return a; 2760 } 2761 if ( bExp == 0 ) --expDiff; 2762 shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 ); 2763 zExp = aExp; 2764 } 2765 else if ( expDiff < 0 ) { 2766 if ( bExp == 0x7FFF ) { 2767 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b ); 2768 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 2769 } 2770 if ( aExp == 0 ) ++expDiff; 2771 shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 ); 2772 zExp = bExp; 2773 } 2774 else { 2775 if ( aExp == 0x7FFF ) { 2776 if ( (bits64) ( ( aSig | bSig )<<1 ) ) { 2777 return propagateFloatx80NaN( a, b ); 2778 } 2779 return a; 2780 } 2781 zSig1 = 0; 2782 zSig0 = aSig + bSig; 2783 if ( aExp == 0 ) { 2784 normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 ); 2785 goto roundAndPack; 2786 } 2787 zExp = aExp; 2788 goto shiftRight1; 2789 } 2790 2791 zSig0 = aSig + bSig; 2792 2793 if ( (sbits64) zSig0 < 0 ) goto roundAndPack; 2794 shiftRight1: 2795 shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 ); 2796 zSig0 |= LIT64( 0x8000000000000000 ); 2797 ++zExp; 2798 roundAndPack: 2799 return 2800 roundAndPackFloatx80( 2801 roundData, zSign, zExp, zSig0, zSig1 ); 2802 2803 } 2804 2805 /* 2806 ------------------------------------------------------------------------------- 2807 Returns the result of subtracting the absolute values of the extended 2808 double-precision floating-point values `a' and `b'. If `zSign' is true, 2809 the difference is negated before being returned. `zSign' is ignored if the 2810 result is a NaN. The subtraction is performed according to the IEC/IEEE 2811 Standard for Binary Floating-point Arithmetic. 2812 ------------------------------------------------------------------------------- 2813 */ 2814 static floatx80 subFloatx80Sigs( struct roundingData *roundData, floatx80 a, floatx80 b, flag zSign ) 2815 { 2816 int32 aExp, bExp, zExp; 2817 bits64 aSig, bSig, zSig0, zSig1; 2818 int32 expDiff; 2819 floatx80 z; 2820 2821 aSig = extractFloatx80Frac( a ); 2822 aExp = extractFloatx80Exp( a ); 2823 bSig = extractFloatx80Frac( b ); 2824 bExp = extractFloatx80Exp( b ); 2825 expDiff = aExp - bExp; 2826 if ( 0 < expDiff ) goto aExpBigger; 2827 if ( expDiff < 0 ) goto bExpBigger; 2828 if ( aExp == 0x7FFF ) { 2829 if ( (bits64) ( ( aSig | bSig )<<1 ) ) { 2830 return propagateFloatx80NaN( a, b ); 2831 } 2832 roundData->exception |= float_flag_invalid; 2833 z.low = floatx80_default_nan_low; 2834 z.high = floatx80_default_nan_high; 2835 z.__padding = 0; 2836 return z; 2837 } 2838 if ( aExp == 0 ) { 2839 aExp = 1; 2840 bExp = 1; 2841 } 2842 zSig1 = 0; 2843 if ( bSig < aSig ) goto aBigger; 2844 if ( aSig < bSig ) goto bBigger; 2845 return packFloatx80( roundData->mode == float_round_down, 0, 0 ); 2846 bExpBigger: 2847 if ( bExp == 0x7FFF ) { 2848 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b ); 2849 return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) ); 2850 } 2851 if ( aExp == 0 ) ++expDiff; 2852 shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 ); 2853 bBigger: 2854 sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 ); 2855 zExp = bExp; 2856 zSign ^= 1; 2857 goto normalizeRoundAndPack; 2858 aExpBigger: 2859 if ( aExp == 0x7FFF ) { 2860 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b ); 2861 return a; 2862 } 2863 if ( bExp == 0 ) --expDiff; 2864 shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 ); 2865 aBigger: 2866 sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 ); 2867 zExp = aExp; 2868 normalizeRoundAndPack: 2869 return 2870 normalizeRoundAndPackFloatx80( 2871 roundData, zSign, zExp, zSig0, zSig1 ); 2872 2873 } 2874 2875 /* 2876 ------------------------------------------------------------------------------- 2877 Returns the result of adding the extended double-precision floating-point 2878 values `a' and `b'. The operation is performed according to the IEC/IEEE 2879 Standard for Binary Floating-point Arithmetic. 2880 ------------------------------------------------------------------------------- 2881 */ 2882 floatx80 floatx80_add( struct roundingData *roundData, floatx80 a, floatx80 b ) 2883 { 2884 flag aSign, bSign; 2885 2886 aSign = extractFloatx80Sign( a ); 2887 bSign = extractFloatx80Sign( b ); 2888 if ( aSign == bSign ) { 2889 return addFloatx80Sigs( roundData, a, b, aSign ); 2890 } 2891 else { 2892 return subFloatx80Sigs( roundData, a, b, aSign ); 2893 } 2894 2895 } 2896 2897 /* 2898 ------------------------------------------------------------------------------- 2899 Returns the result of subtracting the extended double-precision floating- 2900 point values `a' and `b'. The operation is performed according to the 2901 IEC/IEEE Standard for Binary Floating-point Arithmetic. 2902 ------------------------------------------------------------------------------- 2903 */ 2904 floatx80 floatx80_sub( struct roundingData *roundData, floatx80 a, floatx80 b ) 2905 { 2906 flag aSign, bSign; 2907 2908 aSign = extractFloatx80Sign( a ); 2909 bSign = extractFloatx80Sign( b ); 2910 if ( aSign == bSign ) { 2911 return subFloatx80Sigs( roundData, a, b, aSign ); 2912 } 2913 else { 2914 return addFloatx80Sigs( roundData, a, b, aSign ); 2915 } 2916 2917 } 2918 2919 /* 2920 ------------------------------------------------------------------------------- 2921 Returns the result of multiplying the extended double-precision floating- 2922 point values `a' and `b'. The operation is performed according to the 2923 IEC/IEEE Standard for Binary Floating-point Arithmetic. 2924 ------------------------------------------------------------------------------- 2925 */ 2926 floatx80 floatx80_mul( struct roundingData *roundData, floatx80 a, floatx80 b ) 2927 { 2928 flag aSign, bSign, zSign; 2929 int32 aExp, bExp, zExp; 2930 bits64 aSig, bSig, zSig0, zSig1; 2931 floatx80 z; 2932 2933 aSig = extractFloatx80Frac( a ); 2934 aExp = extractFloatx80Exp( a ); 2935 aSign = extractFloatx80Sign( a ); 2936 bSig = extractFloatx80Frac( b ); 2937 bExp = extractFloatx80Exp( b ); 2938 bSign = extractFloatx80Sign( b ); 2939 zSign = aSign ^ bSign; 2940 if ( aExp == 0x7FFF ) { 2941 if ( (bits64) ( aSig<<1 ) 2942 || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) { 2943 return propagateFloatx80NaN( a, b ); 2944 } 2945 if ( ( bExp | bSig ) == 0 ) goto invalid; 2946 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 2947 } 2948 if ( bExp == 0x7FFF ) { 2949 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b ); 2950 if ( ( aExp | aSig ) == 0 ) { 2951 invalid: 2952 roundData->exception |= float_flag_invalid; 2953 z.low = floatx80_default_nan_low; 2954 z.high = floatx80_default_nan_high; 2955 z.__padding = 0; 2956 return z; 2957 } 2958 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 2959 } 2960 if ( aExp == 0 ) { 2961 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 ); 2962 normalizeFloatx80Subnormal( aSig, &aExp, &aSig ); 2963 } 2964 if ( bExp == 0 ) { 2965 if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 ); 2966 normalizeFloatx80Subnormal( bSig, &bExp, &bSig ); 2967 } 2968 zExp = aExp + bExp - 0x3FFE; 2969 mul64To128( aSig, bSig, &zSig0, &zSig1 ); 2970 if ( 0 < (sbits64) zSig0 ) { 2971 shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 ); 2972 --zExp; 2973 } 2974 return 2975 roundAndPackFloatx80( 2976 roundData, zSign, zExp, zSig0, zSig1 ); 2977 2978 } 2979 2980 /* 2981 ------------------------------------------------------------------------------- 2982 Returns the result of dividing the extended double-precision floating-point 2983 value `a' by the corresponding value `b'. The operation is performed 2984 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 2985 ------------------------------------------------------------------------------- 2986 */ 2987 floatx80 floatx80_div( struct roundingData *roundData, floatx80 a, floatx80 b ) 2988 { 2989 flag aSign, bSign, zSign; 2990 int32 aExp, bExp, zExp; 2991 bits64 aSig, bSig, zSig0, zSig1; 2992 bits64 rem0, rem1, rem2, term0, term1, term2; 2993 floatx80 z; 2994 2995 aSig = extractFloatx80Frac( a ); 2996 aExp = extractFloatx80Exp( a ); 2997 aSign = extractFloatx80Sign( a ); 2998 bSig = extractFloatx80Frac( b ); 2999 bExp = extractFloatx80Exp( b ); 3000 bSign = extractFloatx80Sign( b ); 3001 zSign = aSign ^ bSign; 3002 if ( aExp == 0x7FFF ) { 3003 if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b ); 3004 if ( bExp == 0x7FFF ) { 3005 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b ); 3006 goto invalid; 3007 } 3008 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 3009 } 3010 if ( bExp == 0x7FFF ) { 3011 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b ); 3012 return packFloatx80( zSign, 0, 0 ); 3013 } 3014 if ( bExp == 0 ) { 3015 if ( bSig == 0 ) { 3016 if ( ( aExp | aSig ) == 0 ) { 3017 invalid: 3018 roundData->exception |= float_flag_invalid; 3019 z.low = floatx80_default_nan_low; 3020 z.high = floatx80_default_nan_high; 3021 z.__padding = 0; 3022 return z; 3023 } 3024 roundData->exception |= float_flag_divbyzero; 3025 return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) ); 3026 } 3027 normalizeFloatx80Subnormal( bSig, &bExp, &bSig ); 3028 } 3029 if ( aExp == 0 ) { 3030 if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 ); 3031 normalizeFloatx80Subnormal( aSig, &aExp, &aSig ); 3032 } 3033 zExp = aExp - bExp + 0x3FFE; 3034 rem1 = 0; 3035 if ( bSig <= aSig ) { 3036 shift128Right( aSig, 0, 1, &aSig, &rem1 ); 3037 ++zExp; 3038 } 3039 zSig0 = estimateDiv128To64( aSig, rem1, bSig ); 3040 mul64To128( bSig, zSig0, &term0, &term1 ); 3041 sub128( aSig, rem1, term0, term1, &rem0, &rem1 ); 3042 while ( (sbits64) rem0 < 0 ) { 3043 --zSig0; 3044 add128( rem0, rem1, 0, bSig, &rem0, &rem1 ); 3045 } 3046 zSig1 = estimateDiv128To64( rem1, 0, bSig ); 3047 if ( (bits64) ( zSig1<<1 ) <= 8 ) { 3048 mul64To128( bSig, zSig1, &term1, &term2 ); 3049 sub128( rem1, 0, term1, term2, &rem1, &rem2 ); 3050 while ( (sbits64) rem1 < 0 ) { 3051 --zSig1; 3052 add128( rem1, rem2, 0, bSig, &rem1, &rem2 ); 3053 } 3054 zSig1 |= ( ( rem1 | rem2 ) != 0 ); 3055 } 3056 return 3057 roundAndPackFloatx80( 3058 roundData, zSign, zExp, zSig0, zSig1 ); 3059 3060 } 3061 3062 /* 3063 ------------------------------------------------------------------------------- 3064 Returns the remainder of the extended double-precision floating-point value 3065 `a' with respect to the corresponding value `b'. The operation is performed 3066 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 3067 ------------------------------------------------------------------------------- 3068 */ 3069 floatx80 floatx80_rem( struct roundingData *roundData, floatx80 a, floatx80 b ) 3070 { 3071 flag aSign, bSign, zSign; 3072 int32 aExp, bExp, expDiff; 3073 bits64 aSig0, aSig1, bSig; 3074 bits64 q, term0, term1, alternateASig0, alternateASig1; 3075 floatx80 z; 3076 3077 aSig0 = extractFloatx80Frac( a ); 3078 aExp = extractFloatx80Exp( a ); 3079 aSign = extractFloatx80Sign( a ); 3080 bSig = extractFloatx80Frac( b ); 3081 bExp = extractFloatx80Exp( b ); 3082 bSign = extractFloatx80Sign( b ); 3083 if ( aExp == 0x7FFF ) { 3084 if ( (bits64) ( aSig0<<1 ) 3085 || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) { 3086 return propagateFloatx80NaN( a, b ); 3087 } 3088 goto invalid; 3089 } 3090 if ( bExp == 0x7FFF ) { 3091 if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b ); 3092 return a; 3093 } 3094 if ( bExp == 0 ) { 3095 if ( bSig == 0 ) { 3096 invalid: 3097 roundData->exception |= float_flag_invalid; 3098 z.low = floatx80_default_nan_low; 3099 z.high = floatx80_default_nan_high; 3100 z.__padding = 0; 3101 return z; 3102 } 3103 normalizeFloatx80Subnormal( bSig, &bExp, &bSig ); 3104 } 3105 if ( aExp == 0 ) { 3106 if ( (bits64) ( aSig0<<1 ) == 0 ) return a; 3107 normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 ); 3108 } 3109 bSig |= LIT64( 0x8000000000000000 ); 3110 zSign = aSign; 3111 expDiff = aExp - bExp; 3112 aSig1 = 0; 3113 if ( expDiff < 0 ) { 3114 if ( expDiff < -1 ) return a; 3115 shift128Right( aSig0, 0, 1, &aSig0, &aSig1 ); 3116 expDiff = 0; 3117 } 3118 q = ( bSig <= aSig0 ); 3119 if ( q ) aSig0 -= bSig; 3120 expDiff -= 64; 3121 while ( 0 < expDiff ) { 3122 q = estimateDiv128To64( aSig0, aSig1, bSig ); 3123 q = ( 2 < q ) ? q - 2 : 0; 3124 mul64To128( bSig, q, &term0, &term1 ); 3125 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); 3126 shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 ); 3127 expDiff -= 62; 3128 } 3129 expDiff += 64; 3130 if ( 0 < expDiff ) { 3131 q = estimateDiv128To64( aSig0, aSig1, bSig ); 3132 q = ( 2 < q ) ? q - 2 : 0; 3133 q >>= 64 - expDiff; 3134 mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 ); 3135 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); 3136 shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 ); 3137 while ( le128( term0, term1, aSig0, aSig1 ) ) { 3138 ++q; 3139 sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); 3140 } 3141 } 3142 else { 3143 term1 = 0; 3144 term0 = bSig; 3145 } 3146 sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 ); 3147 if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 ) 3148 || ( eq128( alternateASig0, alternateASig1, aSig0, aSig1 ) 3149 && ( q & 1 ) ) 3150 ) { 3151 aSig0 = alternateASig0; 3152 aSig1 = alternateASig1; 3153 zSign = ! zSign; 3154 } 3155 3156 return 3157 normalizeRoundAndPackFloatx80( 3158 roundData, zSign, bExp + expDiff, aSig0, aSig1 ); 3159 3160 } 3161 3162 /* 3163 ------------------------------------------------------------------------------- 3164 Returns the square root of the extended double-precision floating-point 3165 value `a'. The operation is performed according to the IEC/IEEE Standard 3166 for Binary Floating-point Arithmetic. 3167 ------------------------------------------------------------------------------- 3168 */ 3169 floatx80 floatx80_sqrt( struct roundingData *roundData, floatx80 a ) 3170 { 3171 flag aSign; 3172 int32 aExp, zExp; 3173 bits64 aSig0, aSig1, zSig0, zSig1; 3174 bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3; 3175 bits64 shiftedRem0, shiftedRem1; 3176 floatx80 z; 3177 3178 aSig0 = extractFloatx80Frac( a ); 3179 aExp = extractFloatx80Exp( a ); 3180 aSign = extractFloatx80Sign( a ); 3181 if ( aExp == 0x7FFF ) { 3182 if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a ); 3183 if ( ! aSign ) return a; 3184 goto invalid; 3185 } 3186 if ( aSign ) { 3187 if ( ( aExp | aSig0 ) == 0 ) return a; 3188 invalid: 3189 roundData->exception |= float_flag_invalid; 3190 z.low = floatx80_default_nan_low; 3191 z.high = floatx80_default_nan_high; 3192 z.__padding = 0; 3193 return z; 3194 } 3195 if ( aExp == 0 ) { 3196 if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 ); 3197 normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 ); 3198 } 3199 zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF; 3200 zSig0 = estimateSqrt32( aExp, aSig0>>32 ); 3201 zSig0 <<= 31; 3202 aSig1 = 0; 3203 shift128Right( aSig0, 0, ( aExp & 1 ) + 2, &aSig0, &aSig1 ); 3204 zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0 ) + zSig0 + 4; 3205 if ( 0 <= (sbits64) zSig0 ) zSig0 = LIT64( 0xFFFFFFFFFFFFFFFF ); 3206 shortShift128Left( aSig0, aSig1, 2, &aSig0, &aSig1 ); 3207 mul64To128( zSig0, zSig0, &term0, &term1 ); 3208 sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 ); 3209 while ( (sbits64) rem0 < 0 ) { 3210 --zSig0; 3211 shortShift128Left( 0, zSig0, 1, &term0, &term1 ); 3212 term1 |= 1; 3213 add128( rem0, rem1, term0, term1, &rem0, &rem1 ); 3214 } 3215 shortShift128Left( rem0, rem1, 63, &shiftedRem0, &shiftedRem1 ); 3216 zSig1 = estimateDiv128To64( shiftedRem0, shiftedRem1, zSig0 ); 3217 if ( (bits64) ( zSig1<<1 ) <= 10 ) { 3218 if ( zSig1 == 0 ) zSig1 = 1; 3219 mul64To128( zSig0, zSig1, &term1, &term2 ); 3220 shortShift128Left( term1, term2, 1, &term1, &term2 ); 3221 sub128( rem1, 0, term1, term2, &rem1, &rem2 ); 3222 mul64To128( zSig1, zSig1, &term2, &term3 ); 3223 sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 ); 3224 while ( (sbits64) rem1 < 0 ) { 3225 --zSig1; 3226 shortShift192Left( 0, zSig0, zSig1, 1, &term1, &term2, &term3 ); 3227 term3 |= 1; 3228 add192( 3229 rem1, rem2, rem3, term1, term2, term3, &rem1, &rem2, &rem3 ); 3230 } 3231 zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 ); 3232 } 3233 return 3234 roundAndPackFloatx80( 3235 roundData, 0, zExp, zSig0, zSig1 ); 3236 3237 } 3238 3239 /* 3240 ------------------------------------------------------------------------------- 3241 Returns 1 if the extended double-precision floating-point value `a' is 3242 equal to the corresponding value `b', and 0 otherwise. The comparison is 3243 performed according to the IEC/IEEE Standard for Binary Floating-point 3244 Arithmetic. 3245 ------------------------------------------------------------------------------- 3246 */ 3247 flag floatx80_eq( floatx80 a, floatx80 b ) 3248 { 3249 3250 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF ) 3251 && (bits64) ( extractFloatx80Frac( a )<<1 ) ) 3252 || ( ( extractFloatx80Exp( b ) == 0x7FFF ) 3253 && (bits64) ( extractFloatx80Frac( b )<<1 ) ) 3254 ) { 3255 if ( floatx80_is_signaling_nan( a ) 3256 || floatx80_is_signaling_nan( b ) ) { 3257 float_raise( float_flag_invalid ); 3258 } 3259 return 0; 3260 } 3261 return 3262 ( a.low == b.low ) 3263 && ( ( a.high == b.high ) 3264 || ( ( a.low == 0 ) 3265 && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) ) 3266 ); 3267 3268 } 3269 3270 /* 3271 ------------------------------------------------------------------------------- 3272 Returns 1 if the extended double-precision floating-point value `a' is 3273 less than or equal to the corresponding value `b', and 0 otherwise. The 3274 comparison is performed according to the IEC/IEEE Standard for Binary 3275 Floating-point Arithmetic. 3276 ------------------------------------------------------------------------------- 3277 */ 3278 flag floatx80_le( floatx80 a, floatx80 b ) 3279 { 3280 flag aSign, bSign; 3281 3282 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF ) 3283 && (bits64) ( extractFloatx80Frac( a )<<1 ) ) 3284 || ( ( extractFloatx80Exp( b ) == 0x7FFF ) 3285 && (bits64) ( extractFloatx80Frac( b )<<1 ) ) 3286 ) { 3287 float_raise( float_flag_invalid ); 3288 return 0; 3289 } 3290 aSign = extractFloatx80Sign( a ); 3291 bSign = extractFloatx80Sign( b ); 3292 if ( aSign != bSign ) { 3293 return 3294 aSign 3295 || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low ) 3296 == 0 ); 3297 } 3298 return 3299 aSign ? le128( b.high, b.low, a.high, a.low ) 3300 : le128( a.high, a.low, b.high, b.low ); 3301 3302 } 3303 3304 /* 3305 ------------------------------------------------------------------------------- 3306 Returns 1 if the extended double-precision floating-point value `a' is 3307 less than the corresponding value `b', and 0 otherwise. The comparison 3308 is performed according to the IEC/IEEE Standard for Binary Floating-point 3309 Arithmetic. 3310 ------------------------------------------------------------------------------- 3311 */ 3312 flag floatx80_lt( floatx80 a, floatx80 b ) 3313 { 3314 flag aSign, bSign; 3315 3316 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF ) 3317 && (bits64) ( extractFloatx80Frac( a )<<1 ) ) 3318 || ( ( extractFloatx80Exp( b ) == 0x7FFF ) 3319 && (bits64) ( extractFloatx80Frac( b )<<1 ) ) 3320 ) { 3321 float_raise( float_flag_invalid ); 3322 return 0; 3323 } 3324 aSign = extractFloatx80Sign( a ); 3325 bSign = extractFloatx80Sign( b ); 3326 if ( aSign != bSign ) { 3327 return 3328 aSign 3329 && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low ) 3330 != 0 ); 3331 } 3332 return 3333 aSign ? lt128( b.high, b.low, a.high, a.low ) 3334 : lt128( a.high, a.low, b.high, b.low ); 3335 3336 } 3337 3338 /* 3339 ------------------------------------------------------------------------------- 3340 Returns 1 if the extended double-precision floating-point value `a' is equal 3341 to the corresponding value `b', and 0 otherwise. The invalid exception is 3342 raised if either operand is a NaN. Otherwise, the comparison is performed 3343 according to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 3344 ------------------------------------------------------------------------------- 3345 */ 3346 flag floatx80_eq_signaling( floatx80 a, floatx80 b ) 3347 { 3348 3349 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF ) 3350 && (bits64) ( extractFloatx80Frac( a )<<1 ) ) 3351 || ( ( extractFloatx80Exp( b ) == 0x7FFF ) 3352 && (bits64) ( extractFloatx80Frac( b )<<1 ) ) 3353 ) { 3354 float_raise( float_flag_invalid ); 3355 return 0; 3356 } 3357 return 3358 ( a.low == b.low ) 3359 && ( ( a.high == b.high ) 3360 || ( ( a.low == 0 ) 3361 && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) ) 3362 ); 3363 3364 } 3365 3366 /* 3367 ------------------------------------------------------------------------------- 3368 Returns 1 if the extended double-precision floating-point value `a' is less 3369 than or equal to the corresponding value `b', and 0 otherwise. Quiet NaNs 3370 do not cause an exception. Otherwise, the comparison is performed according 3371 to the IEC/IEEE Standard for Binary Floating-point Arithmetic. 3372 ------------------------------------------------------------------------------- 3373 */ 3374 flag floatx80_le_quiet( floatx80 a, floatx80 b ) 3375 { 3376 flag aSign, bSign; 3377 3378 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF ) 3379 && (bits64) ( extractFloatx80Frac( a )<<1 ) ) 3380 || ( ( extractFloatx80Exp( b ) == 0x7FFF ) 3381 && (bits64) ( extractFloatx80Frac( b )<<1 ) ) 3382 ) { 3383 /* Do nothing, even if NaN as we're quiet */ 3384 return 0; 3385 } 3386 aSign = extractFloatx80Sign( a ); 3387 bSign = extractFloatx80Sign( b ); 3388 if ( aSign != bSign ) { 3389 return 3390 aSign 3391 || ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low ) 3392 == 0 ); 3393 } 3394 return 3395 aSign ? le128( b.high, b.low, a.high, a.low ) 3396 : le128( a.high, a.low, b.high, b.low ); 3397 3398 } 3399 3400 /* 3401 ------------------------------------------------------------------------------- 3402 Returns 1 if the extended double-precision floating-point value `a' is less 3403 than the corresponding value `b', and 0 otherwise. Quiet NaNs do not cause 3404 an exception. Otherwise, the comparison is performed according to the 3405 IEC/IEEE Standard for Binary Floating-point Arithmetic. 3406 ------------------------------------------------------------------------------- 3407 */ 3408 flag floatx80_lt_quiet( floatx80 a, floatx80 b ) 3409 { 3410 flag aSign, bSign; 3411 3412 if ( ( ( extractFloatx80Exp( a ) == 0x7FFF ) 3413 && (bits64) ( extractFloatx80Frac( a )<<1 ) ) 3414 || ( ( extractFloatx80Exp( b ) == 0x7FFF ) 3415 && (bits64) ( extractFloatx80Frac( b )<<1 ) ) 3416 ) { 3417 /* Do nothing, even if NaN as we're quiet */ 3418 return 0; 3419 } 3420 aSign = extractFloatx80Sign( a ); 3421 bSign = extractFloatx80Sign( b ); 3422 if ( aSign != bSign ) { 3423 return 3424 aSign 3425 && ( ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low ) 3426 != 0 ); 3427 } 3428 return 3429 aSign ? lt128( b.high, b.low, a.high, a.low ) 3430 : lt128( a.high, a.low, b.high, b.low ); 3431 3432 } 3433 3434 #endif 3435 3436