vcs.maemo.org Git - qemu/blob - qemu/target-arm/nwfpe/softfloat.c

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