lib/libc/softfloat/bits32/softfloat.c

   1 /* $NetBSD: softfloat.c,v 1.3 2013/01/10 08:16:11 matt Exp $ */
   2
   3 /*
   4  * This version hacked for use with gcc -msoft-float by bjh21.
   5  * (Mostly a case of #ifdefing out things GCC doesn't need or provides
   6  *  itself).
   7  */
   8
   9 /*
  10  * Things you may want to define:
  11  *
  12  * SOFTFLOAT_FOR_GCC - build only those functions necessary for GCC (with
  13  *   -msoft-float) to work.  Include "softfloat-for-gcc.h" to get them
  14  *   properly renamed.
  15  */
  16
  17 /*
  18  * This differs from the standard bits32/softfloat.c in that float64
  19  * is defined to be a 64-bit integer rather than a structure.  The
  20  * structure is float64s, with translation between the two going via
  21  * float64u.
  22  */
  23
  24 /*
  25 ===============================================================================
  26
  27 This C source file is part of the SoftFloat IEC/IEEE Floating-Point
  28 Arithmetic Package, Release 2a.
  29
  30 Written by John R. Hauser.  This work was made possible in part by the
  31 International Computer Science Institute, located at Suite 600, 1947 Center
  32 Street, Berkeley, California 94704.  Funding was partially provided by the
  33 National Science Foundation under grant MIP-9311980.  The original version
  34 of this code was written as part of a project to build a fixed-point vector
  35 processor in collaboration with the University of California at Berkeley,
  36 overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
  37 is available through the Web page `http://HTTP.CS.Berkeley.EDU/~jhauser/
  38 arithmetic/SoftFloat.html'.
  39
  40 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
  41 has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
  42 TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
  43 PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
  44 AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
  45
  46 Derivative works are acceptable, even for commercial purposes, so long as
  47 (1) they include prominent notice that the work is derivative, and (2) they
  48 include prominent notice akin to these four paragraphs for those parts of
  49 this code that are retained.
  50
  51 ===============================================================================
  52 */
  53
  54 #include <sys/cdefs.h>
  55 #if defined(LIBC_SCCS) && !defined(lint)
  56 __RCSID("$NetBSD: softfloat.c,v 1.3 2013/01/10 08:16:11 matt Exp $");
  57 #endif /* LIBC_SCCS and not lint */
  58
  59 #ifdef SOFTFLOAT_FOR_GCC
  60 #include "softfloat-for-gcc.h"
  61 #endif
  62
  63 #include "milieu.h"
  64 #include "softfloat.h"
  65
  66 /*
  67  * Conversions between floats as stored in memory and floats as
  68  * SoftFloat uses them
  69  */
  70 #ifndef FLOAT64_DEMANGLE
  71 #define FLOAT64_DEMANGLE(a)     (a)
  72 #endif
  73 #ifndef FLOAT64_MANGLE
  74 #define FLOAT64_MANGLE(a)       (a)
  75 #endif
  76
  77 /*
  78 -------------------------------------------------------------------------------
  79 Floating-point rounding mode and exception flags.
  80 -------------------------------------------------------------------------------
  81 */
  82 #ifndef set_float_rounding_mode
  83 fp_rnd float_rounding_mode = float_round_nearest_even;
  84 fp_except float_exception_flags = 0;
  85 #endif
  86 #ifndef set_float_exception_inexact_flag
  87 #define set_float_exception_inexact_flag() \
  88         ((void)(float_exception_flags |= float_flag_inexact))
  89 #endif
  90
  91 /*
  92 -------------------------------------------------------------------------------
  93 Primitive arithmetic functions, including multi-word arithmetic, and
  94 division and square root approximations.  (Can be specialized to target if
  95 desired.)
  96 -------------------------------------------------------------------------------
  97 */
  98 #include "softfloat-macros"
  99
 100 /*
 101 -------------------------------------------------------------------------------
 102 Functions and definitions to determine:  (1) whether tininess for underflow
 103 is detected before or after rounding by default, (2) what (if anything)
 104 happens when exceptions are raised, (3) how signaling NaNs are distinguished
 105 from quiet NaNs, (4) the default generated quiet NaNs, and (4) how NaNs
 106 are propagated from function inputs to output.  These details are target-
 107 specific.
 108 -------------------------------------------------------------------------------
 109 */
 110 #include "softfloat-specialize"
 111
 112 /*
 113 -------------------------------------------------------------------------------
 114 Returns the fraction bits of the single-precision floating-point value `a'.
 115 -------------------------------------------------------------------------------
 116 */
 117 INLINE bits32 extractFloat32Frac( float32 a )
 118 {
 119
 120     return a & 0x007FFFFF;
 121
 122 }
 123
 124 /*
 125 -------------------------------------------------------------------------------
 126 Returns the exponent bits of the single-precision floating-point value `a'.
 127 -------------------------------------------------------------------------------
 128 */
 129 INLINE int16 extractFloat32Exp( float32 a )
 130 {
 131
 132     return ( a>>23 ) & 0xFF;
 133
 134 }
 135
 136 /*
 137 -------------------------------------------------------------------------------
 138 Returns the sign bit of the single-precision floating-point value `a'.
 139 -------------------------------------------------------------------------------
 140 */
 141 INLINE flag extractFloat32Sign( float32 a )
 142 {
 143
 144     return a>>31;
 145
 146 }
 147
 148 /*
 149 -------------------------------------------------------------------------------
 150 Normalizes the subnormal single-precision floating-point value represented
 151 by the denormalized significand `aSig'.  The normalized exponent and
 152 significand are stored at the locations pointed to by `zExpPtr' and
 153 `zSigPtr', respectively.
 154 -------------------------------------------------------------------------------
 155 */
 156 static void
 157  normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
 158 {
 159     int8 shiftCount;
 160
 161     shiftCount = countLeadingZeros32( aSig ) - 8;
 162     *zSigPtr = aSig<<shiftCount;
 163     *zExpPtr = 1 - shiftCount;
 164
 165 }
 166
 167 /*
 168 -------------------------------------------------------------------------------
 169 Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
 170 single-precision floating-point value, returning the result.  After being
 171 shifted into the proper positions, the three fields are simply added
 172 together to form the result.  This means that any integer portion of `zSig'
 173 will be added into the exponent.  Since a properly normalized significand
 174 will have an integer portion equal to 1, the `zExp' input should be 1 less
 175 than the desired result exponent whenever `zSig' is a complete, normalized
 176 significand.
 177 -------------------------------------------------------------------------------
 178 */
 179 INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
 180 {
 181
 182     return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
 183
 184 }
 185
 186 /*
 187 -------------------------------------------------------------------------------
 188 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 189 and significand `zSig', and returns the proper single-precision floating-
 190 point value corresponding to the abstract input.  Ordinarily, the abstract
 191 value is simply rounded and packed into the single-precision format, with
 192 the inexact exception raised if the abstract input cannot be represented
 193 exactly.  However, if the abstract value is too large, the overflow and
 194 inexact exceptions are raised and an infinity or maximal finite value is
 195 returned.  If the abstract value is too small, the input value is rounded to
 196 a subnormal number, and the underflow and inexact exceptions are raised if
 197 the abstract input cannot be represented exactly as a subnormal single-
 198 precision floating-point number.
 199     The input significand `zSig' has its binary point between bits 30
 200 and 29, which is 7 bits to the left of the usual location.  This shifted
 201 significand must be normalized or smaller.  If `zSig' is not normalized,
 202 `zExp' must be 0; in that case, the result returned is a subnormal number,
 203 and it must not require rounding.  In the usual case that `zSig' is
 204 normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
 205 The handling of underflow and overflow follows the IEC/IEEE Standard for
 206 Binary Floating-Point Arithmetic.
 207 -------------------------------------------------------------------------------
 208 */
 209 static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig )
 210 {
 211     int8 roundingMode;
 212     flag roundNearestEven;
 213     int8 roundIncrement, roundBits;
 214     flag isTiny;
 215
 216     roundingMode = float_rounding_mode;
 217     roundNearestEven = roundingMode == float_round_nearest_even;
 218     roundIncrement = 0x40;
 219     if ( ! roundNearestEven ) {
 220         if ( roundingMode == float_round_to_zero ) {
 221             roundIncrement = 0;
 222         }
 223         else {
 224             roundIncrement = 0x7F;
 225             if ( zSign ) {
 226                 if ( roundingMode == float_round_up ) roundIncrement = 0;
 227             }
 228             else {
 229                 if ( roundingMode == float_round_down ) roundIncrement = 0;
 230             }
 231         }
 232     }
 233     roundBits = zSig & 0x7F;
 234     if ( 0xFD <= (bits16) zExp ) {
 235         if (    ( 0xFD < zExp )
 236              || (    ( zExp == 0xFD )
 237                   && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
 238            ) {
 239             float_raise( float_flag_overflow | float_flag_inexact );
 240             return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
 241         }
 242         if ( zExp < 0 ) {
 243             isTiny =
 244                    ( float_detect_tininess == float_tininess_before_rounding )
 245                 || ( zExp < -1 )
 246                 || ( zSig + roundIncrement < (uint32)0x80000000 );
 247             shift32RightJamming( zSig, - zExp, &zSig );
 248             zExp = 0;
 249             roundBits = zSig & 0x7F;
 250             if ( isTiny && roundBits ) float_raise( float_flag_underflow );
 251         }
 252     }
 253     if ( roundBits ) set_float_exception_inexact_flag();
 254     zSig = ( zSig + roundIncrement )>>7;
 255     zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
 256     if ( zSig == 0 ) zExp = 0;
 257     return packFloat32( zSign, zExp, zSig );
 258
 259 }
 260
 261 /*
 262 -------------------------------------------------------------------------------
 263 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 264 and significand `zSig', and returns the proper single-precision floating-
 265 point value corresponding to the abstract input.  This routine is just like
 266 `roundAndPackFloat32' except that `zSig' does not have to be normalized.
 267 Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
 268 floating-point exponent.
 269 -------------------------------------------------------------------------------
 270 */
 271 static float32
 272  normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig )
 273 {
 274     int8 shiftCount;
 275
 276     shiftCount = countLeadingZeros32( zSig ) - 1;
 277     return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount );
 278
 279 }
 280
 281 /*
 282 -------------------------------------------------------------------------------
 283 Returns the least-significant 32 fraction bits of the double-precision
 284 floating-point value `a'.
 285 -------------------------------------------------------------------------------
 286 */
 287 INLINE bits32 extractFloat64Frac1( float64 a )
 288 {
 289
 290     return (bits32)(FLOAT64_DEMANGLE(a) & LIT64(0x00000000FFFFFFFF));
 291
 292 }
 293
 294 /*
 295 -------------------------------------------------------------------------------
 296 Returns the most-significant 20 fraction bits of the double-precision
 297 floating-point value `a'.
 298 -------------------------------------------------------------------------------
 299 */
 300 INLINE bits32 extractFloat64Frac0( float64 a )
 301 {
 302
 303     return (bits32)((FLOAT64_DEMANGLE(a) >> 32) & 0x000FFFFF);
 304
 305 }
 306
 307 /*
 308 -------------------------------------------------------------------------------
 309 Returns the exponent bits of the double-precision floating-point value `a'.
 310 -------------------------------------------------------------------------------
 311 */
 312 INLINE int16 extractFloat64Exp( float64 a )
 313 {
 314
 315     return (int16)((FLOAT64_DEMANGLE(a) >> 52) & 0x7FF);
 316
 317 }
 318
 319 /*
 320 -------------------------------------------------------------------------------
 321 Returns the sign bit of the double-precision floating-point value `a'.
 322 -------------------------------------------------------------------------------
 323 */
 324 INLINE flag extractFloat64Sign( float64 a )
 325 {
 326
 327     return (flag)(FLOAT64_DEMANGLE(a) >> 63);
 328
 329 }
 330
 331 /*
 332 -------------------------------------------------------------------------------
 333 Normalizes the subnormal double-precision floating-point value represented
 334 by the denormalized significand formed by the concatenation of `aSig0' and
 335 `aSig1'.  The normalized exponent is stored at the location pointed to by
 336 `zExpPtr'.  The most significant 21 bits of the normalized significand are
 337 stored at the location pointed to by `zSig0Ptr', and the least significant
 338 32 bits of the normalized significand are stored at the location pointed to
 339 by `zSig1Ptr'.
 340 -------------------------------------------------------------------------------
 341 */
 342 static void
 343  normalizeFloat64Subnormal(
 344      bits32 aSig0,
 345      bits32 aSig1,
 346      int16 *zExpPtr,
 347      bits32 *zSig0Ptr,
 348      bits32 *zSig1Ptr
 349  )
 350 {
 351     int8 shiftCount;
 352
 353     if ( aSig0 == 0 ) {
 354         shiftCount = countLeadingZeros32( aSig1 ) - 11;
 355         if ( shiftCount < 0 ) {
 356             *zSig0Ptr = aSig1>>( - shiftCount );
 357             *zSig1Ptr = aSig1<<( shiftCount & 31 );
 358         }
 359         else {
 360             *zSig0Ptr = aSig1<<shiftCount;
 361             *zSig1Ptr = 0;
 362         }
 363         *zExpPtr = - shiftCount - 31;
 364     }
 365     else {
 366         shiftCount = countLeadingZeros32( aSig0 ) - 11;
 367         shortShift64Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );
 368         *zExpPtr = 1 - shiftCount;
 369     }
 370
 371 }
 372
 373 /*
 374 -------------------------------------------------------------------------------
 375 Packs the sign `zSign', the exponent `zExp', and the significand formed by
 376 the concatenation of `zSig0' and `zSig1' into a double-precision floating-
 377 point value, returning the result.  After being shifted into the proper
 378 positions, the three fields `zSign', `zExp', and `zSig0' are simply added
 379 together to form the most significant 32 bits of the result.  This means
 380 that any integer portion of `zSig0' will be added into the exponent.  Since
 381 a properly normalized significand will have an integer portion equal to 1,
 382 the `zExp' input should be 1 less than the desired result exponent whenever
 383 `zSig0' and `zSig1' concatenated form a complete, normalized significand.
 384 -------------------------------------------------------------------------------
 385 */
 386 INLINE float64
 387  packFloat64( flag zSign, int16 zExp, bits32 zSig0, bits32 zSig1 )
 388 {
 389
 390     return FLOAT64_MANGLE( ( ( (bits64) zSign )<<63 ) +
 391                            ( ( (bits64) zExp )<<52 ) +
 392                            ( ( (bits64) zSig0 )<<32 ) + zSig1 );
 393
 394
 395 }
 396
 397 /*
 398 -------------------------------------------------------------------------------
 399 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 400 and extended significand formed by the concatenation of `zSig0', `zSig1',
 401 and `zSig2', and returns the proper double-precision floating-point value
 402 corresponding to the abstract input.  Ordinarily, the abstract value is
 403 simply rounded and packed into the double-precision format, with the inexact
 404 exception raised if the abstract input cannot be represented exactly.
 405 However, if the abstract value is too large, the overflow and inexact
 406 exceptions are raised and an infinity or maximal finite value is returned.
 407 If the abstract value is too small, the input value is rounded to a
 408 subnormal number, and the underflow and inexact exceptions are raised if the
 409 abstract input cannot be represented exactly as a subnormal double-precision
 410 floating-point number.
 411     The input significand must be normalized or smaller.  If the input
 412 significand is not normalized, `zExp' must be 0; in that case, the result
 413 returned is a subnormal number, and it must not require rounding.  In the
 414 usual case that the input significand is normalized, `zExp' must be 1 less
 415 than the ``true'' floating-point exponent.  The handling of underflow and
 416 overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 417 -------------------------------------------------------------------------------
 418 */
 419 static float64
 420  roundAndPackFloat64(
 421      flag zSign, int16 zExp, bits32 zSig0, bits32 zSig1, bits32 zSig2 )
 422 {
 423     int8 roundingMode;
 424     flag roundNearestEven, increment, isTiny;
 425
 426     roundingMode = float_rounding_mode;
 427     roundNearestEven = ( roundingMode == float_round_nearest_even );
 428     increment = ( (sbits32) zSig2 < 0 );
 429     if ( ! roundNearestEven ) {
 430         if ( roundingMode == float_round_to_zero ) {
 431             increment = 0;
 432         }
 433         else {
 434             if ( zSign ) {
 435                 increment = ( roundingMode == float_round_down ) && zSig2;
 436             }
 437             else {
 438                 increment = ( roundingMode == float_round_up ) && zSig2;
 439             }
 440         }
 441     }
 442     if ( 0x7FD <= (bits16) zExp ) {
 443         if (    ( 0x7FD < zExp )
 444              || (    ( zExp == 0x7FD )
 445                   && eq64( 0x001FFFFF, 0xFFFFFFFF, zSig0, zSig1 )
 446                   && increment
 447                 )
 448            ) {
 449             float_raise( float_flag_overflow | float_flag_inexact );
 450             if (    ( roundingMode == float_round_to_zero )
 451                  || ( zSign && ( roundingMode == float_round_up ) )
 452                  || ( ! zSign && ( roundingMode == float_round_down ) )
 453                ) {
 454                 return packFloat64( zSign, 0x7FE, 0x000FFFFF, 0xFFFFFFFF );
 455             }
 456             return packFloat64( zSign, 0x7FF, 0, 0 );
 457         }
 458         if ( zExp < 0 ) {
 459             isTiny =
 460                    ( float_detect_tininess == float_tininess_before_rounding )
 461                 || ( zExp < -1 )
 462                 || ! increment
 463                 || lt64( zSig0, zSig1, 0x001FFFFF, 0xFFFFFFFF );
 464             shift64ExtraRightJamming(
 465                 zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );
 466             zExp = 0;
 467             if ( isTiny && zSig2 ) float_raise( float_flag_underflow );
 468             if ( roundNearestEven ) {
 469                 increment = ( (sbits32) zSig2 < 0 );
 470             }
 471             else {
 472                 if ( zSign ) {
 473                     increment = ( roundingMode == float_round_down ) && zSig2;
 474                 }
 475                 else {
 476                     increment = ( roundingMode == float_round_up ) && zSig2;
 477                 }
 478             }
 479         }
 480     }
 481     if ( zSig2 ) set_float_exception_inexact_flag();
 482     if ( increment ) {
 483         add64( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
 484         zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );
 485     }
 486     else {
 487         if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;
 488     }
 489     return packFloat64( zSign, zExp, zSig0, zSig1 );
 490
 491 }
 492
 493 /*
 494 -------------------------------------------------------------------------------
 495 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 496 and significand formed by the concatenation of `zSig0' and `zSig1', and
 497 returns the proper double-precision floating-point value corresponding
 498 to the abstract input.  This routine is just like `roundAndPackFloat64'
 499 except that the input significand has fewer bits and does not have to be
 500 normalized.  In all cases, `zExp' must be 1 less than the ``true'' floating-
 501 point exponent.
 502 -------------------------------------------------------------------------------
 503 */
 504 static float64
 505  normalizeRoundAndPackFloat64(
 506      flag zSign, int16 zExp, bits32 zSig0, bits32 zSig1 )
 507 {
 508     int8 shiftCount;
 509     bits32 zSig2;
 510
 511     if ( zSig0 == 0 ) {
 512         zSig0 = zSig1;
 513         zSig1 = 0;
 514         zExp -= 32;
 515     }
 516     shiftCount = countLeadingZeros32( zSig0 ) - 11;
 517     if ( 0 <= shiftCount ) {
 518         zSig2 = 0;
 519         shortShift64Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
 520     }
 521     else {
 522         shift64ExtraRightJamming(
 523             zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );
 524     }
 525     zExp -= shiftCount;
 526     return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 );
 527
 528 }
 529
 530 /*
 531 -------------------------------------------------------------------------------
 532 Returns the result of converting the 32-bit two's complement integer `a' to
 533 the single-precision floating-point format.  The conversion is performed
 534 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 535 -------------------------------------------------------------------------------
 536 */
 537 float32 int32_to_float32( int32 a )
 538 {
 539     flag zSign;
 540
 541     if ( a == 0 ) return 0;
 542     if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
 543     zSign = ( a < 0 );
 544     return normalizeRoundAndPackFloat32(zSign, 0x9C, (uint32)(zSign ? - a : a));
 545
 546 }
 547
 548 /*
 549 -------------------------------------------------------------------------------
 550 Returns the result of converting the 32-bit two's complement integer `a' to
 551 the double-precision floating-point format.  The conversion is performed
 552 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 553 -------------------------------------------------------------------------------
 554 */
 555 float64 int32_to_float64( int32 a )
 556 {
 557     flag zSign;
 558     bits32 absA;
 559     int8 shiftCount;
 560     bits32 zSig0, zSig1;
 561
 562     if ( a == 0 ) return packFloat64( 0, 0, 0, 0 );
 563     zSign = ( a < 0 );
 564     absA = zSign ? - a : a;
 565     shiftCount = countLeadingZeros32( absA ) - 11;
 566     if ( 0 <= shiftCount ) {
 567         zSig0 = absA<<shiftCount;
 568         zSig1 = 0;
 569     }
 570     else {
 571         shift64Right( absA, 0, - shiftCount, &zSig0, &zSig1 );
 572     }
 573     return packFloat64( zSign, 0x412 - shiftCount, zSig0, zSig1 );
 574
 575 }
 576
 577 #ifndef SOFTFLOAT_FOR_GCC
 578 /*
 579 -------------------------------------------------------------------------------
 580 Returns the result of converting the single-precision floating-point value
 581 `a' to the 32-bit two's complement integer format.  The conversion is
 582 performed according to the IEC/IEEE Standard for Binary Floating-Point
 583 Arithmetic---which means in particular that the conversion is rounded
 584 according to the current rounding mode.  If `a' is a NaN, the largest
 585 positive integer is returned.  Otherwise, if the conversion overflows, the
 586 largest integer with the same sign as `a' is returned.
 587 -------------------------------------------------------------------------------
 588 */
 589 int32 float32_to_int32( float32 a )
 590 {
 591     flag aSign;
 592     int16 aExp, shiftCount;
 593     bits32 aSig, aSigExtra;
 594     int32 z;
 595     int8 roundingMode;
 596
 597     aSig = extractFloat32Frac( a );
 598     aExp = extractFloat32Exp( a );
 599     aSign = extractFloat32Sign( a );
 600     shiftCount = aExp - 0x96;
 601     if ( 0 <= shiftCount ) {
 602         if ( 0x9E <= aExp ) {
 603             if ( a != 0xCF000000 ) {
 604                 float_raise( float_flag_invalid );
 605                 if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
 606                     return 0x7FFFFFFF;
 607                 }
 608             }
 609             return (sbits32) 0x80000000;
 610         }
 611         z = ( aSig | 0x00800000 )<<shiftCount;
 612         if ( aSign ) z = - z;
 613     }
 614     else {
 615         if ( aExp < 0x7E ) {
 616             aSigExtra = aExp | aSig;
 617             z = 0;
 618         }
 619         else {
 620             aSig |= 0x00800000;
 621             aSigExtra = aSig<<( shiftCount & 31 );
 622             z = aSig>>( - shiftCount );
 623         }
 624         if ( aSigExtra ) set_float_exception_inexact_flag();
 625         roundingMode = float_rounding_mode;
 626         if ( roundingMode == float_round_nearest_even ) {
 627             if ( (sbits32) aSigExtra < 0 ) {
 628                 ++z;
 629                 if ( (bits32) ( aSigExtra<<1 ) == 0 ) z &= ~1;
 630             }
 631             if ( aSign ) z = - z;
 632         }
 633         else {
 634             aSigExtra = ( aSigExtra != 0 );
 635             if ( aSign ) {
 636                 z += ( roundingMode == float_round_down ) & aSigExtra;
 637                 z = - z;
 638             }
 639             else {
 640                 z += ( roundingMode == float_round_up ) & aSigExtra;
 641             }
 642         }
 643     }
 644     return z;
 645
 646 }
 647 #endif
 648
 649 /*
 650 -------------------------------------------------------------------------------
 651 Returns the result of converting the single-precision floating-point value
 652 `a' to the 32-bit two's complement integer format.  The conversion is
 653 performed according to the IEC/IEEE Standard for Binary Floating-Point
 654 Arithmetic, except that the conversion is always rounded toward zero.
 655 If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
 656 the conversion overflows, the largest integer with the same sign as `a' is
 657 returned.
 658 -------------------------------------------------------------------------------
 659 */
 660 int32 float32_to_int32_round_to_zero( float32 a )
 661 {
 662     flag aSign;
 663     int16 aExp, shiftCount;
 664     bits32 aSig;
 665     int32 z;
 666
 667     aSig = extractFloat32Frac( a );
 668     aExp = extractFloat32Exp( a );
 669     aSign = extractFloat32Sign( a );
 670     shiftCount = aExp - 0x9E;
 671     if ( 0 <= shiftCount ) {
 672         if ( a != 0xCF000000 ) {
 673             float_raise( float_flag_invalid );
 674             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
 675         }
 676         return (sbits32) 0x80000000;
 677     }
 678     else if ( aExp <= 0x7E ) {
 679         if ( aExp | aSig ) set_float_exception_inexact_flag();
 680         return 0;
 681     }
 682     aSig = ( aSig | 0x00800000 )<<8;
 683     z = aSig>>( - shiftCount );
 684     if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
 685         set_float_exception_inexact_flag();
 686     }
 687     if ( aSign ) z = - z;
 688     return z;
 689
 690 }
 691
 692 /*
 693 -------------------------------------------------------------------------------
 694 Returns the result of converting the single-precision floating-point value
 695 `a' to the double-precision floating-point format.  The conversion is
 696 performed according to the IEC/IEEE Standard for Binary Floating-Point
 697 Arithmetic.
 698 -------------------------------------------------------------------------------
 699 */
 700 float64 float32_to_float64( float32 a )
 701 {
 702     flag aSign;
 703     int16 aExp;
 704     bits32 aSig, zSig0, zSig1;
 705
 706     aSig = extractFloat32Frac( a );
 707     aExp = extractFloat32Exp( a );
 708     aSign = extractFloat32Sign( a );
 709     if ( aExp == 0xFF ) {
 710         if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) );
 711         return packFloat64( aSign, 0x7FF, 0, 0 );
 712     }
 713     if ( aExp == 0 ) {
 714         if ( aSig == 0 ) return packFloat64( aSign, 0, 0, 0 );
 715         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 716         --aExp;
 717     }
 718     shift64Right( aSig, 0, 3, &zSig0, &zSig1 );
 719     return packFloat64( aSign, aExp + 0x380, zSig0, zSig1 );
 720
 721 }
 722
 723 #ifndef SOFTFLOAT_FOR_GCC
 724 /*
 725 -------------------------------------------------------------------------------
 726 Rounds the single-precision floating-point value `a' to an integer,
 727 and returns the result as a single-precision floating-point value.  The
 728 operation is performed according to the IEC/IEEE Standard for Binary
 729 Floating-Point Arithmetic.
 730 -------------------------------------------------------------------------------
 731 */
 732 float32 float32_round_to_int( float32 a )
 733 {
 734     flag aSign;
 735     int16 aExp;
 736     bits32 lastBitMask, roundBitsMask;
 737     int8 roundingMode;
 738     float32 z;
 739
 740     aExp = extractFloat32Exp( a );
 741     if ( 0x96 <= aExp ) {
 742         if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
 743             return propagateFloat32NaN( a, a );
 744         }
 745         return a;
 746     }
 747     if ( aExp <= 0x7E ) {
 748         if ( (bits32) ( a<<1 ) == 0 ) return a;
 749         set_float_exception_inexact_flag();
 750         aSign = extractFloat32Sign( a );
 751         switch ( float_rounding_mode ) {
 752          case float_round_nearest_even:
 753             if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
 754                 return packFloat32( aSign, 0x7F, 0 );
 755             }
 756             break;
 757          case float_round_to_zero:
 758             break;
 759          case float_round_down:
 760             return aSign ? 0xBF800000 : 0;
 761          case float_round_up:
 762             return aSign ? 0x80000000 : 0x3F800000;
 763         }
 764         return packFloat32( aSign, 0, 0 );
 765     }
 766     lastBitMask = 1;
 767     lastBitMask <<= 0x96 - aExp;
 768     roundBitsMask = lastBitMask - 1;
 769     z = a;
 770     roundingMode = float_rounding_mode;
 771     if ( roundingMode == float_round_nearest_even ) {
 772         z += lastBitMask>>1;
 773         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
 774     }
 775     else if ( roundingMode != float_round_to_zero ) {
 776         if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
 777             z += roundBitsMask;
 778         }
 779     }
 780     z &= ~ roundBitsMask;
 781     if ( z != a ) set_float_exception_inexact_flag();
 782     return z;
 783
 784 }
 785 #endif
 786
 787 /*
 788 -------------------------------------------------------------------------------
 789 Returns the result of adding the absolute values of the single-precision
 790 floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
 791 before being returned.  `zSign' is ignored if the result is a NaN.
 792 The addition is performed according to the IEC/IEEE Standard for Binary
 793 Floating-Point Arithmetic.
 794 -------------------------------------------------------------------------------
 795 */
 796 static float32 addFloat32Sigs( float32 a, float32 b, flag zSign )
 797 {
 798     int16 aExp, bExp, zExp;
 799     bits32 aSig, bSig, zSig;
 800     int16 expDiff;
 801
 802     aSig = extractFloat32Frac( a );
 803     aExp = extractFloat32Exp( a );
 804     bSig = extractFloat32Frac( b );
 805     bExp = extractFloat32Exp( b );
 806     expDiff = aExp - bExp;
 807     aSig <<= 6;
 808     bSig <<= 6;
 809     if ( 0 < expDiff ) {
 810         if ( aExp == 0xFF ) {
 811             if ( aSig ) return propagateFloat32NaN( a, b );
 812             return a;
 813         }
 814         if ( bExp == 0 ) {
 815             --expDiff;
 816         }
 817         else {
 818             bSig |= 0x20000000;
 819         }
 820         shift32RightJamming( bSig, expDiff, &bSig );
 821         zExp = aExp;
 822     }
 823     else if ( expDiff < 0 ) {
 824         if ( bExp == 0xFF ) {
 825             if ( bSig ) return propagateFloat32NaN( a, b );
 826             return packFloat32( zSign, 0xFF, 0 );
 827         }
 828         if ( aExp == 0 ) {
 829             ++expDiff;
 830         }
 831         else {
 832             aSig |= 0x20000000;
 833         }
 834         shift32RightJamming( aSig, - expDiff, &aSig );
 835         zExp = bExp;
 836     }
 837     else {
 838         if ( aExp == 0xFF ) {
 839             if ( aSig | bSig ) return propagateFloat32NaN( a, b );
 840             return a;
 841         }
 842         if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
 843         zSig = 0x40000000 + aSig + bSig;
 844         zExp = aExp;
 845         goto roundAndPack;
 846     }
 847     aSig |= 0x20000000;
 848     zSig = ( aSig + bSig )<<1;
 849     --zExp;
 850     if ( (sbits32) zSig < 0 ) {
 851         zSig = aSig + bSig;
 852         ++zExp;
 853     }
 854  roundAndPack:
 855     return roundAndPackFloat32( zSign, zExp, zSig );
 856
 857 }
 858
 859 /*
 860 -------------------------------------------------------------------------------
 861 Returns the result of subtracting the absolute values of the single-
 862 precision floating-point values `a' and `b'.  If `zSign' is 1, the
 863 difference is negated before being returned.  `zSign' is ignored if the
 864 result is a NaN.  The subtraction is performed according to the IEC/IEEE
 865 Standard for Binary Floating-Point Arithmetic.
 866 -------------------------------------------------------------------------------
 867 */
 868 static float32 subFloat32Sigs( float32 a, float32 b, flag zSign )
 869 {
 870     int16 aExp, bExp, zExp;
 871     bits32 aSig, bSig, zSig;
 872     int16 expDiff;
 873
 874     aSig = extractFloat32Frac( a );
 875     aExp = extractFloat32Exp( a );
 876     bSig = extractFloat32Frac( b );
 877     bExp = extractFloat32Exp( b );
 878     expDiff = aExp - bExp;
 879     aSig <<= 7;
 880     bSig <<= 7;
 881     if ( 0 < expDiff ) goto aExpBigger;
 882     if ( expDiff < 0 ) goto bExpBigger;
 883     if ( aExp == 0xFF ) {
 884         if ( aSig | bSig ) return propagateFloat32NaN( a, b );
 885         float_raise( float_flag_invalid );
 886         return float32_default_nan;
 887     }
 888     if ( aExp == 0 ) {
 889         aExp = 1;
 890         bExp = 1;
 891     }
 892     if ( bSig < aSig ) goto aBigger;
 893     if ( aSig < bSig ) goto bBigger;
 894     return packFloat32( float_rounding_mode == float_round_down, 0, 0 );
 895  bExpBigger:
 896     if ( bExp == 0xFF ) {
 897         if ( bSig ) return propagateFloat32NaN( a, b );
 898         return packFloat32( zSign ^ 1, 0xFF, 0 );
 899     }
 900     if ( aExp == 0 ) {
 901         ++expDiff;
 902     }
 903     else {
 904         aSig |= 0x40000000;
 905     }
 906     shift32RightJamming( aSig, - expDiff, &aSig );
 907     bSig |= 0x40000000;
 908  bBigger:
 909     zSig = bSig - aSig;
 910     zExp = bExp;
 911     zSign ^= 1;
 912     goto normalizeRoundAndPack;
 913  aExpBigger:
 914     if ( aExp == 0xFF ) {
 915         if ( aSig ) return propagateFloat32NaN( a, b );
 916         return a;
 917     }
 918     if ( bExp == 0 ) {
 919         --expDiff;
 920     }
 921     else {
 922         bSig |= 0x40000000;
 923     }
 924     shift32RightJamming( bSig, expDiff, &bSig );
 925     aSig |= 0x40000000;
 926  aBigger:
 927     zSig = aSig - bSig;
 928     zExp = aExp;
 929  normalizeRoundAndPack:
 930     --zExp;
 931     return normalizeRoundAndPackFloat32( zSign, zExp, zSig );
 932
 933 }
 934
 935 /*
 936 -------------------------------------------------------------------------------
 937 Returns the result of adding the single-precision floating-point values `a'
 938 and `b'.  The operation is performed according to the IEC/IEEE Standard for
 939 Binary Floating-Point Arithmetic.
 940 -------------------------------------------------------------------------------
 941 */
 942 float32 float32_add( float32 a, float32 b )
 943 {
 944     flag aSign, bSign;
 945
 946     aSign = extractFloat32Sign( a );
 947     bSign = extractFloat32Sign( b );
 948     if ( aSign == bSign ) {
 949         return addFloat32Sigs( a, b, aSign );
 950     }
 951     else {
 952         return subFloat32Sigs( a, b, aSign );
 953     }
 954
 955 }
 956
 957 /*
 958 -------------------------------------------------------------------------------
 959 Returns the result of subtracting the single-precision floating-point values
 960 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 961 for Binary Floating-Point Arithmetic.
 962 -------------------------------------------------------------------------------
 963 */
 964 float32 float32_sub( float32 a, float32 b )
 965 {
 966     flag aSign, bSign;
 967
 968     aSign = extractFloat32Sign( a );
 969     bSign = extractFloat32Sign( b );
 970     if ( aSign == bSign ) {
 971         return subFloat32Sigs( a, b, aSign );
 972     }
 973     else {
 974         return addFloat32Sigs( a, b, aSign );
 975     }
 976
 977 }
 978
 979 /*
 980 -------------------------------------------------------------------------------
 981 Returns the result of multiplying the single-precision floating-point values
 982 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 983 for Binary Floating-Point Arithmetic.
 984 -------------------------------------------------------------------------------
 985 */
 986 float32 float32_mul( float32 a, float32 b )
 987 {
 988     flag aSign, bSign, zSign;
 989     int16 aExp, bExp, zExp;
 990     bits32 aSig, bSig, zSig0, zSig1;
 991
 992     aSig = extractFloat32Frac( a );
 993     aExp = extractFloat32Exp( a );
 994     aSign = extractFloat32Sign( a );
 995     bSig = extractFloat32Frac( b );
 996     bExp = extractFloat32Exp( b );
 997     bSign = extractFloat32Sign( b );
 998     zSign = aSign ^ bSign;
 999     if ( aExp == 0xFF ) {
1000         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1001             return propagateFloat32NaN( a, b );
1002         }
1003         if ( ( bExp | bSig ) == 0 ) {
1004             float_raise( float_flag_invalid );
1005             return float32_default_nan;
1006         }
1007         return packFloat32( zSign, 0xFF, 0 );
1008     }
1009     if ( bExp == 0xFF ) {
1010         if ( bSig ) return propagateFloat32NaN( a, b );
1011         if ( ( aExp | aSig ) == 0 ) {
1012             float_raise( float_flag_invalid );
1013             return float32_default_nan;
1014         }
1015         return packFloat32( zSign, 0xFF, 0 );
1016     }
1017     if ( aExp == 0 ) {
1018         if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1019         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1020     }
1021     if ( bExp == 0 ) {
1022         if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
1023         normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1024     }
1025     zExp = aExp + bExp - 0x7F;
1026     aSig = ( aSig | 0x00800000 )<<7;
1027     bSig = ( bSig | 0x00800000 )<<8;
1028     mul32To64( aSig, bSig, &zSig0, &zSig1 );
1029     zSig0 |= ( zSig1 != 0 );
1030     if ( 0 <= (sbits32) ( zSig0<<1 ) ) {
1031         zSig0 <<= 1;
1032         --zExp;
1033     }
1034     return roundAndPackFloat32( zSign, zExp, zSig0 );
1035
1036 }
1037
1038 /*
1039 -------------------------------------------------------------------------------
1040 Returns the result of dividing the single-precision floating-point value `a'
1041 by the corresponding value `b'.  The operation is performed according to the
1042 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1043 -------------------------------------------------------------------------------
1044 */
1045 float32 float32_div( float32 a, float32 b )
1046 {
1047     flag aSign, bSign, zSign;
1048     int16 aExp, bExp, zExp;
1049     bits32 aSig, bSig, zSig, rem0, rem1, term0, term1;
1050
1051     aSig = extractFloat32Frac( a );
1052     aExp = extractFloat32Exp( a );
1053     aSign = extractFloat32Sign( a );
1054     bSig = extractFloat32Frac( b );
1055     bExp = extractFloat32Exp( b );
1056     bSign = extractFloat32Sign( b );
1057     zSign = aSign ^ bSign;
1058     if ( aExp == 0xFF ) {
1059         if ( aSig ) return propagateFloat32NaN( a, b );
1060         if ( bExp == 0xFF ) {
1061             if ( bSig ) return propagateFloat32NaN( a, b );
1062             float_raise( float_flag_invalid );
1063             return float32_default_nan;
1064         }
1065         return packFloat32( zSign, 0xFF, 0 );
1066     }
1067     if ( bExp == 0xFF ) {
1068         if ( bSig ) return propagateFloat32NaN( a, b );
1069         return packFloat32( zSign, 0, 0 );
1070     }
1071     if ( bExp == 0 ) {
1072         if ( bSig == 0 ) {
1073             if ( ( aExp | aSig ) == 0 ) {
1074                 float_raise( float_flag_invalid );
1075                 return float32_default_nan;
1076             }
1077             float_raise( float_flag_divbyzero );
1078             return packFloat32( zSign, 0xFF, 0 );
1079         }
1080         normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1081     }
1082     if ( aExp == 0 ) {
1083         if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
1084         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1085     }
1086     zExp = aExp - bExp + 0x7D;
1087     aSig = ( aSig | 0x00800000 )<<7;
1088     bSig = ( bSig | 0x00800000 )<<8;
1089     if ( bSig <= ( aSig + aSig ) ) {
1090         aSig >>= 1;
1091         ++zExp;
1092     }
1093     zSig = estimateDiv64To32( aSig, 0, bSig );
1094     if ( ( zSig & 0x3F ) <= 2 ) {
1095         mul32To64( bSig, zSig, &term0, &term1 );
1096         sub64( aSig, 0, term0, term1, &rem0, &rem1 );
1097         while ( (sbits32) rem0 < 0 ) {
1098             --zSig;
1099             add64( rem0, rem1, 0, bSig, &rem0, &rem1 );
1100         }
1101         zSig |= ( rem1 != 0 );
1102     }
1103     return roundAndPackFloat32( zSign, zExp, zSig );
1104
1105 }
1106
1107 #ifndef SOFTFLOAT_FOR_GCC
1108 /*
1109 -------------------------------------------------------------------------------
1110 Returns the remainder of the single-precision floating-point value `a'
1111 with respect to the corresponding value `b'.  The operation is performed
1112 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1113 -------------------------------------------------------------------------------
1114 */
1115 float32 float32_rem( float32 a, float32 b )
1116 {
1117     flag aSign, bSign, zSign;
1118     int16 aExp, bExp, expDiff;
1119     bits32 aSig, bSig, q, allZero, alternateASig;
1120     sbits32 sigMean;
1121
1122     aSig = extractFloat32Frac( a );
1123     aExp = extractFloat32Exp( a );
1124     aSign = extractFloat32Sign( a );
1125     bSig = extractFloat32Frac( b );
1126     bExp = extractFloat32Exp( b );
1127     bSign = extractFloat32Sign( b );
1128     if ( aExp == 0xFF ) {
1129         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
1130             return propagateFloat32NaN( a, b );
1131         }
1132         float_raise( float_flag_invalid );
1133         return float32_default_nan;
1134     }
1135     if ( bExp == 0xFF ) {
1136         if ( bSig ) return propagateFloat32NaN( a, b );
1137         return a;
1138     }
1139     if ( bExp == 0 ) {
1140         if ( bSig == 0 ) {
1141             float_raise( float_flag_invalid );
1142             return float32_default_nan;
1143         }
1144         normalizeFloat32Subnormal( bSig, &bExp, &bSig );
1145     }
1146     if ( aExp == 0 ) {
1147         if ( aSig == 0 ) return a;
1148         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1149     }
1150     expDiff = aExp - bExp;
1151     aSig = ( aSig | 0x00800000 )<<8;
1152     bSig = ( bSig | 0x00800000 )<<8;
1153     if ( expDiff < 0 ) {
1154         if ( expDiff < -1 ) return a;
1155         aSig >>= 1;
1156     }
1157     q = ( bSig <= aSig );
1158     if ( q ) aSig -= bSig;
1159     expDiff -= 32;
1160     while ( 0 < expDiff ) {
1161         q = estimateDiv64To32( aSig, 0, bSig );
1162         q = ( 2 < q ) ? q - 2 : 0;
1163         aSig = - ( ( bSig>>2 ) * q );
1164         expDiff -= 30;
1165     }
1166     expDiff += 32;
1167     if ( 0 < expDiff ) {
1168         q = estimateDiv64To32( aSig, 0, bSig );
1169         q = ( 2 < q ) ? q - 2 : 0;
1170         q >>= 32 - expDiff;
1171         bSig >>= 2;
1172         aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
1173     }
1174     else {
1175         aSig >>= 2;
1176         bSig >>= 2;
1177     }
1178     do {
1179         alternateASig = aSig;
1180         ++q;
1181         aSig -= bSig;
1182     } while ( 0 <= (sbits32) aSig );
1183     sigMean = aSig + alternateASig;
1184     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
1185         aSig = alternateASig;
1186     }
1187     zSign = ( (sbits32) aSig < 0 );
1188     if ( zSign ) aSig = - aSig;
1189     return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig );
1190
1191 }
1192 #endif
1193
1194 #ifndef SOFTFLOAT_FOR_GCC
1195 /*
1196 -------------------------------------------------------------------------------
1197 Returns the square root of the single-precision floating-point value `a'.
1198 The operation is performed according to the IEC/IEEE Standard for Binary
1199 Floating-Point Arithmetic.
1200 -------------------------------------------------------------------------------
1201 */
1202 float32 float32_sqrt( float32 a )
1203 {
1204     flag aSign;
1205     int16 aExp, zExp;
1206     bits32 aSig, zSig, rem0, rem1, term0, term1;
1207
1208     aSig = extractFloat32Frac( a );
1209     aExp = extractFloat32Exp( a );
1210     aSign = extractFloat32Sign( a );
1211     if ( aExp == 0xFF ) {
1212         if ( aSig ) return propagateFloat32NaN( a, 0 );
1213         if ( ! aSign ) return a;
1214         float_raise( float_flag_invalid );
1215         return float32_default_nan;
1216     }
1217     if ( aSign ) {
1218         if ( ( aExp | aSig ) == 0 ) return a;
1219         float_raise( float_flag_invalid );
1220         return float32_default_nan;
1221     }
1222     if ( aExp == 0 ) {
1223         if ( aSig == 0 ) return 0;
1224         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
1225     }
1226     zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
1227     aSig = ( aSig | 0x00800000 )<<8;
1228     zSig = estimateSqrt32( aExp, aSig ) + 2;
1229     if ( ( zSig & 0x7F ) <= 5 ) {
1230         if ( zSig < 2 ) {
1231             zSig = 0x7FFFFFFF;
1232             goto roundAndPack;
1233         }
1234         else {
1235             aSig >>= aExp & 1;
1236             mul32To64( zSig, zSig, &term0, &term1 );
1237             sub64( aSig, 0, term0, term1, &rem0, &rem1 );
1238             while ( (sbits32) rem0 < 0 ) {
1239                 --zSig;
1240                 shortShift64Left( 0, zSig, 1, &term0, &term1 );
1241                 term1 |= 1;
1242                 add64( rem0, rem1, term0, term1, &rem0, &rem1 );
1243             }
1244             zSig |= ( ( rem0 | rem1 ) != 0 );
1245         }
1246     }
1247     shift32RightJamming( zSig, 1, &zSig );
1248  roundAndPack:
1249     return roundAndPackFloat32( 0, zExp, zSig );
1250
1251 }
1252 #endif
1253
1254 /*
1255 -------------------------------------------------------------------------------
1256 Returns 1 if the single-precision floating-point value `a' is equal to
1257 the corresponding value `b', and 0 otherwise.  The comparison is performed
1258 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1259 -------------------------------------------------------------------------------
1260 */
1261 flag float32_eq( float32 a, float32 b )
1262 {
1263
1264     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
1265          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
1266        ) {
1267         if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
1268             float_raise( float_flag_invalid );
1269         }
1270         return 0;
1271     }
1272     return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
1273
1274 }
1275
1276 /*
1277 -------------------------------------------------------------------------------
1278 Returns 1 if the single-precision floating-point value `a' is less than
1279 or equal to the corresponding value `b', and 0 otherwise.  The comparison
1280 is performed according to the IEC/IEEE Standard for Binary Floating-Point
1281 Arithmetic.
1282 -------------------------------------------------------------------------------
1283 */
1284 flag float32_le( float32 a, float32 b )
1285 {
1286     flag aSign, bSign;
1287
1288     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
1289          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
1290        ) {
1291         float_raise( float_flag_invalid );
1292         return 0;
1293     }
1294     aSign = extractFloat32Sign( a );
1295     bSign = extractFloat32Sign( b );
1296     if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
1297     return ( a == b ) || ( aSign ^ ( a < b ) );
1298
1299 }
1300
1301 /*
1302 -------------------------------------------------------------------------------
1303 Returns 1 if the single-precision floating-point value `a' is less than
1304 the corresponding value `b', and 0 otherwise.  The comparison is performed
1305 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1306 -------------------------------------------------------------------------------
1307 */
1308 flag float32_lt( float32 a, float32 b )
1309 {
1310     flag aSign, bSign;
1311
1312     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
1313          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
1314        ) {
1315         float_raise( float_flag_invalid );
1316         return 0;
1317     }
1318     aSign = extractFloat32Sign( a );
1319     bSign = extractFloat32Sign( b );
1320     if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
1321     return ( a != b ) && ( aSign ^ ( a < b ) );
1322
1323 }
1324
1325 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
1326 /*
1327 -------------------------------------------------------------------------------
1328 Returns 1 if the single-precision floating-point value `a' is equal to
1329 the corresponding value `b', and 0 otherwise.  The invalid exception is
1330 raised if either operand is a NaN.  Otherwise, the comparison is performed
1331 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1332 -------------------------------------------------------------------------------
1333 */
1334 flag float32_eq_signaling( float32 a, float32 b )
1335 {
1336
1337     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
1338          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
1339        ) {
1340         float_raise( float_flag_invalid );
1341         return 0;
1342     }
1343     return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
1344
1345 }
1346
1347 /*
1348 -------------------------------------------------------------------------------
1349 Returns 1 if the single-precision floating-point value `a' is less than or
1350 equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
1351 cause an exception.  Otherwise, the comparison is performed according to the
1352 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1353 -------------------------------------------------------------------------------
1354 */
1355 flag float32_le_quiet( float32 a, float32 b )
1356 {
1357     flag aSign, bSign;
1358     int16 aExp, bExp;
1359
1360     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
1361          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
1362        ) {
1363         if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
1364             float_raise( float_flag_invalid );
1365         }
1366         return 0;
1367     }
1368     aSign = extractFloat32Sign( a );
1369     bSign = extractFloat32Sign( b );
1370     if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
1371     return ( a == b ) || ( aSign ^ ( a < b ) );
1372
1373 }
1374
1375 /*
1376 -------------------------------------------------------------------------------
1377 Returns 1 if the single-precision floating-point value `a' is less than
1378 the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
1379 exception.  Otherwise, the comparison is performed according to the IEC/IEEE
1380 Standard for Binary Floating-Point Arithmetic.
1381 -------------------------------------------------------------------------------
1382 */
1383 flag float32_lt_quiet( float32 a, float32 b )
1384 {
1385     flag aSign, bSign;
1386
1387     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
1388          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
1389        ) {
1390         if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
1391             float_raise( float_flag_invalid );
1392         }
1393         return 0;
1394     }
1395     aSign = extractFloat32Sign( a );
1396     bSign = extractFloat32Sign( b );
1397     if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
1398     return ( a != b ) && ( aSign ^ ( a < b ) );
1399
1400 }
1401 #endif /* !SOFTFLOAT_FOR_GCC */
1402
1403 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
1404 /*
1405 -------------------------------------------------------------------------------
1406 Returns the result of converting the double-precision floating-point value
1407 `a' to the 32-bit two's complement integer format.  The conversion is
1408 performed according to the IEC/IEEE Standard for Binary Floating-Point
1409 Arithmetic---which means in particular that the conversion is rounded
1410 according to the current rounding mode.  If `a' is a NaN, the largest
1411 positive integer is returned.  Otherwise, if the conversion overflows, the
1412 largest integer with the same sign as `a' is returned.
1413 -------------------------------------------------------------------------------
1414 */
1415 int32 float64_to_int32( float64 a )
1416 {
1417     flag aSign;
1418     int16 aExp, shiftCount;
1419     bits32 aSig0, aSig1, absZ, aSigExtra;
1420     int32 z;
1421     int8 roundingMode;
1422
1423     aSig1 = extractFloat64Frac1( a );
1424     aSig0 = extractFloat64Frac0( a );
1425     aExp = extractFloat64Exp( a );
1426     aSign = extractFloat64Sign( a );
1427     shiftCount = aExp - 0x413;
1428     if ( 0 <= shiftCount ) {
1429         if ( 0x41E < aExp ) {
1430             if ( ( aExp == 0x7FF ) && ( aSig0 | aSig1 ) ) aSign = 0;
1431             goto invalid;
1432         }
1433         shortShift64Left(
1434             aSig0 | 0x00100000, aSig1, shiftCount, &absZ, &aSigExtra );
1435         if ( 0x80000000 < absZ ) goto invalid;
1436     }
1437     else {
1438         aSig1 = ( aSig1 != 0 );
1439         if ( aExp < 0x3FE ) {
1440             aSigExtra = aExp | aSig0 | aSig1;
1441             absZ = 0;
1442         }
1443         else {
1444             aSig0 |= 0x00100000;
1445             aSigExtra = ( aSig0<<( shiftCount & 31 ) ) | aSig1;
1446             absZ = aSig0>>( - shiftCount );
1447         }
1448     }
1449     roundingMode = float_rounding_mode;
1450     if ( roundingMode == float_round_nearest_even ) {
1451         if ( (sbits32) aSigExtra < 0 ) {
1452             ++absZ;
1453             if ( (bits32) ( aSigExtra<<1 ) == 0 ) absZ &= ~1;
1454         }
1455         z = aSign ? - absZ : absZ;
1456     }
1457     else {
1458         aSigExtra = ( aSigExtra != 0 );
1459         if ( aSign ) {
1460             z = - (   absZ
1461                     + ( ( roundingMode == float_round_down ) & aSigExtra ) );
1462         }
1463         else {
1464             z = absZ + ( ( roundingMode == float_round_up ) & aSigExtra );
1465         }
1466     }
1467     if ( ( aSign ^ ( z < 0 ) ) && z ) {
1468  invalid:
1469         float_raise( float_flag_invalid );
1470         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
1471     }
1472     if ( aSigExtra ) set_float_exception_inexact_flag();
1473     return z;
1474
1475 }
1476 #endif /* !SOFTFLOAT_FOR_GCC */
1477
1478 /*
1479 -------------------------------------------------------------------------------
1480 Returns the result of converting the double-precision floating-point value
1481 `a' to the 32-bit two's complement integer format.  The conversion is
1482 performed according to the IEC/IEEE Standard for Binary Floating-Point
1483 Arithmetic, except that the conversion is always rounded toward zero.
1484 If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
1485 the conversion overflows, the largest integer with the same sign as `a' is
1486 returned.
1487 -------------------------------------------------------------------------------
1488 */
1489 int32 float64_to_int32_round_to_zero( float64 a )
1490 {
1491     flag aSign;
1492     int16 aExp, shiftCount;
1493     bits32 aSig0, aSig1, absZ, aSigExtra;
1494     int32 z;
1495
1496     aSig1 = extractFloat64Frac1( a );
1497     aSig0 = extractFloat64Frac0( a );
1498     aExp = extractFloat64Exp( a );
1499     aSign = extractFloat64Sign( a );
1500     shiftCount = aExp - 0x413;
1501     if ( 0 <= shiftCount ) {
1502         if ( 0x41E < aExp ) {
1503             if ( ( aExp == 0x7FF ) && ( aSig0 | aSig1 ) ) aSign = 0;
1504             goto invalid;
1505         }
1506         shortShift64Left(
1507             aSig0 | 0x00100000, aSig1, shiftCount, &absZ, &aSigExtra );
1508     }
1509     else {
1510         if ( aExp < 0x3FF ) {
1511             if ( aExp | aSig0 | aSig1 ) {
1512                 set_float_exception_inexact_flag();
1513             }
1514             return 0;
1515         }
1516         aSig0 |= 0x00100000;
1517         aSigExtra = ( aSig0<<( shiftCount & 31 ) ) | aSig1;
1518         absZ = aSig0>>( - shiftCount );
1519     }
1520     z = aSign ? - absZ : absZ;
1521     if ( ( aSign ^ ( z < 0 ) ) && z ) {
1522  invalid:
1523         float_raise( float_flag_invalid );
1524         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
1525     }
1526     if ( aSigExtra ) set_float_exception_inexact_flag();
1527     return z;
1528
1529 }
1530
1531 /*
1532 -------------------------------------------------------------------------------
1533 Returns the result of converting the double-precision floating-point value
1534 `a' to the single-precision floating-point format.  The conversion is
1535 performed according to the IEC/IEEE Standard for Binary Floating-Point
1536 Arithmetic.
1537 -------------------------------------------------------------------------------
1538 */
1539 float32 float64_to_float32( float64 a )
1540 {
1541     flag aSign;
1542     int16 aExp;
1543     bits32 aSig0, aSig1, zSig;
1544     bits32 allZero;
1545
1546     aSig1 = extractFloat64Frac1( a );
1547     aSig0 = extractFloat64Frac0( a );
1548     aExp = extractFloat64Exp( a );
1549     aSign = extractFloat64Sign( a );
1550     if ( aExp == 0x7FF ) {
1551         if ( aSig0 | aSig1 ) {
1552             return commonNaNToFloat32( float64ToCommonNaN( a ) );
1553         }
1554         return packFloat32( aSign, 0xFF, 0 );
1555     }
1556     shift64RightJamming( aSig0, aSig1, 22, &allZero, &zSig );
1557     if ( aExp ) zSig |= 0x40000000;
1558     return roundAndPackFloat32( aSign, aExp - 0x381, zSig );
1559
1560 }
1561
1562 #ifndef SOFTFLOAT_FOR_GCC
1563 /*
1564 -------------------------------------------------------------------------------
1565 Rounds the double-precision floating-point value `a' to an integer,
1566 and returns the result as a double-precision floating-point value.  The
1567 operation is performed according to the IEC/IEEE Standard for Binary
1568 Floating-Point Arithmetic.
1569 -------------------------------------------------------------------------------
1570 */
1571 float64 float64_round_to_int( float64 a )
1572 {
1573     flag aSign;
1574     int16 aExp;
1575     bits32 lastBitMask, roundBitsMask;
1576     int8 roundingMode;
1577     float64 z;
1578
1579     aExp = extractFloat64Exp( a );
1580     if ( 0x413 <= aExp ) {
1581         if ( 0x433 <= aExp ) {
1582             if (    ( aExp == 0x7FF )
1583                  && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) ) {
1584                 return propagateFloat64NaN( a, a );
1585             }
1586             return a;
1587         }
1588         lastBitMask = 1;
1589         lastBitMask = ( lastBitMask<<( 0x432 - aExp ) )<<1;
1590         roundBitsMask = lastBitMask - 1;
1591         z = a;
1592         roundingMode = float_rounding_mode;
1593         if ( roundingMode == float_round_nearest_even ) {
1594             if ( lastBitMask ) {
1595                 add64( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
1596                 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
1597             }
1598             else {
1599                 if ( (sbits32) z.low < 0 ) {
1600                     ++z.high;
1601                     if ( (bits32) ( z.low<<1 ) == 0 ) z.high &= ~1;
1602                 }
1603             }
1604         }
1605         else if ( roundingMode != float_round_to_zero ) {
1606             if (   extractFloat64Sign( z )
1607                  ^ ( roundingMode == float_round_up ) ) {
1608                 add64( z.high, z.low, 0, roundBitsMask, &z.high, &z.low );
1609             }
1610         }
1611         z.low &= ~ roundBitsMask;
1612     }
1613     else {
1614         if ( aExp <= 0x3FE ) {
1615             if ( ( ( (bits32) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
1616             set_float_exception_inexact_flag();
1617             aSign = extractFloat64Sign( a );
1618             switch ( float_rounding_mode ) {
1619              case float_round_nearest_even:
1620                 if (    ( aExp == 0x3FE )
1621                      && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) )
1622                    ) {
1623                     return packFloat64( aSign, 0x3FF, 0, 0 );
1624                 }
1625                 break;
1626              case float_round_down:
1627                 return
1628                       aSign ? packFloat64( 1, 0x3FF, 0, 0 )
1629                     : packFloat64( 0, 0, 0, 0 );
1630              case float_round_up:
1631                 return
1632                       aSign ? packFloat64( 1, 0, 0, 0 )
1633                     : packFloat64( 0, 0x3FF, 0, 0 );
1634             }
1635             return packFloat64( aSign, 0, 0, 0 );
1636         }
1637         lastBitMask = 1;
1638         lastBitMask <<= 0x413 - aExp;
1639         roundBitsMask = lastBitMask - 1;
1640         z.low = 0;
1641         z.high = a.high;
1642         roundingMode = float_rounding_mode;
1643         if ( roundingMode == float_round_nearest_even ) {
1644             z.high += lastBitMask>>1;
1645             if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
1646                 z.high &= ~ lastBitMask;
1647             }
1648         }
1649         else if ( roundingMode != float_round_to_zero ) {
1650             if (   extractFloat64Sign( z )
1651                  ^ ( roundingMode == float_round_up ) ) {
1652                 z.high |= ( a.low != 0 );
1653                 z.high += roundBitsMask;
1654             }
1655         }
1656         z.high &= ~ roundBitsMask;
1657     }
1658     if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
1659         set_float_exception_inexact_flag();
1660     }
1661     return z;
1662
1663 }
1664 #endif
1665
1666 /*
1667 -------------------------------------------------------------------------------
1668 Returns the result of adding the absolute values of the double-precision
1669 floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
1670 before being returned.  `zSign' is ignored if the result is a NaN.
1671 The addition is performed according to the IEC/IEEE Standard for Binary
1672 Floating-Point Arithmetic.
1673 -------------------------------------------------------------------------------
1674 */
1675 static float64 addFloat64Sigs( float64 a, float64 b, flag zSign )
1676 {
1677     int16 aExp, bExp, zExp;
1678     bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
1679     int16 expDiff;
1680
1681     aSig1 = extractFloat64Frac1( a );
1682     aSig0 = extractFloat64Frac0( a );
1683     aExp = extractFloat64Exp( a );
1684     bSig1 = extractFloat64Frac1( b );
1685     bSig0 = extractFloat64Frac0( b );
1686     bExp = extractFloat64Exp( b );
1687     expDiff = aExp - bExp;
1688     if ( 0 < expDiff ) {
1689         if ( aExp == 0x7FF ) {
1690             if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, b );
1691             return a;
1692         }
1693         if ( bExp == 0 ) {
1694             --expDiff;
1695         }
1696         else {
1697             bSig0 |= 0x00100000;
1698         }
1699         shift64ExtraRightJamming(
1700             bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
1701         zExp = aExp;
1702     }
1703     else if ( expDiff < 0 ) {
1704         if ( bExp == 0x7FF ) {
1705             if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b );
1706             return packFloat64( zSign, 0x7FF, 0, 0 );
1707         }
1708         if ( aExp == 0 ) {
1709             ++expDiff;
1710         }
1711         else {
1712             aSig0 |= 0x00100000;
1713         }
1714         shift64ExtraRightJamming(
1715             aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
1716         zExp = bExp;
1717     }
1718     else {
1719         if ( aExp == 0x7FF ) {
1720             if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
1721                 return propagateFloat64NaN( a, b );
1722             }
1723             return a;
1724         }
1725         add64( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
1726         if ( aExp == 0 ) return packFloat64( zSign, 0, zSig0, zSig1 );
1727         zSig2 = 0;
1728         zSig0 |= 0x00200000;
1729         zExp = aExp;
1730         goto shiftRight1;
1731     }
1732     aSig0 |= 0x00100000;
1733     add64( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
1734     --zExp;
1735     if ( zSig0 < 0x00200000 ) goto roundAndPack;
1736     ++zExp;
1737  shiftRight1:
1738     shift64ExtraRightJamming( zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
1739  roundAndPack:
1740     return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 );
1741
1742 }
1743
1744 /*
1745 -------------------------------------------------------------------------------
1746 Returns the result of subtracting the absolute values of the double-
1747 precision floating-point values `a' and `b'.  If `zSign' is 1, the
1748 difference is negated before being returned.  `zSign' is ignored if the
1749 result is a NaN.  The subtraction is performed according to the IEC/IEEE
1750 Standard for Binary Floating-Point Arithmetic.
1751 -------------------------------------------------------------------------------
1752 */
1753 static float64 subFloat64Sigs( float64 a, float64 b, flag zSign )
1754 {
1755     int16 aExp, bExp, zExp;
1756     bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
1757     int16 expDiff;
1758
1759     aSig1 = extractFloat64Frac1( a );
1760     aSig0 = extractFloat64Frac0( a );
1761     aExp = extractFloat64Exp( a );
1762     bSig1 = extractFloat64Frac1( b );
1763     bSig0 = extractFloat64Frac0( b );
1764     bExp = extractFloat64Exp( b );
1765     expDiff = aExp - bExp;
1766     shortShift64Left( aSig0, aSig1, 10, &aSig0, &aSig1 );
1767     shortShift64Left( bSig0, bSig1, 10, &bSig0, &bSig1 );
1768     if ( 0 < expDiff ) goto aExpBigger;
1769     if ( expDiff < 0 ) goto bExpBigger;
1770     if ( aExp == 0x7FF ) {
1771         if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
1772             return propagateFloat64NaN( a, b );
1773         }
1774         float_raise( float_flag_invalid );
1775         return float64_default_nan;
1776     }
1777     if ( aExp == 0 ) {
1778         aExp = 1;
1779         bExp = 1;
1780     }
1781     if ( bSig0 < aSig0 ) goto aBigger;
1782     if ( aSig0 < bSig0 ) goto bBigger;
1783     if ( bSig1 < aSig1 ) goto aBigger;
1784     if ( aSig1 < bSig1 ) goto bBigger;
1785     return packFloat64( float_rounding_mode == float_round_down, 0, 0, 0 );
1786  bExpBigger:
1787     if ( bExp == 0x7FF ) {
1788         if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b );
1789         return packFloat64( zSign ^ 1, 0x7FF, 0, 0 );
1790     }
1791     if ( aExp == 0 ) {
1792         ++expDiff;
1793     }
1794     else {
1795         aSig0 |= 0x40000000;
1796     }
1797     shift64RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
1798     bSig0 |= 0x40000000;
1799  bBigger:
1800     sub64( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
1801     zExp = bExp;
1802     zSign ^= 1;
1803     goto normalizeRoundAndPack;
1804  aExpBigger:
1805     if ( aExp == 0x7FF ) {
1806         if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, b );
1807         return a;
1808     }
1809     if ( bExp == 0 ) {
1810         --expDiff;
1811     }
1812     else {
1813         bSig0 |= 0x40000000;
1814     }
1815     shift64RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
1816     aSig0 |= 0x40000000;
1817  aBigger:
1818     sub64( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
1819     zExp = aExp;
1820  normalizeRoundAndPack:
1821     --zExp;
1822     return normalizeRoundAndPackFloat64( zSign, zExp - 10, zSig0, zSig1 );
1823
1824 }
1825
1826 /*
1827 -------------------------------------------------------------------------------
1828 Returns the result of adding the double-precision floating-point values `a'
1829 and `b'.  The operation is performed according to the IEC/IEEE Standard for
1830 Binary Floating-Point Arithmetic.
1831 -------------------------------------------------------------------------------
1832 */
1833 float64 float64_add( float64 a, float64 b )
1834 {
1835     flag aSign, bSign;
1836
1837     aSign = extractFloat64Sign( a );
1838     bSign = extractFloat64Sign( b );
1839     if ( aSign == bSign ) {
1840         return addFloat64Sigs( a, b, aSign );
1841     }
1842     else {
1843         return subFloat64Sigs( a, b, aSign );
1844     }
1845
1846 }
1847
1848 /*
1849 -------------------------------------------------------------------------------
1850 Returns the result of subtracting the double-precision floating-point values
1851 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
1852 for Binary Floating-Point Arithmetic.
1853 -------------------------------------------------------------------------------
1854 */
1855 float64 float64_sub( float64 a, float64 b )
1856 {
1857     flag aSign, bSign;
1858
1859     aSign = extractFloat64Sign( a );
1860     bSign = extractFloat64Sign( b );
1861     if ( aSign == bSign ) {
1862         return subFloat64Sigs( a, b, aSign );
1863     }
1864     else {
1865         return addFloat64Sigs( a, b, aSign );
1866     }
1867
1868 }
1869
1870 /*
1871 -------------------------------------------------------------------------------
1872 Returns the result of multiplying the double-precision floating-point values
1873 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
1874 for Binary Floating-Point Arithmetic.
1875 -------------------------------------------------------------------------------
1876 */
1877 float64 float64_mul( float64 a, float64 b )
1878 {
1879     flag aSign, bSign, zSign;
1880     int16 aExp, bExp, zExp;
1881     bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
1882
1883     aSig1 = extractFloat64Frac1( a );
1884     aSig0 = extractFloat64Frac0( a );
1885     aExp = extractFloat64Exp( a );
1886     aSign = extractFloat64Sign( a );
1887     bSig1 = extractFloat64Frac1( b );
1888     bSig0 = extractFloat64Frac0( b );
1889     bExp = extractFloat64Exp( b );
1890     bSign = extractFloat64Sign( b );
1891     zSign = aSign ^ bSign;
1892     if ( aExp == 0x7FF ) {
1893         if (    ( aSig0 | aSig1 )
1894              || ( ( bExp == 0x7FF ) && ( bSig0 | bSig1 ) ) ) {
1895             return propagateFloat64NaN( a, b );
1896         }
1897         if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
1898         return packFloat64( zSign, 0x7FF, 0, 0 );
1899     }
1900     if ( bExp == 0x7FF ) {
1901         if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b );
1902         if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
1903  invalid:
1904             float_raise( float_flag_invalid );
1905             return float64_default_nan;
1906         }
1907         return packFloat64( zSign, 0x7FF, 0, 0 );
1908     }
1909     if ( aExp == 0 ) {
1910         if ( ( aSig0 | aSig1 ) == 0 ) return packFloat64( zSign, 0, 0, 0 );
1911         normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
1912     }
1913     if ( bExp == 0 ) {
1914         if ( ( bSig0 | bSig1 ) == 0 ) return packFloat64( zSign, 0, 0, 0 );
1915         normalizeFloat64Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
1916     }
1917     zExp = aExp + bExp - 0x400;
1918     aSig0 |= 0x00100000;
1919     shortShift64Left( bSig0, bSig1, 12, &bSig0, &bSig1 );
1920     mul64To128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
1921     add64( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
1922     zSig2 |= ( zSig3 != 0 );
1923     if ( 0x00200000 <= zSig0 ) {
1924         shift64ExtraRightJamming(
1925             zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
1926         ++zExp;
1927     }
1928     return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 );
1929
1930 }
1931
1932 /*
1933 -------------------------------------------------------------------------------
1934 Returns the result of dividing the double-precision floating-point value `a'
1935 by the corresponding value `b'.  The operation is performed according to the
1936 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
1937 -------------------------------------------------------------------------------
1938 */
1939 float64 float64_div( float64 a, float64 b )
1940 {
1941     flag aSign, bSign, zSign;
1942     int16 aExp, bExp, zExp;
1943     bits32 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
1944     bits32 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
1945
1946     aSig1 = extractFloat64Frac1( a );
1947     aSig0 = extractFloat64Frac0( a );
1948     aExp = extractFloat64Exp( a );
1949     aSign = extractFloat64Sign( a );
1950     bSig1 = extractFloat64Frac1( b );
1951     bSig0 = extractFloat64Frac0( b );
1952     bExp = extractFloat64Exp( b );
1953     bSign = extractFloat64Sign( b );
1954     zSign = aSign ^ bSign;
1955     if ( aExp == 0x7FF ) {
1956         if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, b );
1957         if ( bExp == 0x7FF ) {
1958             if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b );
1959             goto invalid;
1960         }
1961         return packFloat64( zSign, 0x7FF, 0, 0 );
1962     }
1963     if ( bExp == 0x7FF ) {
1964         if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b );
1965         return packFloat64( zSign, 0, 0, 0 );
1966     }
1967     if ( bExp == 0 ) {
1968         if ( ( bSig0 | bSig1 ) == 0 ) {
1969             if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
1970  invalid:
1971                 float_raise( float_flag_invalid );
1972                 return float64_default_nan;
1973             }
1974             float_raise( float_flag_divbyzero );
1975             return packFloat64( zSign, 0x7FF, 0, 0 );
1976         }
1977         normalizeFloat64Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
1978     }
1979     if ( aExp == 0 ) {
1980         if ( ( aSig0 | aSig1 ) == 0 ) return packFloat64( zSign, 0, 0, 0 );
1981         normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
1982     }
1983     zExp = aExp - bExp + 0x3FD;
1984     shortShift64Left( aSig0 | 0x00100000, aSig1, 11, &aSig0, &aSig1 );
1985     shortShift64Left( bSig0 | 0x00100000, bSig1, 11, &bSig0, &bSig1 );
1986     if ( le64( bSig0, bSig1, aSig0, aSig1 ) ) {
1987         shift64Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
1988         ++zExp;
1989     }
1990     zSig0 = estimateDiv64To32( aSig0, aSig1, bSig0 );
1991     mul64By32To96( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
1992     sub96( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
1993     while ( (sbits32) rem0 < 0 ) {
1994         --zSig0;
1995         add96( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
1996     }
1997     zSig1 = estimateDiv64To32( rem1, rem2, bSig0 );
1998     if ( ( zSig1 & 0x3FF ) <= 4 ) {
1999         mul64By32To96( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
2000         sub96( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
2001         while ( (sbits32) rem1 < 0 ) {
2002             --zSig1;
2003             add96( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
2004         }
2005         zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
2006     }
2007     shift64ExtraRightJamming( zSig0, zSig1, 0, 11, &zSig0, &zSig1, &zSig2 );
2008     return roundAndPackFloat64( zSign, zExp, zSig0, zSig1, zSig2 );
2009
2010 }
2011
2012 #ifndef SOFTFLOAT_FOR_GCC
2013 /*
2014 -------------------------------------------------------------------------------
2015 Returns the remainder of the double-precision floating-point value `a'
2016 with respect to the corresponding value `b'.  The operation is performed
2017 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2018 -------------------------------------------------------------------------------
2019 */
2020 float64 float64_rem( float64 a, float64 b )
2021 {
2022     flag aSign, bSign, zSign;
2023     int16 aExp, bExp, expDiff;
2024     bits32 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
2025     bits32 allZero, alternateASig0, alternateASig1, sigMean1;
2026     sbits32 sigMean0;
2027     float64 z;
2028
2029     aSig1 = extractFloat64Frac1( a );
2030     aSig0 = extractFloat64Frac0( a );
2031     aExp = extractFloat64Exp( a );
2032     aSign = extractFloat64Sign( a );
2033     bSig1 = extractFloat64Frac1( b );
2034     bSig0 = extractFloat64Frac0( b );
2035     bExp = extractFloat64Exp( b );
2036     bSign = extractFloat64Sign( b );
2037     if ( aExp == 0x7FF ) {
2038         if (    ( aSig0 | aSig1 )
2039              || ( ( bExp == 0x7FF ) && ( bSig0 | bSig1 ) ) ) {
2040             return propagateFloat64NaN( a, b );
2041         }
2042         goto invalid;
2043     }
2044     if ( bExp == 0x7FF ) {
2045         if ( bSig0 | bSig1 ) return propagateFloat64NaN( a, b );
2046         return a;
2047     }
2048     if ( bExp == 0 ) {
2049         if ( ( bSig0 | bSig1 ) == 0 ) {
2050  invalid:
2051             float_raise( float_flag_invalid );
2052             return float64_default_nan;
2053         }
2054         normalizeFloat64Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
2055     }
2056     if ( aExp == 0 ) {
2057         if ( ( aSig0 | aSig1 ) == 0 ) return a;
2058         normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
2059     }
2060     expDiff = aExp - bExp;
2061     if ( expDiff < -1 ) return a;
2062     shortShift64Left(
2063         aSig0 | 0x00100000, aSig1, 11 - ( expDiff < 0 ), &aSig0, &aSig1 );
2064     shortShift64Left( bSig0 | 0x00100000, bSig1, 11, &bSig0, &bSig1 );
2065     q = le64( bSig0, bSig1, aSig0, aSig1 );
2066     if ( q ) sub64( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
2067     expDiff -= 32;
2068     while ( 0 < expDiff ) {
2069         q = estimateDiv64To32( aSig0, aSig1, bSig0 );
2070         q = ( 4 < q ) ? q - 4 : 0;
2071         mul64By32To96( bSig0, bSig1, q, &term0, &term1, &term2 );
2072         shortShift96Left( term0, term1, term2, 29, &term1, &term2, &allZero );
2073         shortShift64Left( aSig0, aSig1, 29, &aSig0, &allZero );
2074         sub64( aSig0, 0, term1, term2, &aSig0, &aSig1 );
2075         expDiff -= 29;
2076     }
2077     if ( -32 < expDiff ) {
2078         q = estimateDiv64To32( aSig0, aSig1, bSig0 );
2079         q = ( 4 < q ) ? q - 4 : 0;
2080         q >>= - expDiff;
2081         shift64Right( bSig0, bSig1, 8, &bSig0, &bSig1 );
2082         expDiff += 24;
2083         if ( expDiff < 0 ) {
2084             shift64Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
2085         }
2086         else {
2087             shortShift64Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 );
2088         }
2089         mul64By32To96( bSig0, bSig1, q, &term0, &term1, &term2 );
2090         sub64( aSig0, aSig1, term1, term2, &aSig0, &aSig1 );
2091     }
2092     else {
2093         shift64Right( aSig0, aSig1, 8, &aSig0, &aSig1 );
2094         shift64Right( bSig0, bSig1, 8, &bSig0, &bSig1 );
2095     }
2096     do {
2097         alternateASig0 = aSig0;
2098         alternateASig1 = aSig1;
2099         ++q;
2100         sub64( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
2101     } while ( 0 <= (sbits32) aSig0 );
2102     add64(
2103         aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 );
2104     if (    ( sigMean0 < 0 )
2105          || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
2106         aSig0 = alternateASig0;
2107         aSig1 = alternateASig1;
2108     }
2109     zSign = ( (sbits32) aSig0 < 0 );
2110     if ( zSign ) sub64( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
2111     return
2112         normalizeRoundAndPackFloat64( aSign ^ zSign, bExp - 4, aSig0, aSig1 );
2113
2114 }
2115 #endif
2116
2117 #ifndef SOFTFLOAT_FOR_GCC
2118 /*
2119 -------------------------------------------------------------------------------
2120 Returns the square root of the double-precision floating-point value `a'.
2121 The operation is performed according to the IEC/IEEE Standard for Binary
2122 Floating-Point Arithmetic.
2123 -------------------------------------------------------------------------------
2124 */
2125 float64 float64_sqrt( float64 a )
2126 {
2127     flag aSign;
2128     int16 aExp, zExp;
2129     bits32 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
2130     bits32 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
2131     float64 z;
2132
2133     aSig1 = extractFloat64Frac1( a );
2134     aSig0 = extractFloat64Frac0( a );
2135     aExp = extractFloat64Exp( a );
2136     aSign = extractFloat64Sign( a );
2137     if ( aExp == 0x7FF ) {
2138         if ( aSig0 | aSig1 ) return propagateFloat64NaN( a, a );
2139         if ( ! aSign ) return a;
2140         goto invalid;
2141     }
2142     if ( aSign ) {
2143         if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a;
2144  invalid:
2145         float_raise( float_flag_invalid );
2146         return float64_default_nan;
2147     }
2148     if ( aExp == 0 ) {
2149         if ( ( aSig0 | aSig1 ) == 0 ) return packFloat64( 0, 0, 0, 0 );
2150         normalizeFloat64Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
2151     }
2152     zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
2153     aSig0 |= 0x00100000;
2154     shortShift64Left( aSig0, aSig1, 11, &term0, &term1 );
2155     zSig0 = ( estimateSqrt32( aExp, term0 )>>1 ) + 1;
2156     if ( zSig0 == 0 ) zSig0 = 0x7FFFFFFF;
2157     doubleZSig0 = zSig0 + zSig0;
2158     shortShift64Left( aSig0, aSig1, 9 - ( aExp & 1 ), &aSig0, &aSig1 );
2159     mul32To64( zSig0, zSig0, &term0, &term1 );
2160     sub64( aSig0, aSig1, term0, term1, &rem0, &rem1 );
2161     while ( (sbits32) rem0 < 0 ) {
2162         --zSig0;
2163         doubleZSig0 -= 2;
2164         add64( rem0, rem1, 0, doubleZSig0 | 1, &rem0, &rem1 );
2165     }
2166     zSig1 = estimateDiv64To32( rem1, 0, doubleZSig0 );
2167     if ( ( zSig1 & 0x1FF ) <= 5 ) {
2168         if ( zSig1 == 0 ) zSig1 = 1;
2169         mul32To64( doubleZSig0, zSig1, &term1, &term2 );
2170         sub64( rem1, 0, term1, term2, &rem1, &rem2 );
2171         mul32To64( zSig1, zSig1, &term2, &term3 );
2172         sub96( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
2173         while ( (sbits32) rem1 < 0 ) {
2174             --zSig1;
2175             shortShift64Left( 0, zSig1, 1, &term2, &term3 );
2176             term3 |= 1;
2177             term2 |= doubleZSig0;
2178             add96( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
2179         }
2180         zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
2181     }
2182     shift64ExtraRightJamming( zSig0, zSig1, 0, 10, &zSig0, &zSig1, &zSig2 );
2183     return roundAndPackFloat64( 0, zExp, zSig0, zSig1, zSig2 );
2184
2185 }
2186 #endif
2187
2188 /*
2189 -------------------------------------------------------------------------------
2190 Returns 1 if the double-precision floating-point value `a' is equal to
2191 the corresponding value `b', and 0 otherwise.  The comparison is performed
2192 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2193 -------------------------------------------------------------------------------
2194 */
2195 flag float64_eq( float64 a, float64 b )
2196 {
2197
2198     if (    (    ( extractFloat64Exp( a ) == 0x7FF )
2199               && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) )
2200          || (    ( extractFloat64Exp( b ) == 0x7FF )
2201               && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) )
2202        ) {
2203         if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
2204             float_raise( float_flag_invalid );
2205         }
2206         return 0;
2207     }
2208     return ( a == b ) ||
2209         ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) == 0 );
2210
2211 }
2212
2213 /*
2214 -------------------------------------------------------------------------------
2215 Returns 1 if the double-precision floating-point value `a' is less than
2216 or equal to the corresponding value `b', and 0 otherwise.  The comparison
2217 is performed according to the IEC/IEEE Standard for Binary Floating-Point
2218 Arithmetic.
2219 -------------------------------------------------------------------------------
2220 */
2221 flag float64_le( float64 a, float64 b )
2222 {
2223     flag aSign, bSign;
2224
2225     if (    (    ( extractFloat64Exp( a ) == 0x7FF )
2226               && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) )
2227          || (    ( extractFloat64Exp( b ) == 0x7FF )
2228               && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) )
2229        ) {
2230         float_raise( float_flag_invalid );
2231         return 0;
2232     }
2233     aSign = extractFloat64Sign( a );
2234     bSign = extractFloat64Sign( b );
2235     if ( aSign != bSign )
2236         return aSign ||
2237             ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) ==
2238               0 );
2239     return ( a == b ) ||
2240         ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) );
2241 }
2242
2243 /*
2244 -------------------------------------------------------------------------------
2245 Returns 1 if the double-precision floating-point value `a' is less than
2246 the corresponding value `b', and 0 otherwise.  The comparison is performed
2247 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2248 -------------------------------------------------------------------------------
2249 */
2250 flag float64_lt( float64 a, float64 b )
2251 {
2252     flag aSign, bSign;
2253
2254     if (    (    ( extractFloat64Exp( a ) == 0x7FF )
2255               && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) )
2256          || (    ( extractFloat64Exp( b ) == 0x7FF )
2257               && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) )
2258        ) {
2259         float_raise( float_flag_invalid );
2260         return 0;
2261     }
2262     aSign = extractFloat64Sign( a );
2263     bSign = extractFloat64Sign( b );
2264     if ( aSign != bSign )
2265         return aSign &&
2266             ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) !=
2267               0 );
2268     return ( a != b ) &&
2269         ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) );
2270
2271 }
2272
2273 #ifndef SOFTFLOAT_FOR_GCC
2274 /*
2275 -------------------------------------------------------------------------------
2276 Returns 1 if the double-precision floating-point value `a' is equal to
2277 the corresponding value `b', and 0 otherwise.  The invalid exception is
2278 raised if either operand is a NaN.  Otherwise, the comparison is performed
2279 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2280 -------------------------------------------------------------------------------
2281 */
2282 flag float64_eq_signaling( float64 a, float64 b )
2283 {
2284
2285     if (    (    ( extractFloat64Exp( a ) == 0x7FF )
2286               && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) )
2287          || (    ( extractFloat64Exp( b ) == 0x7FF )
2288               && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) )
2289        ) {
2290         float_raise( float_flag_invalid );
2291         return 0;
2292     }
2293     return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
2294
2295 }
2296
2297 /*
2298 -------------------------------------------------------------------------------
2299 Returns 1 if the double-precision floating-point value `a' is less than or
2300 equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
2301 cause an exception.  Otherwise, the comparison is performed according to the
2302 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
2303 -------------------------------------------------------------------------------
2304 */
2305 flag float64_le_quiet( float64 a, float64 b )
2306 {
2307     flag aSign, bSign;
2308
2309     if (    (    ( extractFloat64Exp( a ) == 0x7FF )
2310               && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) )
2311          || (    ( extractFloat64Exp( b ) == 0x7FF )
2312               && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) )
2313        ) {
2314         if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
2315             float_raise( float_flag_invalid );
2316         }
2317         return 0;
2318     }
2319     aSign = extractFloat64Sign( a );
2320     bSign = extractFloat64Sign( b );
2321     if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
2322     return ( a == b ) || ( aSign ^ ( a < b ) );
2323
2324 }
2325
2326 /*
2327 -------------------------------------------------------------------------------
2328 Returns 1 if the double-precision floating-point value `a' is less than
2329 the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
2330 exception.  Otherwise, the comparison is performed according to the IEC/IEEE
2331 Standard for Binary Floating-Point Arithmetic.
2332 -------------------------------------------------------------------------------
2333 */
2334 flag float64_lt_quiet( float64 a, float64 b )
2335 {
2336     flag aSign, bSign;
2337
2338     if (    (    ( extractFloat64Exp( a ) == 0x7FF )
2339               && ( extractFloat64Frac0( a ) | extractFloat64Frac1( a ) ) )
2340          || (    ( extractFloat64Exp( b ) == 0x7FF )
2341               && ( extractFloat64Frac0( b ) | extractFloat64Frac1( b ) ) )
2342        ) {
2343         if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
2344             float_raise( float_flag_invalid );
2345         }
2346         return 0;
2347     }
2348     aSign = extractFloat64Sign( a );
2349     bSign = extractFloat64Sign( b );
2350     if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
2351     return ( a != b ) && ( aSign ^ ( a < b ) );
2352
2353 }
2354
2355 #endif