root/lib/libkern/softfloat.c

DEFINITIONS

1. roundAndPackInt32
2. roundAndPackInt64
3. extractFloat32Frac
4. extractFloat32Exp
5. extractFloat32Sign
6. normalizeFloat32Subnormal
7. packFloat32
8. roundAndPackFloat32
9. normalizeRoundAndPackFloat32
10. extractFloat64Frac
11. extractFloat64Exp
12. extractFloat64Sign
13. normalizeFloat64Subnormal
14. packFloat64
15. roundAndPackFloat64
16. normalizeRoundAndPackFloat64
17. extractFloatx80Frac
18. extractFloatx80Exp
19. extractFloatx80Sign
20. normalizeFloatx80Subnormal
21. packFloatx80
22. roundAndPackFloatx80
23. normalizeRoundAndPackFloatx80
24. extractFloat128Frac1
25. extractFloat128Frac0
26. extractFloat128Exp
27. extractFloat128Sign
28. normalizeFloat128Subnormal
29. packFloat128
30. roundAndPackFloat128
31. normalizeRoundAndPackFloat128
32. int32_to_float32
33. int32_to_float64
34. int32_to_floatx80
35. int32_to_float128
36. int64_to_float32
37. int64_to_float64
38. int64_to_floatx80
39. int64_to_float128
40. float32_to_int32
41. float32_to_int32_round_to_zero
42. float32_to_int64
43. float32_to_int64_round_to_zero
44. float32_to_float64
45. float32_to_floatx80
46. float32_to_float128
47. float32_round_to_int
48. addFloat32Sigs
49. subFloat32Sigs
50. float32_add
51. float32_sub
52. float32_mul
53. float32_div
54. float32_rem
55. float32_sqrt
56. float32_eq
57. float32_le
58. float32_lt
59. float32_eq_signaling
60. float32_le_quiet
61. float32_lt_quiet
62. float64_to_int32
63. float64_to_int32_round_to_zero
64. float64_to_int64
65. float64_to_int64_round_to_zero
66. float64_to_float32
67. float64_to_floatx80
68. float64_to_float128
69. float64_round_to_int
70. addFloat64Sigs
71. subFloat64Sigs
72. float64_add
73. float64_sub
74. float64_mul
75. float64_div
76. float64_rem
77. float64_sqrt
78. float64_eq
79. float64_le
80. float64_lt
81. float64_eq_signaling
82. float64_le_quiet
83. float64_lt_quiet
84. floatx80_to_int32
85. floatx80_to_int32_round_to_zero
86. floatx80_to_int64
87. floatx80_to_int64_round_to_zero
88. floatx80_to_float32
89. floatx80_to_float64
90. floatx80_to_float128
91. floatx80_round_to_int
92. addFloatx80Sigs
93. subFloatx80Sigs
94. floatx80_add
95. floatx80_sub
96. floatx80_mul
97. floatx80_div
98. floatx80_rem
99. floatx80_sqrt
100. floatx80_eq
101. floatx80_le
102. floatx80_lt
103. floatx80_eq_signaling
104. floatx80_le_quiet
105. floatx80_lt_quiet
106. float128_to_int32
107. float128_to_int32_round_to_zero
108. float128_to_int64
109. float128_to_int64_round_to_zero
110. float128_to_float32
111. float128_to_float64
112. float128_to_floatx80
113. float128_round_to_int
114. addFloat128Sigs
115. subFloat128Sigs
116. float128_add
117. float128_sub
118. float128_mul
119. float128_div
120. float128_rem
121. float128_sqrt
122. float128_eq
123. float128_le
124. float128_lt
125. float128_eq_signaling
126. float128_le_quiet
127. float128_lt_quiet
128. float64_to_uint32_round_to_zero
129. float32_to_uint32_round_to_zero

    1 /*      $OpenBSD: softfloat.c,v 1.1 2002/04/28 20:55:14 pvalchev Exp $  */
    2 /*      $NetBSD: softfloat.c,v 1.1 2001/04/26 03:10:47 ross Exp $       */
    4 /*
    5  * This version hacked for use with gcc -msoft-float by bjh21.
    6  * (Mostly a case of #ifdefing out things GCC doesn't need or provides
    7  *  itself).
    8  */
   10 /*
   11  * Things you may want to define:
   12  *
   13  * SOFTFLOAT_FOR_GCC - build only those functions necessary for GCC (with
   14  *   -msoft-float) to work.  Include "softfloat-for-gcc.h" to get them
   15  *   properly renamed.
   16  */
   18 /*
   19 ===============================================================================
   20 
   21 This C source file is part of the SoftFloat IEC/IEEE Floating-point
   22 Arithmetic Package, Release 2a.
   23 
   24 Written by John R. Hauser.  This work was made possible in part by the
   25 International Computer Science Institute, located at Suite 600, 1947 Center
   26 Street, Berkeley, California 94704.  Funding was partially provided by the
   27 National Science Foundation under grant MIP-9311980.  The original version
   28 of this code was written as part of a project to build a fixed-point vector
   29 processor in collaboration with the University of California at Berkeley,
   30 overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
   31 is available through the Web page `http://HTTP.CS.Berkeley.EDU/~jhauser/
   32 arithmetic/SoftFloat.html'.
   33 
   34 THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable
   35 effort has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT
   36 WILL AT TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS
   37 RESTRICTED TO PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL
   38 RESPONSIBILITY FOR ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM
   39 THEIR OWN USE OF THE SOFTWARE, AND WHO ALSO EFFECTIVELY INDEMNIFY
   40 (possibly via similar legal warning) JOHN HAUSER AND THE INTERNATIONAL
   41 COMPUTER SCIENCE INSTITUTE AGAINST ALL LOSSES, COSTS, OR OTHER PROBLEMS
   42 ARISING FROM THE USE OF THE SOFTWARE BY THEIR CUSTOMERS AND CLIENTS.
   43 
   44 Derivative works are acceptable, even for commercial purposes, so long as
   45 (1) they include prominent notice that the work is derivative, and (2) they
   46 include prominent notice akin to these four paragraphs for those parts of
   47 this code that are retained.
   48 
   49 ===============================================================================
   50 */
   52 #ifndef NO_IEEE
   54 #include <sys/cdefs.h>
   55 #if defined(LIBC_SCCS) && !defined(lint)
   56 __RCSID("$NetBSD: softfloat.c,v 1.1 2001/04/26 03:10:47 ross Exp $");
   57 #endif /* LIBC_SCCS and not lint */
   59 #ifdef SOFTFLOAT_FOR_GCC
   60 #include "softfloat-for-gcc.h"
   61 #endif
   63 #include "milieu.h"
   64 #include "softfloat.h"
   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
   77 /*
   78 -------------------------------------------------------------------------------
   79 Floating-point rounding mode, extended double-precision rounding precision,
   80 and exception flags.
   81 -------------------------------------------------------------------------------
   82 */
   84 /*
   85  * XXX: This may cause options-MULTIPROCESSOR or thread problems someday.
   86  *      Right now, it does not.  I've removed all other dynamic global
   87  *      variables. [ross]
   88  */
   89 #ifdef FLOATX80
int8 floatx80_rounding_precision = 80;
   91 #endif
   93 /*
   94 -------------------------------------------------------------------------------
   95 Primitive arithmetic functions, including multi-word arithmetic, and
   96 division and square root approximations.  (Can be specialized to target if
   97 desired.)
   98 -------------------------------------------------------------------------------
   99 */
  100 #include "softfloat-macros.h"
  102 /*
  103 -------------------------------------------------------------------------------
  104 Functions and definitions to determine:  (1) whether tininess for underflow
  105 is detected before or after rounding by default, (2) what (if anything)
  106 happens when exceptions are raised, (3) how signaling NaNs are distinguished
  107 from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
  108 are propagated from function inputs to output.  These details are target-
  109 specific.
  110 -------------------------------------------------------------------------------
  111 */
  112 #include "softfloat-specialize.h"
  114 #ifndef SOFTFLOAT_FOR_GCC /* Not used */
  115 /*
  116 -------------------------------------------------------------------------------
  117 Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
  118 and 7, and returns the properly rounded 32-bit integer corresponding to the
  119 input.  If `zSign' is 1, the input is negated before being converted to an
  120 integer.  Bit 63 of `absZ' must be zero.  Ordinarily, the fixed-point input
  121 is simply rounded to an integer, with the inexact exception raised if the
  122 input cannot be represented exactly as an integer.  However, if the fixed-
  123 point input is too large, the invalid exception is raised and the largest
  124 positive or negative integer is returned.
  125 -------------------------------------------------------------------------------
  126 */
  127 static int32 roundAndPackInt32( flag zSign, bits64 absZ )
  128 {
int8 roundingMode;
flag roundNearestEven;
int8 roundIncrement, roundBits;
int32 z;
roundingMode = float_rounding_mode();
roundNearestEven = ( roundingMode == float_round_nearest_even );
roundIncrement = 0x40;
  137     if ( ! roundNearestEven ) {
  138         if ( roundingMode == float_round_to_zero ) {
roundIncrement = 0;
  140         }
  141         else {
roundIncrement = 0x7F;
  143             if ( zSign ) {
  144                 if ( roundingMode == float_round_up ) roundIncrement = 0;
  145             }
  146             else {
  147                 if ( roundingMode == float_round_down ) roundIncrement = 0;
  148             }
  149         }
  150     }
roundBits = absZ & 0x7F;
absZ = ( absZ + roundIncrement )>>7;
absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
z = absZ;
  155     if ( zSign ) z = - z;
  156     if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
float_raise( float_flag_invalid );
  158         return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
  159     }
  160     if ( roundBits ) float_set_inexact();
  161     return z;
  163 }
  165 /*
  166 -------------------------------------------------------------------------------
  167 Takes the 128-bit fixed-point value formed by concatenating `absZ0' and
  168 `absZ1', with binary point between bits 63 and 64 (between the input words),
  169 and returns the properly rounded 64-bit integer corresponding to the input.
  170 If `zSign' is 1, the input is negated before being converted to an integer.
  171 Ordinarily, the fixed-point input is simply rounded to an integer, with
  172 the inexact exception raised if the input cannot be represented exactly as
  173 an integer.  However, if the fixed-point input is too large, the invalid
  174 exception is raised and the largest positive or negative integer is
  175 returned.
  176 -------------------------------------------------------------------------------
  177 */
  178 static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 )
  179 {
int8 roundingMode;
flag roundNearestEven, increment;
int64 z;
roundingMode = float_rounding_mode();
roundNearestEven = ( roundingMode == float_round_nearest_even );
increment = ( (sbits64) absZ1 < 0 );
  187     if ( ! roundNearestEven ) {
  188         if ( roundingMode == float_round_to_zero ) {
increment = 0;
  190         }
  191         else {
  192             if ( zSign ) {
increment = ( roundingMode == float_round_down ) && absZ1;
  194             }
  195             else {
increment = ( roundingMode == float_round_up ) && absZ1;
  197             }
  198         }
  199     }
  200     if ( increment ) {
  201         ++absZ0;
  202         if ( absZ0 == 0 ) goto overflow;
absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );
  204     }
z = absZ0;
  206     if ( zSign ) z = - z;
  207     if ( z && ( ( z < 0 ) ^ zSign ) ) {
overflow:
float_raise( float_flag_invalid );
  210         return
zSign ? (sbits64) LIT64( 0x8000000000000000 )
  212             : LIT64( 0x7FFFFFFFFFFFFFFF );
  213     }
  214     if ( absZ1 ) float_set_inexact();
  215     return z;
  217 }
  218 #endif
  220 /*
  221 -------------------------------------------------------------------------------
  222 Returns the fraction bits of the single-precision floating-point value `a'.
  223 -------------------------------------------------------------------------------
  224 */
INLINE bits32 extractFloat32Frac( float32 a )
  226 {
  228     return a & 0x007FFFFF;
  230 }
  232 /*
  233 -------------------------------------------------------------------------------
  234 Returns the exponent bits of the single-precision floating-point value `a'.
  235 -------------------------------------------------------------------------------
  236 */
INLINE int16 extractFloat32Exp( float32 a )
  238 {
  240     return ( a>>23 ) & 0xFF;
  242 }
  244 /*
  245 -------------------------------------------------------------------------------
  246 Returns the sign bit of the single-precision floating-point value `a'.
  247 -------------------------------------------------------------------------------
  248 */
INLINE flag extractFloat32Sign( float32 a )
  250 {
  252     return a>>31;
  254 }
  256 /*
  257 -------------------------------------------------------------------------------
  258 Normalizes the subnormal single-precision floating-point value represented
  259 by the denormalized significand `aSig'.  The normalized exponent and
  260 significand are stored at the locations pointed to by `zExpPtr' and
  261 `zSigPtr', respectively.
  262 -------------------------------------------------------------------------------
  263 */
  264 static void
normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
  266 {
int8 shiftCount;
shiftCount = countLeadingZeros32( aSig ) - 8;
  270     *zSigPtr = aSig<<shiftCount;
  271     *zExpPtr = 1 - shiftCount;
  273 }
  275 /*
  276 -------------------------------------------------------------------------------
  277 Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
  278 single-precision floating-point value, returning the result.  After being
  279 shifted into the proper positions, the three fields are simply added
  280 together to form the result.  This means that any integer portion of `zSig'
  281 will be added into the exponent.  Since a properly normalized significand
  282 will have an integer portion equal to 1, the `zExp' input should be 1 less
  283 than the desired result exponent whenever `zSig' is a complete, normalized
  284 significand.
  285 -------------------------------------------------------------------------------
  286 */
INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
  288 {
  290     return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
  292 }
  294 /*
  295 -------------------------------------------------------------------------------
  296 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
  297 and significand `zSig', and returns the proper single-precision floating-
  298 point value corresponding to the abstract input.  Ordinarily, the abstract
  299 value is simply rounded and packed into the single-precision format, with
  300 the inexact exception raised if the abstract input cannot be represented
  301 exactly.  However, if the abstract value is too large, the overflow and
  302 inexact exceptions are raised and an infinity or maximal finite value is
  303 returned.  If the abstract value is too small, the input value is rounded to
  304 a subnormal number, and the underflow and inexact exceptions are raised if
  305 the abstract input cannot be represented exactly as a subnormal single-
  306 precision floating-point number.
  307     The input significand `zSig' has its binary point between bits 30
  308 and 29, which is 7 bits to the left of the usual location.  This shifted
  309 significand must be normalized or smaller.  If `zSig' is not normalized,
  310 `zExp' must be 0; in that case, the result returned is a subnormal number,
  311 and it must not require rounding.  In the usual case that `zSig' is
  312 normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
  313 The handling of underflow and overflow follows the IEC/IEEE Standard for
  314 Binary Floating-Point Arithmetic.
  315 -------------------------------------------------------------------------------
  316 */
  317 static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig )
  318 {
int8 roundingMode;
flag roundNearestEven;
int8 roundIncrement, roundBits;
flag isTiny;
roundingMode = float_rounding_mode();
roundNearestEven = ( roundingMode == float_round_nearest_even );
roundIncrement = 0x40;
  327     if ( ! roundNearestEven ) {
  328         if ( roundingMode == float_round_to_zero ) {
roundIncrement = 0;
  330         }
  331         else {
roundIncrement = 0x7F;
  333             if ( zSign ) {
  334                 if ( roundingMode == float_round_up ) roundIncrement = 0;
  335             }
  336             else {
  337                 if ( roundingMode == float_round_down ) roundIncrement = 0;
  338             }
  339         }
  340     }
roundBits = zSig & 0x7F;
  342     if ( 0xFD <= (bits16) zExp ) {
  343         if (    ( 0xFD < zExp )
  344              || (    ( zExp == 0xFD )
  345                   && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
  346            ) {
float_raise( float_flag_overflow | float_flag_inexact );
  348             return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
  349         }
  350         if ( zExp < 0 ) {
isTiny =
  352                    ( float_detect_tininess == float_tininess_before_rounding )
  353                 || ( zExp < -1 )
  354                 || ( zSig + roundIncrement < 0x80000000 );
shift32RightJamming( zSig, - zExp, &zSig );
zExp = 0;
roundBits = zSig & 0x7F;
  358             if ( isTiny && roundBits ) float_raise( float_flag_underflow );
  359         }
  360     }
  361     if ( roundBits ) float_set_inexact();
zSig = ( zSig + roundIncrement )>>7;
zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
  364     if ( zSig == 0 ) zExp = 0;
  365     return packFloat32( zSign, zExp, zSig );
  367 }
  369 /*
  370 -------------------------------------------------------------------------------
  371 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
  372 and significand `zSig', and returns the proper single-precision floating-
  373 point value corresponding to the abstract input.  This routine is just like
  374 `roundAndPackFloat32' except that `zSig' does not have to be normalized.
  375 Bit 31 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
  376 floating-point exponent.
  377 -------------------------------------------------------------------------------
  378 */
  379 static float32
normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig )
  381 {
int8 shiftCount;
shiftCount = countLeadingZeros32( zSig ) - 1;
  385     return roundAndPackFloat32( zSign, zExp - shiftCount, zSig<<shiftCount );
  387 }
  389 /*
  390 -------------------------------------------------------------------------------
  391 Returns the fraction bits of the double-precision floating-point value `a'.
  392 -------------------------------------------------------------------------------
  393 */
INLINE bits64 extractFloat64Frac( float64 a )
  395 {
  397     return FLOAT64_DEMANGLE(a) & LIT64( 0x000FFFFFFFFFFFFF );
  399 }
  401 /*
  402 -------------------------------------------------------------------------------
  403 Returns the exponent bits of the double-precision floating-point value `a'.
  404 -------------------------------------------------------------------------------
  405 */
INLINE int16 extractFloat64Exp( float64 a )
  407 {
  409     return ( FLOAT64_DEMANGLE(a)>>52 ) & 0x7FF;
  411 }
  413 /*
  414 -------------------------------------------------------------------------------
  415 Returns the sign bit of the double-precision floating-point value `a'.
  416 -------------------------------------------------------------------------------
  417 */
INLINE flag extractFloat64Sign( float64 a )
  419 {
  421     return FLOAT64_DEMANGLE(a)>>63;
  423 }
  425 /*
  426 -------------------------------------------------------------------------------
  427 Normalizes the subnormal double-precision floating-point value represented
  428 by the denormalized significand `aSig'.  The normalized exponent and
  429 significand are stored at the locations pointed to by `zExpPtr' and
  430 `zSigPtr', respectively.
  431 -------------------------------------------------------------------------------
  432 */
  433 static void
normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
  435 {
int8 shiftCount;
shiftCount = countLeadingZeros64( aSig ) - 11;
  439     *zSigPtr = aSig<<shiftCount;
  440     *zExpPtr = 1 - shiftCount;
  442 }
  444 /*
  445 -------------------------------------------------------------------------------
  446 Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
  447 double-precision floating-point value, returning the result.  After being
  448 shifted into the proper positions, the three fields are simply added
  449 together to form the result.  This means that any integer portion of `zSig'
  450 will be added into the exponent.  Since a properly normalized significand
  451 will have an integer portion equal to 1, the `zExp' input should be 1 less
  452 than the desired result exponent whenever `zSig' is a complete, normalized
  453 significand.
  454 -------------------------------------------------------------------------------
  455 */
INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
  457 {
  459     return FLOAT64_MANGLE( ( ( (bits64) zSign )<<63 ) +
  460                            ( ( (bits64) zExp )<<52 ) + zSig );
  462 }
  464 /*
  465 -------------------------------------------------------------------------------
  466 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
  467 and significand `zSig', and returns the proper double-precision floating-
  468 point value corresponding to the abstract input.  Ordinarily, the abstract
  469 value is simply rounded and packed into the double-precision format, with
  470 the inexact exception raised if the abstract input cannot be represented
  471 exactly.  However, if the abstract value is too large, the overflow and
  472 inexact exceptions are raised and an infinity or maximal finite value is
  473 returned.  If the abstract value is too small, the input value is rounded to
  474 a subnormal number, and the underflow and inexact exceptions are raised if
  475 the abstract input cannot be represented exactly as a subnormal double-
  476 precision floating-point number.
  477     The input significand `zSig' has its binary point between bits 62
  478 and 61, which is 10 bits to the left of the usual location.  This shifted
  479 significand must be normalized or smaller.  If `zSig' is not normalized,
  480 `zExp' must be 0; in that case, the result returned is a subnormal number,
  481 and it must not require rounding.  In the usual case that `zSig' is
  482 normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
  483 The handling of underflow and overflow follows the IEC/IEEE Standard for
  484 Binary Floating-Point Arithmetic.
  485 -------------------------------------------------------------------------------
  486 */
  487 static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig )
  488 {
int8 roundingMode;
flag roundNearestEven;
int16 roundIncrement, roundBits;
flag isTiny;
roundingMode = float_rounding_mode();
roundNearestEven = ( roundingMode == float_round_nearest_even );
roundIncrement = 0x200;
  497     if ( ! roundNearestEven ) {
  498         if ( roundingMode == float_round_to_zero ) {
roundIncrement = 0;
  500         }
  501         else {
roundIncrement = 0x3FF;
  503             if ( zSign ) {
  504                 if ( roundingMode == float_round_up ) roundIncrement = 0;
  505             }
  506             else {
  507                 if ( roundingMode == float_round_down ) roundIncrement = 0;
  508             }
  509         }
  510     }
roundBits = zSig & 0x3FF;
  512     if ( 0x7FD <= (bits16) zExp ) {
  513         if (    ( 0x7FD < zExp )
  514              || (    ( zExp == 0x7FD )
  515                   && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
  516            ) {
float_raise( float_flag_overflow | float_flag_inexact );
  518             return FLOAT64_MANGLE(
FLOAT64_DEMANGLE(packFloat64( zSign, 0x7FF, 0 )) -
  520                 ( roundIncrement == 0 ));
  521         }
  522         if ( zExp < 0 ) {
isTiny =
  524                    ( float_detect_tininess == float_tininess_before_rounding )
  525                 || ( zExp < -1 )
  526                 || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
shift64RightJamming( zSig, - zExp, &zSig );
zExp = 0;
roundBits = zSig & 0x3FF;
  530             if ( isTiny && roundBits ) float_raise( float_flag_underflow );
  531         }
  532     }
  533     if ( roundBits ) float_set_inexact();
zSig = ( zSig + roundIncrement )>>10;
zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
  536     if ( zSig == 0 ) zExp = 0;
  537     return packFloat64( zSign, zExp, zSig );
  539 }
  541 /*
  542 -------------------------------------------------------------------------------
  543 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
  544 and significand `zSig', and returns the proper double-precision floating-
  545 point value corresponding to the abstract input.  This routine is just like
  546 `roundAndPackFloat64' except that `zSig' does not have to be normalized.
  547 Bit 63 of `zSig' must be zero, and `zExp' must be 1 less than the ``true''
  548 floating-point exponent.
  549 -------------------------------------------------------------------------------
  550 */
  551 static float64
normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig )
  553 {
int8 shiftCount;
shiftCount = countLeadingZeros64( zSig ) - 1;
  557     return roundAndPackFloat64( zSign, zExp - shiftCount, zSig<<shiftCount );
  559 }
  561 #ifdef FLOATX80
  563 /*
  564 -------------------------------------------------------------------------------
  565 Returns the fraction bits of the extended double-precision floating-point
  566 value `a'.
  567 -------------------------------------------------------------------------------
  568 */
INLINE bits64 extractFloatx80Frac( floatx80 a )
  570 {
  572     return a.low;
  574 }
  576 /*
  577 -------------------------------------------------------------------------------
  578 Returns the exponent bits of the extended double-precision floating-point
  579 value `a'.
  580 -------------------------------------------------------------------------------
  581 */
INLINE int32 extractFloatx80Exp( floatx80 a )
  583 {
  585     return a.high & 0x7FFF;
  587 }
  589 /*
  590 -------------------------------------------------------------------------------
  591 Returns the sign bit of the extended double-precision floating-point value
  592 `a'.
  593 -------------------------------------------------------------------------------
  594 */
INLINE flag extractFloatx80Sign( floatx80 a )
  596 {
  598     return a.high>>15;
  600 }
  602 /*
  603 -------------------------------------------------------------------------------
  604 Normalizes the subnormal extended double-precision floating-point value
  605 represented by the denormalized significand `aSig'.  The normalized exponent
  606 and significand are stored at the locations pointed to by `zExpPtr' and
  607 `zSigPtr', respectively.
  608 -------------------------------------------------------------------------------
  609 */
  610 static void
normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
  612 {
int8 shiftCount;
shiftCount = countLeadingZeros64( aSig );
  616     *zSigPtr = aSig<<shiftCount;
  617     *zExpPtr = 1 - shiftCount;
  619 }
  621 /*
  622 -------------------------------------------------------------------------------
  623 Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
  624 extended double-precision floating-point value, returning the result.
  625 -------------------------------------------------------------------------------
  626 */
INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
  628 {
floatx80 z;
z.low = zSig;
z.high = ( ( (bits16) zSign )<<15 ) + zExp;
  633     return z;
  635 }
  637 /*
  638 -------------------------------------------------------------------------------
  639 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
  640 and extended significand formed by the concatenation of `zSig0' and `zSig1',
  641 and returns the proper extended double-precision floating-point value
  642 corresponding to the abstract input.  Ordinarily, the abstract value is
  643 rounded and packed into the extended double-precision format, with the
  644 inexact exception raised if the abstract input cannot be represented
  645 exactly.  However, if the abstract value is too large, the overflow and
  646 inexact exceptions are raised and an infinity or maximal finite value is
  647 returned.  If the abstract value is too small, the input value is rounded to
  648 a subnormal number, and the underflow and inexact exceptions are raised if
  649 the abstract input cannot be represented exactly as a subnormal extended
  650 double-precision floating-point number.
  651     If `roundingPrecision' is 32 or 64, the result is rounded to the same
  652 number of bits as single or double precision, respectively.  Otherwise, the
  653 result is rounded to the full precision of the extended double-precision
  654 format.
  655     The input significand must be normalized or smaller.  If the input
  656 significand is not normalized, `zExp' must be 0; in that case, the result
  657 returned is a subnormal number, and it must not require rounding.  The
  658 handling of underflow and overflow follows the IEC/IEEE Standard for Binary
  659 Floating-Point Arithmetic.
  660 -------------------------------------------------------------------------------
  661 */
  662 static floatx80
roundAndPackFloatx80(
int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
  665  )
  666 {
int8 roundingMode;
flag roundNearestEven, increment, isTiny;
int64 roundIncrement, roundMask, roundBits;
roundingMode = float_rounding_mode();
roundNearestEven = ( roundingMode == float_round_nearest_even );
  673     if ( roundingPrecision == 80 ) goto precision80;
  674     if ( roundingPrecision == 64 ) {
roundIncrement = LIT64( 0x0000000000000400 );
roundMask = LIT64( 0x00000000000007FF );
  677     }
  678     else if ( roundingPrecision == 32 ) {
roundIncrement = LIT64( 0x0000008000000000 );
roundMask = LIT64( 0x000000FFFFFFFFFF );
  681     }
  682     else {
  683         goto precision80;
  684     }
zSig0 |= ( zSig1 != 0 );
  686     if ( ! roundNearestEven ) {
  687         if ( roundingMode == float_round_to_zero ) {
roundIncrement = 0;
  689         }
  690         else {
roundIncrement = roundMask;
  692             if ( zSign ) {
  693                 if ( roundingMode == float_round_up ) roundIncrement = 0;
  694             }
  695             else {
  696                 if ( roundingMode == float_round_down ) roundIncrement = 0;
  697             }
  698         }
  699     }
roundBits = zSig0 & roundMask;
  701     if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
  702         if (    ( 0x7FFE < zExp )
  703              || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
  704            ) {
  705             goto overflow;
  706         }
  707         if ( zExp <= 0 ) {
isTiny =
  709                    ( float_detect_tininess == float_tininess_before_rounding )
  710                 || ( zExp < 0 )
  711                 || ( zSig0 <= zSig0 + roundIncrement );
shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
zExp = 0;
roundBits = zSig0 & roundMask;
  715             if ( isTiny && roundBits ) float_raise( float_flag_underflow );
  716             if ( roundBits ) float_set_inexact();
zSig0 += roundIncrement;
  718             if ( (sbits64) zSig0 < 0 ) zExp = 1;
roundIncrement = roundMask + 1;
  720             if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
roundMask |= roundIncrement;
  722             }
zSig0 &= ~ roundMask;
  724             return packFloatx80( zSign, zExp, zSig0 );
  725         }
  726     }
  727     if ( roundBits ) float_set_inexact();
zSig0 += roundIncrement;
  729     if ( zSig0 < roundIncrement ) {
  730         ++zExp;
zSig0 = LIT64( 0x8000000000000000 );
  732     }
roundIncrement = roundMask + 1;
  734     if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
roundMask |= roundIncrement;
  736     }
zSig0 &= ~ roundMask;
  738     if ( zSig0 == 0 ) zExp = 0;
  739     return packFloatx80( zSign, zExp, zSig0 );
precision80:
increment = ( (sbits64) zSig1 < 0 );
  742     if ( ! roundNearestEven ) {
  743         if ( roundingMode == float_round_to_zero ) {
increment = 0;
  745         }
  746         else {
  747             if ( zSign ) {
increment = ( roundingMode == float_round_down ) && zSig1;
  749             }
  750             else {
increment = ( roundingMode == float_round_up ) && zSig1;
  752             }
  753         }
  754     }
  755     if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
  756         if (    ( 0x7FFE < zExp )
  757              || (    ( zExp == 0x7FFE )
  758                   && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
  759                   && increment
  760                 )
  761            ) {
roundMask = 0;
overflow:
float_raise( float_flag_overflow | float_flag_inexact );
  765             if (    ( roundingMode == float_round_to_zero )
  766                  || ( zSign && ( roundingMode == float_round_up ) )
  767                  || ( ! zSign && ( roundingMode == float_round_down ) )
  768                ) {
  769                 return packFloatx80( zSign, 0x7FFE, ~ roundMask );
  770             }
  771             return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
  772         }
  773         if ( zExp <= 0 ) {
isTiny =
  775                    ( float_detect_tininess == float_tininess_before_rounding )
  776                 || ( zExp < 0 )
  777                 || ! increment
  778                 || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
zExp = 0;
  781             if ( isTiny && zSig1 ) float_raise( float_flag_underflow );
  782             if ( zSig1 ) float_set_inexact();
  783             if ( roundNearestEven ) {
increment = ( (sbits64) zSig1 < 0 );
  785             }
  786             else {
  787                 if ( zSign ) {
increment = ( roundingMode == float_round_down ) && zSig1;
  789                 }
  790                 else {
increment = ( roundingMode == float_round_up ) && zSig1;
  792                 }
  793             }
  794             if ( increment ) {
  795                 ++zSig0;
zSig0 &=
  797                     ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
  798                 if ( (sbits64) zSig0 < 0 ) zExp = 1;
  799             }
  800             return packFloatx80( zSign, zExp, zSig0 );
  801         }
  802     }
  803     if ( zSig1 ) float_set_inexact();
  804     if ( increment ) {
  805         ++zSig0;
  806         if ( zSig0 == 0 ) {
  807             ++zExp;
zSig0 = LIT64( 0x8000000000000000 );
  809         }
  810         else {
zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
  812         }
  813     }
  814     else {
  815         if ( zSig0 == 0 ) zExp = 0;
  816     }
  817     return packFloatx80( zSign, zExp, zSig0 );
  819 }
  821 /*
  822 -------------------------------------------------------------------------------
  823 Takes an abstract floating-point value having sign `zSign', exponent
  824 `zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
  825 and returns the proper extended double-precision floating-point value
  826 corresponding to the abstract input.  This routine is just like
  827 `roundAndPackFloatx80' except that the input significand does not have to be
  828 normalized.
  829 -------------------------------------------------------------------------------
  830 */
  831 static floatx80
normalizeRoundAndPackFloatx80(
int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
  834  )
  835 {
int8 shiftCount;
  838     if ( zSig0 == 0 ) {
zSig0 = zSig1;
zSig1 = 0;
zExp -= 64;
  842     }
shiftCount = countLeadingZeros64( zSig0 );
shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
zExp -= shiftCount;
  846     return
roundAndPackFloatx80( roundingPrecision, zSign, zExp, zSig0, zSig1 );
  849 }
  851 #endif
  853 #ifdef FLOAT128
  855 /*
  856 -------------------------------------------------------------------------------
  857 Returns the least-significant 64 fraction bits of the quadruple-precision
  858 floating-point value `a'.
  859 -------------------------------------------------------------------------------
  860 */
INLINE bits64 extractFloat128Frac1( float128 a )
  862 {
  864     return a.low;
  866 }
  868 /*
  869 -------------------------------------------------------------------------------
  870 Returns the most-significant 48 fraction bits of the quadruple-precision
  871 floating-point value `a'.
  872 -------------------------------------------------------------------------------
  873 */
INLINE bits64 extractFloat128Frac0( float128 a )
  875 {
  877     return a.high & LIT64( 0x0000FFFFFFFFFFFF );
  879 }
  881 /*
  882 -------------------------------------------------------------------------------
  883 Returns the exponent bits of the quadruple-precision floating-point value
  884 `a'.
  885 -------------------------------------------------------------------------------
  886 */
INLINE int32 extractFloat128Exp( float128 a )
  888 {
  890     return ( a.high>>48 ) & 0x7FFF;
  892 }
  894 /*
  895 -------------------------------------------------------------------------------
  896 Returns the sign bit of the quadruple-precision floating-point value `a'.
  897 -------------------------------------------------------------------------------
  898 */
INLINE flag extractFloat128Sign( float128 a )
  900 {
  902     return a.high>>63;
  904 }
  906 /*
  907 -------------------------------------------------------------------------------
  908 Normalizes the subnormal quadruple-precision floating-point value
  909 represented by the denormalized significand formed by the concatenation of
  910 `aSig0' and `aSig1'.  The normalized exponent is stored at the location
  911 pointed to by `zExpPtr'.  The most significant 49 bits of the normalized
  912 significand are stored at the location pointed to by `zSig0Ptr', and the
  913 least significant 64 bits of the normalized significand are stored at the
  914 location pointed to by `zSig1Ptr'.
  915 -------------------------------------------------------------------------------
  916 */
  917 static void
normalizeFloat128Subnormal(
bits64 aSig0,
bits64 aSig1,
int32 *zExpPtr,
bits64 *zSig0Ptr,
bits64 *zSig1Ptr
  924  )
  925 {
int8 shiftCount;
  928     if ( aSig0 == 0 ) {
shiftCount = countLeadingZeros64( aSig1 ) - 15;
  930         if ( shiftCount < 0 ) {
  931             *zSig0Ptr = aSig1>>( - shiftCount );
  932             *zSig1Ptr = aSig1<<( shiftCount & 63 );
  933         }
  934         else {
  935             *zSig0Ptr = aSig1<<shiftCount;
  936             *zSig1Ptr = 0;
  937         }
  938         *zExpPtr = - shiftCount - 63;
  939     }
  940     else {
shiftCount = countLeadingZeros64( aSig0 ) - 15;
shortShift128Left( aSig0, aSig1, shiftCount, zSig0Ptr, zSig1Ptr );
  943         *zExpPtr = 1 - shiftCount;
  944     }
  946 }
  948 /*
  949 -------------------------------------------------------------------------------
  950 Packs the sign `zSign', the exponent `zExp', and the significand formed
  951 by the concatenation of `zSig0' and `zSig1' into a quadruple-precision
  952 floating-point value, returning the result.  After being shifted into the
  953 proper positions, the three fields `zSign', `zExp', and `zSig0' are simply
  954 added together to form the most significant 32 bits of the result.  This
  955 means that any integer portion of `zSig0' will be added into the exponent.
  956 Since a properly normalized significand will have an integer portion equal
  957 to 1, the `zExp' input should be 1 less than the desired result exponent
  958 whenever `zSig0' and `zSig1' concatenated form a complete, normalized
  959 significand.
  960 -------------------------------------------------------------------------------
  961 */
INLINE float128
packFloat128( flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
  964 {
float128 z;
z.low = zSig1;
z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;
  969     return z;
  971 }
  973 /*
  974 -------------------------------------------------------------------------------
  975 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
  976 and extended significand formed by the concatenation of `zSig0', `zSig1',
  977 and `zSig2', and returns the proper quadruple-precision floating-point value
  978 corresponding to the abstract input.  Ordinarily, the abstract value is
  979 simply rounded and packed into the quadruple-precision format, with the
  980 inexact exception raised if the abstract input cannot be represented
  981 exactly.  However, if the abstract value is too large, the overflow and
  982 inexact exceptions are raised and an infinity or maximal finite value is
  983 returned.  If the abstract value is too small, the input value is rounded to
  984 a subnormal number, and the underflow and inexact exceptions are raised if
  985 the abstract input cannot be represented exactly as a subnormal quadruple-
  986 precision floating-point number.
  987     The input significand must be normalized or smaller.  If the input
  988 significand is not normalized, `zExp' must be 0; in that case, the result
  989 returned is a subnormal number, and it must not require rounding.  In the
  990 usual case that the input significand is normalized, `zExp' must be 1 less
  991 than the ``true'' floating-point exponent.  The handling of underflow and
  992 overflow follows the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
  993 -------------------------------------------------------------------------------
  994 */
  995 static float128
roundAndPackFloat128(
flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 )
  998 {
int8 roundingMode;
flag roundNearestEven, increment, isTiny;
roundingMode = float_rounding_mode();
roundNearestEven = ( roundingMode == float_round_nearest_even );
increment = ( (sbits64) zSig2 < 0 );
 1005     if ( ! roundNearestEven ) {
 1006         if ( roundingMode == float_round_to_zero ) {
increment = 0;
 1008         }
 1009         else {
 1010             if ( zSign ) {
increment = ( roundingMode == float_round_down ) && zSig2;
 1012             }
 1013             else {
increment = ( roundingMode == float_round_up ) && zSig2;
 1015             }
 1016         }
 1017     }
 1018     if ( 0x7FFD <= (bits32) zExp ) {
 1019         if (    ( 0x7FFD < zExp )
 1020              || (    ( zExp == 0x7FFD )
 1021                   && eq128(
LIT64( 0x0001FFFFFFFFFFFF ),
LIT64( 0xFFFFFFFFFFFFFFFF ),
zSig0,
zSig1
 1026                      )
 1027                   && increment
 1028                 )
 1029            ) {
float_raise( float_flag_overflow | float_flag_inexact );
 1031             if (    ( roundingMode == float_round_to_zero )
 1032                  || ( zSign && ( roundingMode == float_round_up ) )
 1033                  || ( ! zSign && ( roundingMode == float_round_down ) )
 1034                ) {
 1035                 return
packFloat128(
zSign,
 1038                         0x7FFE,
LIT64( 0x0000FFFFFFFFFFFF ),
LIT64( 0xFFFFFFFFFFFFFFFF )
 1041                     );
 1042             }
 1043             return packFloat128( zSign, 0x7FFF, 0, 0 );
 1044         }
 1045         if ( zExp < 0 ) {
isTiny =
 1047                    ( float_detect_tininess == float_tininess_before_rounding )
 1048                 || ( zExp < -1 )
 1049                 || ! increment
 1050                 || lt128(
zSig0,
zSig1,
LIT64( 0x0001FFFFFFFFFFFF ),
LIT64( 0xFFFFFFFFFFFFFFFF )
 1055                    );
shift128ExtraRightJamming(
zSig0, zSig1, zSig2, - zExp, &zSig0, &zSig1, &zSig2 );
zExp = 0;
 1059             if ( isTiny && zSig2 ) float_raise( float_flag_underflow );
 1060             if ( roundNearestEven ) {
increment = ( (sbits64) zSig2 < 0 );
 1062             }
 1063             else {
 1064                 if ( zSign ) {
increment = ( roundingMode == float_round_down ) && zSig2;
 1066                 }
 1067                 else {
increment = ( roundingMode == float_round_up ) && zSig2;
 1069                 }
 1070             }
 1071         }
 1072     }
 1073     if ( zSig2 ) float_set_inexact();
 1074     if ( increment ) {
add128( zSig0, zSig1, 0, 1, &zSig0, &zSig1 );
zSig1 &= ~ ( ( zSig2 + zSig2 == 0 ) & roundNearestEven );
 1077     }
 1078     else {
 1079         if ( ( zSig0 | zSig1 ) == 0 ) zExp = 0;
 1080     }
 1081     return packFloat128( zSign, zExp, zSig0, zSig1 );
 1083 }
 1085 /*
 1086 -------------------------------------------------------------------------------
 1087 Takes an abstract floating-point value having sign `zSign', exponent `zExp',
 1088 and significand formed by the concatenation of `zSig0' and `zSig1', and
 1089 returns the proper quadruple-precision floating-point value corresponding
 1090 to the abstract input.  This routine is just like `roundAndPackFloat128'
 1091 except that the input significand has fewer bits and does not have to be
 1092 normalized.  In all cases, `zExp' must be 1 less than the ``true'' floating-
 1093 point exponent.
 1094 -------------------------------------------------------------------------------
 1095 */
 1096 static float128
normalizeRoundAndPackFloat128(
flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
 1099 {
int8 shiftCount;
bits64 zSig2;
 1103     if ( zSig0 == 0 ) {
zSig0 = zSig1;
zSig1 = 0;
zExp -= 64;
 1107     }
shiftCount = countLeadingZeros64( zSig0 ) - 15;
 1109     if ( 0 <= shiftCount ) {
zSig2 = 0;
shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
 1112     }
 1113     else {
shift128ExtraRightJamming(
zSig0, zSig1, 0, - shiftCount, &zSig0, &zSig1, &zSig2 );
 1116     }
zExp -= shiftCount;
 1118     return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
 1120 }
 1122 #endif
 1124 /*
 1125 -------------------------------------------------------------------------------
 1126 Returns the result of converting the 32-bit two's complement integer `a'
 1127 to the single-precision floating-point format.  The conversion is performed
 1128 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1129 -------------------------------------------------------------------------------
 1130 */
float32 int32_to_float32( int32 a )
 1132 {
flag zSign;
 1135     if ( a == 0 ) return 0;
 1136     if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
zSign = ( a < 0 );
 1138     return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a );
 1140 }
 1142 /*
 1143 -------------------------------------------------------------------------------
 1144 Returns the result of converting the 32-bit two's complement integer `a'
 1145 to the double-precision floating-point format.  The conversion is performed
 1146 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1147 -------------------------------------------------------------------------------
 1148 */
float64 int32_to_float64( int32 a )
 1150 {
flag zSign;
uint32 absA;
int8 shiftCount;
bits64 zSig;
 1156     if ( a == 0 ) return 0;
zSign = ( a < 0 );
absA = zSign ? - a : a;
shiftCount = countLeadingZeros32( absA ) + 21;
zSig = absA;
 1161     return packFloat64( zSign, 0x432 - shiftCount, zSig<<shiftCount );
 1163 }
 1165 #ifdef FLOATX80
 1167 /*
 1168 -------------------------------------------------------------------------------
 1169 Returns the result of converting the 32-bit two's complement integer `a'
 1170 to the extended double-precision floating-point format.  The conversion
 1171 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 1172 Arithmetic.
 1173 -------------------------------------------------------------------------------
 1174 */
floatx80 int32_to_floatx80( int32 a )
 1176 {
flag zSign;
uint32 absA;
int8 shiftCount;
bits64 zSig;
 1182     if ( a == 0 ) return packFloatx80( 0, 0, 0 );
zSign = ( a < 0 );
absA = zSign ? - a : a;
shiftCount = countLeadingZeros32( absA ) + 32;
zSig = absA;
 1187     return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
 1189 }
 1191 #endif
 1193 #ifdef FLOAT128
 1195 /*
 1196 -------------------------------------------------------------------------------
 1197 Returns the result of converting the 32-bit two's complement integer `a' to
 1198 the quadruple-precision floating-point format.  The conversion is performed
 1199 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1200 -------------------------------------------------------------------------------
 1201 */
float128 int32_to_float128( int32 a )
 1203 {
flag zSign;
uint32 absA;
int8 shiftCount;
bits64 zSig0;
 1209     if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
zSign = ( a < 0 );
absA = zSign ? - a : a;
shiftCount = countLeadingZeros32( absA ) + 17;
zSig0 = absA;
 1214     return packFloat128( zSign, 0x402E - shiftCount, zSig0<<shiftCount, 0 );
 1216 }
 1218 #endif
 1220 #ifndef SOFTFLOAT_FOR_GCC /* __floatdi?f is in libgcc2.c */
 1221 /*
 1222 -------------------------------------------------------------------------------
 1223 Returns the result of converting the 64-bit two's complement integer `a'
 1224 to the single-precision floating-point format.  The conversion is performed
 1225 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1226 -------------------------------------------------------------------------------
 1227 */
float32 int64_to_float32( int64 a )
 1229 {
flag zSign;
uint64 absA;
int8 shiftCount;
 1234     if ( a == 0 ) return 0;
zSign = ( a < 0 );
absA = zSign ? - a : a;
shiftCount = countLeadingZeros64( absA ) - 40;
 1238     if ( 0 <= shiftCount ) {
 1239         return packFloat32( zSign, 0x95 - shiftCount, absA<<shiftCount );
 1240     }
 1241     else {
shiftCount += 7;
 1243         if ( shiftCount < 0 ) {
shift64RightJamming( absA, - shiftCount, &absA );
 1245         }
 1246         else {
absA <<= shiftCount;
 1248         }
 1249         return roundAndPackFloat32( zSign, 0x9C - shiftCount, absA );
 1250     }
 1252 }
 1254 /*
 1255 -------------------------------------------------------------------------------
 1256 Returns the result of converting the 64-bit two's complement integer `a'
 1257 to the double-precision floating-point format.  The conversion is performed
 1258 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1259 -------------------------------------------------------------------------------
 1260 */
float64 int64_to_float64( int64 a )
 1262 {
flag zSign;
 1265     if ( a == 0 ) return 0;
 1266     if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
 1267         return packFloat64( 1, 0x43E, 0 );
 1268     }
zSign = ( a < 0 );
 1270     return normalizeRoundAndPackFloat64( zSign, 0x43C, zSign ? - a : a );
 1272 }
 1274 #ifdef FLOATX80
 1276 /*
 1277 -------------------------------------------------------------------------------
 1278 Returns the result of converting the 64-bit two's complement integer `a'
 1279 to the extended double-precision floating-point format.  The conversion
 1280 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 1281 Arithmetic.
 1282 -------------------------------------------------------------------------------
 1283 */
floatx80 int64_to_floatx80( int64 a )
 1285 {
flag zSign;
uint64 absA;
int8 shiftCount;
 1290     if ( a == 0 ) return packFloatx80( 0, 0, 0 );
zSign = ( a < 0 );
absA = zSign ? - a : a;
shiftCount = countLeadingZeros64( absA );
 1294     return packFloatx80( zSign, 0x403E - shiftCount, absA<<shiftCount );
 1296 }
 1298 #endif
 1300 #ifdef FLOAT128
 1302 /*
 1303 -------------------------------------------------------------------------------
 1304 Returns the result of converting the 64-bit two's complement integer `a' to
 1305 the quadruple-precision floating-point format.  The conversion is performed
 1306 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1307 -------------------------------------------------------------------------------
 1308 */
float128 int64_to_float128( int64 a )
 1310 {
flag zSign;
uint64 absA;
int8 shiftCount;
int32 zExp;
bits64 zSig0, zSig1;
 1317     if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
zSign = ( a < 0 );
absA = zSign ? - a : a;
shiftCount = countLeadingZeros64( absA ) + 49;
zExp = 0x406E - shiftCount;
 1322     if ( 64 <= shiftCount ) {
zSig1 = 0;
zSig0 = absA;
shiftCount -= 64;
 1326     }
 1327     else {
zSig1 = absA;
zSig0 = 0;
 1330     }
shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
 1332     return packFloat128( zSign, zExp, zSig0, zSig1 );
 1334 }
 1336 #endif
 1337 #endif /* !SOFTFLOAT_FOR_GCC */
 1339 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 1340 /*
 1341 -------------------------------------------------------------------------------
 1342 Returns the result of converting the single-precision floating-point value
 1343 `a' to the 32-bit two's complement integer format.  The conversion is
 1344 performed according to the IEC/IEEE Standard for Binary Floating-Point
 1345 Arithmetic---which means in particular that the conversion is rounded
 1346 according to the current rounding mode.  If `a' is a NaN, the largest
 1347 positive integer is returned.  Otherwise, if the conversion overflows, the
 1348 largest integer with the same sign as `a' is returned.
 1349 -------------------------------------------------------------------------------
 1350 */
int32 float32_to_int32( float32 a )
 1352 {
flag aSign;
int16 aExp, shiftCount;
bits32 aSig;
bits64 aSig64;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
 1361     if ( ( aExp == 0xFF ) && aSig ) aSign = 0;
 1362     if ( aExp ) aSig |= 0x00800000;
shiftCount = 0xAF - aExp;
aSig64 = aSig;
aSig64 <<= 32;
 1366     if ( 0 < shiftCount ) shift64RightJamming( aSig64, shiftCount, &aSig64 );
 1367     return roundAndPackInt32( aSign, aSig64 );
 1369 }
 1370 #endif /* !SOFTFLOAT_FOR_GCC */
 1372 /*
 1373 -------------------------------------------------------------------------------
 1374 Returns the result of converting the single-precision floating-point value
 1375 `a' to the 32-bit two's complement integer format.  The conversion is
 1376 performed according to the IEC/IEEE Standard for Binary Floating-Point
 1377 Arithmetic, except that the conversion is always rounded toward zero.
 1378 If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
 1379 the conversion overflows, the largest integer with the same sign as `a' is
 1380 returned.
 1381 -------------------------------------------------------------------------------
 1382 */
int32 float32_to_int32_round_to_zero( float32 a )
 1384 {
flag aSign;
int16 aExp, shiftCount;
bits32 aSig;
int32 z;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
shiftCount = aExp - 0x9E;
 1394     if ( 0 <= shiftCount ) {
 1395         if ( a != 0xCF000000 ) {
float_raise( float_flag_invalid );
 1397             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
 1398         }
 1399         return (sbits32) 0x80000000;
 1400     }
 1401     else if ( aExp <= 0x7E ) {
 1402         if ( aExp | aSig ) float_set_inexact();
 1403         return 0;
 1404     }
aSig = ( aSig | 0x00800000 )<<8;
z = aSig>>( - shiftCount );
 1407     if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
float_set_inexact();
 1409     }
 1410     if ( aSign ) z = - z;
 1411     return z;
 1413 }
 1415 #ifndef SOFTFLOAT_FOR_GCC /* __fix?fdi provided by libgcc2.c */
 1416 /*
 1417 -------------------------------------------------------------------------------
 1418 Returns the result of converting the single-precision floating-point value
 1419 `a' to the 64-bit two's complement integer format.  The conversion is
 1420 performed according to the IEC/IEEE Standard for Binary Floating-Point
 1421 Arithmetic---which means in particular that the conversion is rounded
 1422 according to the current rounding mode.  If `a' is a NaN, the largest
 1423 positive integer is returned.  Otherwise, if the conversion overflows, the
 1424 largest integer with the same sign as `a' is returned.
 1425 -------------------------------------------------------------------------------
 1426 */
int64 float32_to_int64( float32 a )
 1428 {
flag aSign;
int16 aExp, shiftCount;
bits32 aSig;
bits64 aSig64, aSigExtra;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
shiftCount = 0xBE - aExp;
 1438     if ( shiftCount < 0 ) {
float_raise( float_flag_invalid );
 1440         if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
 1441             return LIT64( 0x7FFFFFFFFFFFFFFF );
 1442         }
 1443         return (sbits64) LIT64( 0x8000000000000000 );
 1444     }
 1445     if ( aExp ) aSig |= 0x00800000;
aSig64 = aSig;
aSig64 <<= 40;
shift64ExtraRightJamming( aSig64, 0, shiftCount, &aSig64, &aSigExtra );
 1449     return roundAndPackInt64( aSign, aSig64, aSigExtra );
 1451 }
 1453 /*
 1454 -------------------------------------------------------------------------------
 1455 Returns the result of converting the single-precision floating-point value
 1456 `a' to the 64-bit two's complement integer format.  The conversion is
 1457 performed according to the IEC/IEEE Standard for Binary Floating-Point
 1458 Arithmetic, except that the conversion is always rounded toward zero.  If
 1459 `a' is a NaN, the largest positive integer is returned.  Otherwise, if the
 1460 conversion overflows, the largest integer with the same sign as `a' is
 1461 returned.
 1462 -------------------------------------------------------------------------------
 1463 */
int64 float32_to_int64_round_to_zero( float32 a )
 1465 {
flag aSign;
int16 aExp, shiftCount;
bits32 aSig;
bits64 aSig64;
int64 z;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
shiftCount = aExp - 0xBE;
 1476     if ( 0 <= shiftCount ) {
 1477         if ( a != 0xDF000000 ) {
float_raise( float_flag_invalid );
 1479             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
 1480                 return LIT64( 0x7FFFFFFFFFFFFFFF );
 1481             }
 1482         }
 1483         return (sbits64) LIT64( 0x8000000000000000 );
 1484     }
 1485     else if ( aExp <= 0x7E ) {
 1486         if ( aExp | aSig ) float_set_inexact();
 1487         return 0;
 1488     }
aSig64 = aSig | 0x00800000;
aSig64 <<= 40;
z = aSig64>>( - shiftCount );
 1492     if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {
float_set_inexact();
 1494     }
 1495     if ( aSign ) z = - z;
 1496     return z;
 1498 }
 1499 #endif /* !SOFTFLOAT_FOR_GCC */
 1501 /*
 1502 -------------------------------------------------------------------------------
 1503 Returns the result of converting the single-precision floating-point value
 1504 `a' to the double-precision floating-point format.  The conversion is
 1505 performed according to the IEC/IEEE Standard for Binary Floating-Point
 1506 Arithmetic.
 1507 -------------------------------------------------------------------------------
 1508 */
float64 float32_to_float64( float32 a )
 1510 {
flag aSign;
int16 aExp;
bits32 aSig;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
 1518     if ( aExp == 0xFF ) {
 1519         if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) );
 1520         return packFloat64( aSign, 0x7FF, 0 );
 1521     }
 1522     if ( aExp == 0 ) {
 1523         if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 1525         --aExp;
 1526     }
 1527     return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
 1529 }
 1531 #ifdef FLOATX80
 1533 /*
 1534 -------------------------------------------------------------------------------
 1535 Returns the result of converting the single-precision floating-point value
 1536 `a' to the extended double-precision floating-point format.  The conversion
 1537 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 1538 Arithmetic.
 1539 -------------------------------------------------------------------------------
 1540 */
floatx80 float32_to_floatx80( float32 a )
 1542 {
flag aSign;
int16 aExp;
bits32 aSig;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
 1550     if ( aExp == 0xFF ) {
 1551         if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a ) );
 1552         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 1553     }
 1554     if ( aExp == 0 ) {
 1555         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 1557     }
aSig |= 0x00800000;
 1559     return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
 1561 }
 1563 #endif
 1565 #ifdef FLOAT128
 1567 /*
 1568 -------------------------------------------------------------------------------
 1569 Returns the result of converting the single-precision floating-point value
 1570 `a' to the double-precision floating-point format.  The conversion is
 1571 performed according to the IEC/IEEE Standard for Binary Floating-Point
 1572 Arithmetic.
 1573 -------------------------------------------------------------------------------
 1574 */
float128 float32_to_float128( float32 a )
 1576 {
flag aSign;
int16 aExp;
bits32 aSig;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
 1584     if ( aExp == 0xFF ) {
 1585         if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a ) );
 1586         return packFloat128( aSign, 0x7FFF, 0, 0 );
 1587     }
 1588     if ( aExp == 0 ) {
 1589         if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 1591         --aExp;
 1592     }
 1593     return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );
 1595 }
 1597 #endif
 1599 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 1600 /*
 1601 -------------------------------------------------------------------------------
 1602 Rounds the single-precision floating-point value `a' to an integer, and
 1603 returns the result as a single-precision floating-point value.  The
 1604 operation is performed according to the IEC/IEEE Standard for Binary
 1605 Floating-Point Arithmetic.
 1606 -------------------------------------------------------------------------------
 1607 */
float32 float32_round_to_int( float32 a )
 1609 {
flag aSign;
int16 aExp;
bits32 lastBitMask, roundBitsMask;
int8 roundingMode;
float32 z;
aExp = extractFloat32Exp( a );
 1617     if ( 0x96 <= aExp ) {
 1618         if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
 1619             return propagateFloat32NaN( a, a );
 1620         }
 1621         return a;
 1622     }
 1623     if ( aExp <= 0x7E ) {
 1624         if ( (bits32) ( a<<1 ) == 0 ) return a;
float_set_inexact();
aSign = extractFloat32Sign( a );
 1627         switch ( float_rounding_mode() ) {
 1628          case float_round_nearest_even:
 1629             if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
 1630                 return packFloat32( aSign, 0x7F, 0 );
 1631             }
 1632             break;
 1633          case float_round_down:
 1634             return aSign ? 0xBF800000 : 0;
 1635          case float_round_up:
 1636             return aSign ? 0x80000000 : 0x3F800000;
 1637         }
 1638         return packFloat32( aSign, 0, 0 );
 1639     }
lastBitMask = 1;
lastBitMask <<= 0x96 - aExp;
roundBitsMask = lastBitMask - 1;
z = a;
roundingMode = float_rounding_mode();
 1645     if ( roundingMode == float_round_nearest_even ) {
z += lastBitMask>>1;
 1647         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
 1648     }
 1649     else if ( roundingMode != float_round_to_zero ) {
 1650         if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
z += roundBitsMask;
 1652         }
 1653     }
z &= ~ roundBitsMask;
 1655     if ( z != a ) float_set_inexact();
 1656     return z;
 1658 }
 1659 #endif /* !SOFTFLOAT_FOR_GCC */
 1661 /*
 1662 -------------------------------------------------------------------------------
 1663 Returns the result of adding the absolute values of the single-precision
 1664 floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
 1665 before being returned.  `zSign' is ignored if the result is a NaN.
 1666 The addition is performed according to the IEC/IEEE Standard for Binary
 1667 Floating-Point Arithmetic.
 1668 -------------------------------------------------------------------------------
 1669 */
 1670 static float32 addFloat32Sigs( float32 a, float32 b, flag zSign )
 1671 {
int16 aExp, bExp, zExp;
bits32 aSig, bSig, zSig;
int16 expDiff;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
bSig = extractFloat32Frac( b );
bExp = extractFloat32Exp( b );
expDiff = aExp - bExp;
aSig <<= 6;
bSig <<= 6;
 1683     if ( 0 < expDiff ) {
 1684         if ( aExp == 0xFF ) {
 1685             if ( aSig ) return propagateFloat32NaN( a, b );
 1686             return a;
 1687         }
 1688         if ( bExp == 0 ) {
 1689             --expDiff;
 1690         }
 1691         else {
bSig |= 0x20000000;
 1693         }
shift32RightJamming( bSig, expDiff, &bSig );
zExp = aExp;
 1696     }
 1697     else if ( expDiff < 0 ) {
 1698         if ( bExp == 0xFF ) {
 1699             if ( bSig ) return propagateFloat32NaN( a, b );
 1700             return packFloat32( zSign, 0xFF, 0 );
 1701         }
 1702         if ( aExp == 0 ) {
 1703             ++expDiff;
 1704         }
 1705         else {
aSig |= 0x20000000;
 1707         }
shift32RightJamming( aSig, - expDiff, &aSig );
zExp = bExp;
 1710     }
 1711     else {
 1712         if ( aExp == 0xFF ) {
 1713             if ( aSig | bSig ) return propagateFloat32NaN( a, b );
 1714             return a;
 1715         }
 1716         if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
zSig = 0x40000000 + aSig + bSig;
zExp = aExp;
 1719         goto roundAndPack;
 1720     }
aSig |= 0x20000000;
zSig = ( aSig + bSig )<<1;
 1723     --zExp;
 1724     if ( (sbits32) zSig < 0 ) {
zSig = aSig + bSig;
 1726         ++zExp;
 1727     }
roundAndPack:
 1729     return roundAndPackFloat32( zSign, zExp, zSig );
 1731 }
 1733 /*
 1734 -------------------------------------------------------------------------------
 1735 Returns the result of subtracting the absolute values of the single-
 1736 precision floating-point values `a' and `b'.  If `zSign' is 1, the
 1737 difference is negated before being returned.  `zSign' is ignored if the
 1738 result is a NaN.  The subtraction is performed according to the IEC/IEEE
 1739 Standard for Binary Floating-Point Arithmetic.
 1740 -------------------------------------------------------------------------------
 1741 */
 1742 static float32 subFloat32Sigs( float32 a, float32 b, flag zSign )
 1743 {
int16 aExp, bExp, zExp;
bits32 aSig, bSig, zSig;
int16 expDiff;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
bSig = extractFloat32Frac( b );
bExp = extractFloat32Exp( b );
expDiff = aExp - bExp;
aSig <<= 7;
bSig <<= 7;
 1755     if ( 0 < expDiff ) goto aExpBigger;
 1756     if ( expDiff < 0 ) goto bExpBigger;
 1757     if ( aExp == 0xFF ) {
 1758         if ( aSig | bSig ) return propagateFloat32NaN( a, b );
float_raise( float_flag_invalid );
 1760         return float32_default_nan;
 1761     }
 1762     if ( aExp == 0 ) {
aExp = 1;
bExp = 1;
 1765     }
 1766     if ( bSig < aSig ) goto aBigger;
 1767     if ( aSig < bSig ) goto bBigger;
 1768     return packFloat32( float_rounding_mode() == float_round_down, 0, 0 );
bExpBigger:
 1770     if ( bExp == 0xFF ) {
 1771         if ( bSig ) return propagateFloat32NaN( a, b );
 1772         return packFloat32( zSign ^ 1, 0xFF, 0 );
 1773     }
 1774     if ( aExp == 0 ) {
 1775         ++expDiff;
 1776     }
 1777     else {
aSig |= 0x40000000;
 1779     }
shift32RightJamming( aSig, - expDiff, &aSig );
bSig |= 0x40000000;
bBigger:
zSig = bSig - aSig;
zExp = bExp;
zSign ^= 1;
 1786     goto normalizeRoundAndPack;
aExpBigger:
 1788     if ( aExp == 0xFF ) {
 1789         if ( aSig ) return propagateFloat32NaN( a, b );
 1790         return a;
 1791     }
 1792     if ( bExp == 0 ) {
 1793         --expDiff;
 1794     }
 1795     else {
bSig |= 0x40000000;
 1797     }
shift32RightJamming( bSig, expDiff, &bSig );
aSig |= 0x40000000;
aBigger:
zSig = aSig - bSig;
zExp = aExp;
normalizeRoundAndPack:
 1804     --zExp;
 1805     return normalizeRoundAndPackFloat32( zSign, zExp, zSig );
 1807 }
 1809 /*
 1810 -------------------------------------------------------------------------------
 1811 Returns the result of adding the single-precision floating-point values `a'
 1812 and `b'.  The operation is performed according to the IEC/IEEE Standard for
 1813 Binary Floating-Point Arithmetic.
 1814 -------------------------------------------------------------------------------
 1815 */
float32 float32_add( float32 a, float32 b )
 1817 {
flag aSign, bSign;
aSign = extractFloat32Sign( a );
bSign = extractFloat32Sign( b );
 1822     if ( aSign == bSign ) {
 1823         return addFloat32Sigs( a, b, aSign );
 1824     }
 1825     else {
 1826         return subFloat32Sigs( a, b, aSign );
 1827     }
 1829 }
 1831 /*
 1832 -------------------------------------------------------------------------------
 1833 Returns the result of subtracting the single-precision floating-point values
 1834 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 1835 for Binary Floating-Point Arithmetic.
 1836 -------------------------------------------------------------------------------
 1837 */
float32 float32_sub( float32 a, float32 b )
 1839 {
flag aSign, bSign;
aSign = extractFloat32Sign( a );
bSign = extractFloat32Sign( b );
 1844     if ( aSign == bSign ) {
 1845         return subFloat32Sigs( a, b, aSign );
 1846     }
 1847     else {
 1848         return addFloat32Sigs( a, b, aSign );
 1849     }
 1851 }
 1853 /*
 1854 -------------------------------------------------------------------------------
 1855 Returns the result of multiplying the single-precision floating-point values
 1856 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 1857 for Binary Floating-Point Arithmetic.
 1858 -------------------------------------------------------------------------------
 1859 */
float32 float32_mul( float32 a, float32 b )
 1861 {
flag aSign, bSign, zSign;
int16 aExp, bExp, zExp;
bits32 aSig, bSig;
bits64 zSig64;
bits32 zSig;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
bSig = extractFloat32Frac( b );
bExp = extractFloat32Exp( b );
bSign = extractFloat32Sign( b );
zSign = aSign ^ bSign;
 1875     if ( aExp == 0xFF ) {
 1876         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
 1877             return propagateFloat32NaN( a, b );
 1878         }
 1879         if ( ( bExp | bSig ) == 0 ) {
float_raise( float_flag_invalid );
 1881             return float32_default_nan;
 1882         }
 1883         return packFloat32( zSign, 0xFF, 0 );
 1884     }
 1885     if ( bExp == 0xFF ) {
 1886         if ( bSig ) return propagateFloat32NaN( a, b );
 1887         if ( ( aExp | aSig ) == 0 ) {
float_raise( float_flag_invalid );
 1889             return float32_default_nan;
 1890         }
 1891         return packFloat32( zSign, 0xFF, 0 );
 1892     }
 1893     if ( aExp == 0 ) {
 1894         if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 1896     }
 1897     if ( bExp == 0 ) {
 1898         if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
normalizeFloat32Subnormal( bSig, &bExp, &bSig );
 1900     }
zExp = aExp + bExp - 0x7F;
aSig = ( aSig | 0x00800000 )<<7;
bSig = ( bSig | 0x00800000 )<<8;
shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
zSig = zSig64;
 1906     if ( 0 <= (sbits32) ( zSig<<1 ) ) {
zSig <<= 1;
 1908         --zExp;
 1909     }
 1910     return roundAndPackFloat32( zSign, zExp, zSig );
 1912 }
 1914 /*
 1915 -------------------------------------------------------------------------------
 1916 Returns the result of dividing the single-precision floating-point value `a'
 1917 by the corresponding value `b'.  The operation is performed according to the
 1918 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1919 -------------------------------------------------------------------------------
 1920 */
float32 float32_div( float32 a, float32 b )
 1922 {
flag aSign, bSign, zSign;
int16 aExp, bExp, zExp;
bits32 aSig, bSig, zSig;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
bSig = extractFloat32Frac( b );
bExp = extractFloat32Exp( b );
bSign = extractFloat32Sign( b );
zSign = aSign ^ bSign;
 1934     if ( aExp == 0xFF ) {
 1935         if ( aSig ) return propagateFloat32NaN( a, b );
 1936         if ( bExp == 0xFF ) {
 1937             if ( bSig ) return propagateFloat32NaN( a, b );
float_raise( float_flag_invalid );
 1939             return float32_default_nan;
 1940         }
 1941         return packFloat32( zSign, 0xFF, 0 );
 1942     }
 1943     if ( bExp == 0xFF ) {
 1944         if ( bSig ) return propagateFloat32NaN( a, b );
 1945         return packFloat32( zSign, 0, 0 );
 1946     }
 1947     if ( bExp == 0 ) {
 1948         if ( bSig == 0 ) {
 1949             if ( ( aExp | aSig ) == 0 ) {
float_raise( float_flag_invalid );
 1951                 return float32_default_nan;
 1952             }
float_raise( float_flag_divbyzero );
 1954             return packFloat32( zSign, 0xFF, 0 );
 1955         }
normalizeFloat32Subnormal( bSig, &bExp, &bSig );
 1957     }
 1958     if ( aExp == 0 ) {
 1959         if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 1961     }
zExp = aExp - bExp + 0x7D;
aSig = ( aSig | 0x00800000 )<<7;
bSig = ( bSig | 0x00800000 )<<8;
 1965     if ( bSig <= ( aSig + aSig ) ) {
aSig >>= 1;
 1967         ++zExp;
 1968     }
zSig = ( ( (bits64) aSig )<<32 ) / bSig;
 1970     if ( ( zSig & 0x3F ) == 0 ) {
zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );
 1972     }
 1973     return roundAndPackFloat32( zSign, zExp, zSig );
 1975 }
 1977 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 1978 /*
 1979 -------------------------------------------------------------------------------
 1980 Returns the remainder of the single-precision floating-point value `a'
 1981 with respect to the corresponding value `b'.  The operation is performed
 1982 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 1983 -------------------------------------------------------------------------------
 1984 */
float32 float32_rem( float32 a, float32 b )
 1986 {
flag aSign, bSign, zSign;
int16 aExp, bExp, expDiff;
bits32 aSig, bSig;
bits32 q;
bits64 aSig64, bSig64, q64;
bits32 alternateASig;
sbits32 sigMean;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
bSig = extractFloat32Frac( b );
bExp = extractFloat32Exp( b );
bSign = extractFloat32Sign( b );
 2001     if ( aExp == 0xFF ) {
 2002         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
 2003             return propagateFloat32NaN( a, b );
 2004         }
float_raise( float_flag_invalid );
 2006         return float32_default_nan;
 2007     }
 2008     if ( bExp == 0xFF ) {
 2009         if ( bSig ) return propagateFloat32NaN( a, b );
 2010         return a;
 2011     }
 2012     if ( bExp == 0 ) {
 2013         if ( bSig == 0 ) {
float_raise( float_flag_invalid );
 2015             return float32_default_nan;
 2016         }
normalizeFloat32Subnormal( bSig, &bExp, &bSig );
 2018     }
 2019     if ( aExp == 0 ) {
 2020         if ( aSig == 0 ) return a;
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 2022     }
expDiff = aExp - bExp;
aSig |= 0x00800000;
bSig |= 0x00800000;
 2026     if ( expDiff < 32 ) {
aSig <<= 8;
bSig <<= 8;
 2029         if ( expDiff < 0 ) {
 2030             if ( expDiff < -1 ) return a;
aSig >>= 1;
 2032         }
q = ( bSig <= aSig );
 2034         if ( q ) aSig -= bSig;
 2035         if ( 0 < expDiff ) {
q = ( ( (bits64) aSig )<<32 ) / bSig;
q >>= 32 - expDiff;
bSig >>= 2;
aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
 2040         }
 2041         else {
aSig >>= 2;
bSig >>= 2;
 2044         }
 2045     }
 2046     else {
 2047         if ( bSig <= aSig ) aSig -= bSig;
aSig64 = ( (bits64) aSig )<<40;
bSig64 = ( (bits64) bSig )<<40;
expDiff -= 64;
 2051         while ( 0 < expDiff ) {
q64 = estimateDiv128To64( aSig64, 0, bSig64 );
q64 = ( 2 < q64 ) ? q64 - 2 : 0;
aSig64 = - ( ( bSig * q64 )<<38 );
expDiff -= 62;
 2056         }
expDiff += 64;
q64 = estimateDiv128To64( aSig64, 0, bSig64 );
q64 = ( 2 < q64 ) ? q64 - 2 : 0;
q = q64>>( 64 - expDiff );
bSig <<= 6;
aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;
 2063     }
 2064     do {
alternateASig = aSig;
 2066         ++q;
aSig -= bSig;
 2068     } while ( 0 <= (sbits32) aSig );
sigMean = aSig + alternateASig;
 2070     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
aSig = alternateASig;
 2072     }
zSign = ( (sbits32) aSig < 0 );
 2074     if ( zSign ) aSig = - aSig;
 2075     return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig );
 2077 }
 2078 #endif /* !SOFTFLOAT_FOR_GCC */
 2080 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 2081 /*
 2082 -------------------------------------------------------------------------------
 2083 Returns the square root of the single-precision floating-point value `a'.
 2084 The operation is performed according to the IEC/IEEE Standard for Binary
 2085 Floating-Point Arithmetic.
 2086 -------------------------------------------------------------------------------
 2087 */
float32 float32_sqrt( float32 a )
 2089 {
flag aSign;
int16 aExp, zExp;
bits32 aSig, zSig;
bits64 rem, term;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
 2098     if ( aExp == 0xFF ) {
 2099         if ( aSig ) return propagateFloat32NaN( a, 0 );
 2100         if ( ! aSign ) return a;
float_raise( float_flag_invalid );
 2102         return float32_default_nan;
 2103     }
 2104     if ( aSign ) {
 2105         if ( ( aExp | aSig ) == 0 ) return a;
float_raise( float_flag_invalid );
 2107         return float32_default_nan;
 2108     }
 2109     if ( aExp == 0 ) {
 2110         if ( aSig == 0 ) return 0;
normalizeFloat32Subnormal( aSig, &aExp, &aSig );
 2112     }
zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
aSig = ( aSig | 0x00800000 )<<8;
zSig = estimateSqrt32( aExp, aSig ) + 2;
 2116     if ( ( zSig & 0x7F ) <= 5 ) {
 2117         if ( zSig < 2 ) {
zSig = 0x7FFFFFFF;
 2119             goto roundAndPack;
 2120         }
aSig >>= aExp & 1;
term = ( (bits64) zSig ) * zSig;
rem = ( ( (bits64) aSig )<<32 ) - term;
 2124         while ( (sbits64) rem < 0 ) {
 2125             --zSig;
rem += ( ( (bits64) zSig )<<1 ) | 1;
 2127         }
zSig |= ( rem != 0 );
 2129     }
shift32RightJamming( zSig, 1, &zSig );
roundAndPack:
 2132     return roundAndPackFloat32( 0, zExp, zSig );
 2134 }
 2135 #endif /* !SOFTFLOAT_FOR_GCC */
 2137 /*
 2138 -------------------------------------------------------------------------------
 2139 Returns 1 if the single-precision floating-point value `a' is equal to
 2140 the corresponding value `b', and 0 otherwise.  The comparison is performed
 2141 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 2142 -------------------------------------------------------------------------------
 2143 */
flag float32_eq( float32 a, float32 b )
 2145 {
 2147     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
 2148          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
 2149        ) {
 2150         if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 2152         }
 2153         return 0;
 2154     }
 2155     return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
 2157 }
 2159 /*
 2160 -------------------------------------------------------------------------------
 2161 Returns 1 if the single-precision floating-point value `a' is less than
 2162 or equal to the corresponding value `b', and 0 otherwise.  The comparison
 2163 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 2164 Arithmetic.
 2165 -------------------------------------------------------------------------------
 2166 */
flag float32_le( float32 a, float32 b )
 2168 {
flag aSign, bSign;
 2171     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
 2172          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
 2173        ) {
float_raise( float_flag_invalid );
 2175         return 0;
 2176     }
aSign = extractFloat32Sign( a );
bSign = extractFloat32Sign( b );
 2179     if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
 2180     return ( a == b ) || ( aSign ^ ( a < b ) );
 2182 }
 2184 /*
 2185 -------------------------------------------------------------------------------
 2186 Returns 1 if the single-precision floating-point value `a' is less than
 2187 the corresponding value `b', and 0 otherwise.  The comparison is performed
 2188 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 2189 -------------------------------------------------------------------------------
 2190 */
flag float32_lt( float32 a, float32 b )
 2192 {
flag aSign, bSign;
 2195     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
 2196          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
 2197        ) {
float_raise( float_flag_invalid );
 2199         return 0;
 2200     }
aSign = extractFloat32Sign( a );
bSign = extractFloat32Sign( b );
 2203     if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
 2204     return ( a != b ) && ( aSign ^ ( a < b ) );
 2206 }
 2208 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 2209 /*
 2210 -------------------------------------------------------------------------------
 2211 Returns 1 if the single-precision floating-point value `a' is equal to
 2212 the corresponding value `b', and 0 otherwise.  The invalid exception is
 2213 raised if either operand is a NaN.  Otherwise, the comparison is performed
 2214 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 2215 -------------------------------------------------------------------------------
 2216 */
flag float32_eq_signaling( float32 a, float32 b )
 2218 {
 2220     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
 2221          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
 2222        ) {
float_raise( float_flag_invalid );
 2224         return 0;
 2225     }
 2226     return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
 2228 }
 2230 /*
 2231 -------------------------------------------------------------------------------
 2232 Returns 1 if the single-precision floating-point value `a' is less than or
 2233 equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
 2234 cause an exception.  Otherwise, the comparison is performed according to the
 2235 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 2236 -------------------------------------------------------------------------------
 2237 */
flag float32_le_quiet( float32 a, float32 b )
 2239 {
flag aSign, bSign;
 2242     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
 2243          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
 2244        ) {
 2245         if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 2247         }
 2248         return 0;
 2249     }
aSign = extractFloat32Sign( a );
bSign = extractFloat32Sign( b );
 2252     if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
 2253     return ( a == b ) || ( aSign ^ ( a < b ) );
 2255 }
 2257 /*
 2258 -------------------------------------------------------------------------------
 2259 Returns 1 if the single-precision floating-point value `a' is less than
 2260 the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
 2261 exception.  Otherwise, the comparison is performed according to the IEC/IEEE
 2262 Standard for Binary Floating-Point Arithmetic.
 2263 -------------------------------------------------------------------------------
 2264 */
flag float32_lt_quiet( float32 a, float32 b )
 2266 {
flag aSign, bSign;
 2269     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
 2270          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
 2271        ) {
 2272         if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 2274         }
 2275         return 0;
 2276     }
aSign = extractFloat32Sign( a );
bSign = extractFloat32Sign( b );
 2279     if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
 2280     return ( a != b ) && ( aSign ^ ( a < b ) );
 2282 }
 2283 #endif /* !SOFTFLOAT_FOR_GCC */
 2285 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 2286 /*
 2287 -------------------------------------------------------------------------------
 2288 Returns the result of converting the double-precision floating-point value
 2289 `a' to the 32-bit two's complement integer format.  The conversion is
 2290 performed according to the IEC/IEEE Standard for Binary Floating-Point
 2291 Arithmetic---which means in particular that the conversion is rounded
 2292 according to the current rounding mode.  If `a' is a NaN, the largest
 2293 positive integer is returned.  Otherwise, if the conversion overflows, the
 2294 largest integer with the same sign as `a' is returned.
 2295 -------------------------------------------------------------------------------
 2296 */
int32 float64_to_int32( float64 a )
 2298 {
flag aSign;
int16 aExp, shiftCount;
bits64 aSig;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2306     if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
 2307     if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
shiftCount = 0x42C - aExp;
 2309     if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
 2310     return roundAndPackInt32( aSign, aSig );
 2312 }
 2313 #endif /* !SOFTFLOAT_FOR_GCC */
 2315 /*
 2316 -------------------------------------------------------------------------------
 2317 Returns the result of converting the double-precision floating-point value
 2318 `a' to the 32-bit two's complement integer format.  The conversion is
 2319 performed according to the IEC/IEEE Standard for Binary Floating-Point
 2320 Arithmetic, except that the conversion is always rounded toward zero.
 2321 If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
 2322 the conversion overflows, the largest integer with the same sign as `a' is
 2323 returned.
 2324 -------------------------------------------------------------------------------
 2325 */
int32 float64_to_int32_round_to_zero( float64 a )
 2327 {
flag aSign;
int16 aExp, shiftCount;
bits64 aSig, savedASig;
int32 z;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2336     if ( 0x41E < aExp ) {
 2337         if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
 2338         goto invalid;
 2339     }
 2340     else if ( aExp < 0x3FF ) {
 2341         if ( aExp || aSig ) float_set_inexact();
 2342         return 0;
 2343     }
aSig |= LIT64( 0x0010000000000000 );
shiftCount = 0x433 - aExp;
savedASig = aSig;
aSig >>= shiftCount;
z = aSig;
 2349     if ( aSign ) z = - z;
 2350     if ( ( z < 0 ) ^ aSign ) {
invalid:
float_raise( float_flag_invalid );
 2353         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
 2354     }
 2355     if ( ( aSig<<shiftCount ) != savedASig ) {
float_set_inexact();
 2357     }
 2358     return z;
 2360 }
 2362 #ifndef SOFTFLOAT_FOR_GCC /* Not needed */
 2363 /*
 2364 -------------------------------------------------------------------------------
 2365 Returns the result of converting the double-precision floating-point value
 2366 `a' to the 64-bit two's complement integer format.  The conversion is
 2367 performed according to the IEC/IEEE Standard for Binary Floating-Point
 2368 Arithmetic---which means in particular that the conversion is rounded
 2369 according to the current rounding mode.  If `a' is a NaN, the largest
 2370 positive integer is returned.  Otherwise, if the conversion overflows, the
 2371 largest integer with the same sign as `a' is returned.
 2372 -------------------------------------------------------------------------------
 2373 */
int64 float64_to_int64( float64 a )
 2375 {
flag aSign;
int16 aExp, shiftCount;
bits64 aSig, aSigExtra;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2383     if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
shiftCount = 0x433 - aExp;
 2385     if ( shiftCount <= 0 ) {
 2386         if ( 0x43E < aExp ) {
float_raise( float_flag_invalid );
 2388             if (    ! aSign
 2389                  || (    ( aExp == 0x7FF )
 2390                       && ( aSig != LIT64( 0x0010000000000000 ) ) )
 2391                ) {
 2392                 return LIT64( 0x7FFFFFFFFFFFFFFF );
 2393             }
 2394             return (sbits64) LIT64( 0x8000000000000000 );
 2395         }
aSigExtra = 0;
aSig <<= - shiftCount;
 2398     }
 2399     else {
shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
 2401     }
 2402     return roundAndPackInt64( aSign, aSig, aSigExtra );
 2404 }
 2406 /*
 2407 -------------------------------------------------------------------------------
 2408 Returns the result of converting the double-precision floating-point value
 2409 `a' to the 64-bit two's complement integer format.  The conversion is
 2410 performed according to the IEC/IEEE Standard for Binary Floating-Point
 2411 Arithmetic, except that the conversion is always rounded toward zero.
 2412 If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
 2413 the conversion overflows, the largest integer with the same sign as `a' is
 2414 returned.
 2415 -------------------------------------------------------------------------------
 2416 */
int64 float64_to_int64_round_to_zero( float64 a )
 2418 {
flag aSign;
int16 aExp, shiftCount;
bits64 aSig;
int64 z;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2427     if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
shiftCount = aExp - 0x433;
 2429     if ( 0 <= shiftCount ) {
 2430         if ( 0x43E <= aExp ) {
 2431             if ( a != LIT64( 0xC3E0000000000000 ) ) {
float_raise( float_flag_invalid );
 2433                 if (    ! aSign
 2434                      || (    ( aExp == 0x7FF )
 2435                           && ( aSig != LIT64( 0x0010000000000000 ) ) )
 2436                    ) {
 2437                     return LIT64( 0x7FFFFFFFFFFFFFFF );
 2438                 }
 2439             }
 2440             return (sbits64) LIT64( 0x8000000000000000 );
 2441         }
z = aSig<<shiftCount;
 2443     }
 2444     else {
 2445         if ( aExp < 0x3FE ) {
 2446             if ( aExp | aSig ) float_set_inexact();
 2447             return 0;
 2448         }
z = aSig>>( - shiftCount );
 2450         if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
float_set_inexact();
 2452         }
 2453     }
 2454     if ( aSign ) z = - z;
 2455     return z;
 2457 }
 2458 #endif /* !SOFTFLOAT_FOR_GCC */
 2460 /*
 2461 -------------------------------------------------------------------------------
 2462 Returns the result of converting the double-precision floating-point value
 2463 `a' to the single-precision floating-point format.  The conversion is
 2464 performed according to the IEC/IEEE Standard for Binary Floating-Point
 2465 Arithmetic.
 2466 -------------------------------------------------------------------------------
 2467 */
float32 float64_to_float32( float64 a )
 2469 {
flag aSign;
int16 aExp;
bits64 aSig;
bits32 zSig;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2478     if ( aExp == 0x7FF ) {
 2479         if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a ) );
 2480         return packFloat32( aSign, 0xFF, 0 );
 2481     }
shift64RightJamming( aSig, 22, &aSig );
zSig = aSig;
 2484     if ( aExp || zSig ) {
zSig |= 0x40000000;
aExp -= 0x381;
 2487     }
 2488     return roundAndPackFloat32( aSign, aExp, zSig );
 2490 }
 2492 #ifdef FLOATX80
 2494 /*
 2495 -------------------------------------------------------------------------------
 2496 Returns the result of converting the double-precision floating-point value
 2497 `a' to the extended double-precision floating-point format.  The conversion
 2498 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 2499 Arithmetic.
 2500 -------------------------------------------------------------------------------
 2501 */
floatx80 float64_to_floatx80( float64 a )
 2503 {
flag aSign;
int16 aExp;
bits64 aSig;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2511     if ( aExp == 0x7FF ) {
 2512         if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a ) );
 2513         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 2514     }
 2515     if ( aExp == 0 ) {
 2516         if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
normalizeFloat64Subnormal( aSig, &aExp, &aSig );
 2518     }
 2519     return
packFloatx80(
aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
 2523 }
 2525 #endif
 2527 #ifdef FLOAT128
 2529 /*
 2530 -------------------------------------------------------------------------------
 2531 Returns the result of converting the double-precision floating-point value
 2532 `a' to the quadruple-precision floating-point format.  The conversion is
 2533 performed according to the IEC/IEEE Standard for Binary Floating-Point
 2534 Arithmetic.
 2535 -------------------------------------------------------------------------------
 2536 */
float128 float64_to_float128( float64 a )
 2538 {
flag aSign;
int16 aExp;
bits64 aSig, zSig0, zSig1;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 2546     if ( aExp == 0x7FF ) {
 2547         if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a ) );
 2548         return packFloat128( aSign, 0x7FFF, 0, 0 );
 2549     }
 2550     if ( aExp == 0 ) {
 2551         if ( aSig == 0 ) return packFloat128( aSign, 0, 0, 0 );
normalizeFloat64Subnormal( aSig, &aExp, &aSig );
 2553         --aExp;
 2554     }
shift128Right( aSig, 0, 4, &zSig0, &zSig1 );
 2556     return packFloat128( aSign, aExp + 0x3C00, zSig0, zSig1 );
 2558 }
 2560 #endif
 2562 #ifndef SOFTFLOAT_FOR_GCC
 2563 /*
 2564 -------------------------------------------------------------------------------
 2565 Rounds the double-precision floating-point value `a' to an integer, and
 2566 returns the result as a double-precision floating-point value.  The
 2567 operation is performed according to the IEC/IEEE Standard for Binary
 2568 Floating-Point Arithmetic.
 2569 -------------------------------------------------------------------------------
 2570 */
float64 float64_round_to_int( float64 a )
 2572 {
flag aSign;
int16 aExp;
bits64 lastBitMask, roundBitsMask;
int8 roundingMode;
float64 z;
aExp = extractFloat64Exp( a );
 2580     if ( 0x433 <= aExp ) {
 2581         if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {
 2582             return propagateFloat64NaN( a, a );
 2583         }
 2584         return a;
 2585     }
 2586     if ( aExp < 0x3FF ) {
 2587         if ( (bits64) ( a<<1 ) == 0 ) return a;
float_set_inexact();
aSign = extractFloat64Sign( a );
 2590         switch ( float_rounding_mode() ) {
 2591          case float_round_nearest_even:
 2592             if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {
 2593                 return packFloat64( aSign, 0x3FF, 0 );
 2594             }
 2595             break;
 2596          case float_round_down:
 2597             return aSign ? LIT64( 0xBFF0000000000000 ) : 0;
 2598          case float_round_up:
 2599             return
aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );
 2601         }
 2602         return packFloat64( aSign, 0, 0 );
 2603     }
lastBitMask = 1;
lastBitMask <<= 0x433 - aExp;
roundBitsMask = lastBitMask - 1;
z = a;
roundingMode = float_rounding_mode();
 2609     if ( roundingMode == float_round_nearest_even ) {
z += lastBitMask>>1;
 2611         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
 2612     }
 2613     else if ( roundingMode != float_round_to_zero ) {
 2614         if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {
z += roundBitsMask;
 2616         }
 2617     }
z &= ~ roundBitsMask;
 2619     if ( z != a ) float_set_inexact();
 2620     return z;
 2622 }
 2623 #endif
 2625 /*
 2626 -------------------------------------------------------------------------------
 2627 Returns the result of adding the absolute values of the double-precision
 2628 floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
 2629 before being returned.  `zSign' is ignored if the result is a NaN.
 2630 The addition is performed according to the IEC/IEEE Standard for Binary
 2631 Floating-Point Arithmetic.
 2632 -------------------------------------------------------------------------------
 2633 */
 2634 static float64 addFloat64Sigs( float64 a, float64 b, flag zSign )
 2635 {
int16 aExp, bExp, zExp;
bits64 aSig, bSig, zSig;
int16 expDiff;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
bSig = extractFloat64Frac( b );
bExp = extractFloat64Exp( b );
expDiff = aExp - bExp;
aSig <<= 9;
bSig <<= 9;
 2647     if ( 0 < expDiff ) {
 2648         if ( aExp == 0x7FF ) {
 2649             if ( aSig ) return propagateFloat64NaN( a, b );
 2650             return a;
 2651         }
 2652         if ( bExp == 0 ) {
 2653             --expDiff;
 2654         }
 2655         else {
bSig |= LIT64( 0x2000000000000000 );
 2657         }
shift64RightJamming( bSig, expDiff, &bSig );
zExp = aExp;
 2660     }
 2661     else if ( expDiff < 0 ) {
 2662         if ( bExp == 0x7FF ) {
 2663             if ( bSig ) return propagateFloat64NaN( a, b );
 2664             return packFloat64( zSign, 0x7FF, 0 );
 2665         }
 2666         if ( aExp == 0 ) {
 2667             ++expDiff;
 2668         }
 2669         else {
aSig |= LIT64( 0x2000000000000000 );
 2671         }
shift64RightJamming( aSig, - expDiff, &aSig );
zExp = bExp;
 2674     }
 2675     else {
 2676         if ( aExp == 0x7FF ) {
 2677             if ( aSig | bSig ) return propagateFloat64NaN( a, b );
 2678             return a;
 2679         }
 2680         if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
zExp = aExp;
 2683         goto roundAndPack;
 2684     }
aSig |= LIT64( 0x2000000000000000 );
zSig = ( aSig + bSig )<<1;
 2687     --zExp;
 2688     if ( (sbits64) zSig < 0 ) {
zSig = aSig + bSig;
 2690         ++zExp;
 2691     }
roundAndPack:
 2693     return roundAndPackFloat64( zSign, zExp, zSig );
 2695 }
 2697 /*
 2698 -------------------------------------------------------------------------------
 2699 Returns the result of subtracting the absolute values of the double-
 2700 precision floating-point values `a' and `b'.  If `zSign' is 1, the
 2701 difference is negated before being returned.  `zSign' is ignored if the
 2702 result is a NaN.  The subtraction is performed according to the IEC/IEEE
 2703 Standard for Binary Floating-Point Arithmetic.
 2704 -------------------------------------------------------------------------------
 2705 */
 2706 static float64 subFloat64Sigs( float64 a, float64 b, flag zSign )
 2707 {
int16 aExp, bExp, zExp;
bits64 aSig, bSig, zSig;
int16 expDiff;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
bSig = extractFloat64Frac( b );
bExp = extractFloat64Exp( b );
expDiff = aExp - bExp;
aSig <<= 10;
bSig <<= 10;
 2719     if ( 0 < expDiff ) goto aExpBigger;
 2720     if ( expDiff < 0 ) goto bExpBigger;
 2721     if ( aExp == 0x7FF ) {
 2722         if ( aSig | bSig ) return propagateFloat64NaN( a, b );
float_raise( float_flag_invalid );
 2724         return float64_default_nan;
 2725     }
 2726     if ( aExp == 0 ) {
aExp = 1;
bExp = 1;
 2729     }
 2730     if ( bSig < aSig ) goto aBigger;
 2731     if ( aSig < bSig ) goto bBigger;
 2732     return packFloat64( float_rounding_mode() == float_round_down, 0, 0 );
bExpBigger:
 2734     if ( bExp == 0x7FF ) {
 2735         if ( bSig ) return propagateFloat64NaN( a, b );
 2736         return packFloat64( zSign ^ 1, 0x7FF, 0 );
 2737     }
 2738     if ( aExp == 0 ) {
 2739         ++expDiff;
 2740     }
 2741     else {
aSig |= LIT64( 0x4000000000000000 );
 2743     }
shift64RightJamming( aSig, - expDiff, &aSig );
bSig |= LIT64( 0x4000000000000000 );
bBigger:
zSig = bSig - aSig;
zExp = bExp;
zSign ^= 1;
 2750     goto normalizeRoundAndPack;
aExpBigger:
 2752     if ( aExp == 0x7FF ) {
 2753         if ( aSig ) return propagateFloat64NaN( a, b );
 2754         return a;
 2755     }
 2756     if ( bExp == 0 ) {
 2757         --expDiff;
 2758     }
 2759     else {
bSig |= LIT64( 0x4000000000000000 );
 2761     }
shift64RightJamming( bSig, expDiff, &bSig );
aSig |= LIT64( 0x4000000000000000 );
aBigger:
zSig = aSig - bSig;
zExp = aExp;
normalizeRoundAndPack:
 2768     --zExp;
 2769     return normalizeRoundAndPackFloat64( zSign, zExp, zSig );
 2771 }
 2773 /*
 2774 -------------------------------------------------------------------------------
 2775 Returns the result of adding the double-precision floating-point values `a'
 2776 and `b'.  The operation is performed according to the IEC/IEEE Standard for
 2777 Binary Floating-Point Arithmetic.
 2778 -------------------------------------------------------------------------------
 2779 */
float64 float64_add( float64 a, float64 b )
 2781 {
flag aSign, bSign;
aSign = extractFloat64Sign( a );
bSign = extractFloat64Sign( b );
 2786     if ( aSign == bSign ) {
 2787         return addFloat64Sigs( a, b, aSign );
 2788     }
 2789     else {
 2790         return subFloat64Sigs( a, b, aSign );
 2791     }
 2793 }
 2795 /*
 2796 -------------------------------------------------------------------------------
 2797 Returns the result of subtracting the double-precision floating-point values
 2798 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 2799 for Binary Floating-Point Arithmetic.
 2800 -------------------------------------------------------------------------------
 2801 */
float64 float64_sub( float64 a, float64 b )
 2803 {
flag aSign, bSign;
aSign = extractFloat64Sign( a );
bSign = extractFloat64Sign( b );
 2808     if ( aSign == bSign ) {
 2809         return subFloat64Sigs( a, b, aSign );
 2810     }
 2811     else {
 2812         return addFloat64Sigs( a, b, aSign );
 2813     }
 2815 }
 2817 /*
 2818 -------------------------------------------------------------------------------
 2819 Returns the result of multiplying the double-precision floating-point values
 2820 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 2821 for Binary Floating-Point Arithmetic.
 2822 -------------------------------------------------------------------------------
 2823 */
float64 float64_mul( float64 a, float64 b )
 2825 {
flag aSign, bSign, zSign;
int16 aExp, bExp, zExp;
bits64 aSig, bSig, zSig0, zSig1;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
bSig = extractFloat64Frac( b );
bExp = extractFloat64Exp( b );
bSign = extractFloat64Sign( b );
zSign = aSign ^ bSign;
 2837     if ( aExp == 0x7FF ) {
 2838         if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
 2839             return propagateFloat64NaN( a, b );
 2840         }
 2841         if ( ( bExp | bSig ) == 0 ) {
float_raise( float_flag_invalid );
 2843             return float64_default_nan;
 2844         }
 2845         return packFloat64( zSign, 0x7FF, 0 );
 2846     }
 2847     if ( bExp == 0x7FF ) {
 2848         if ( bSig ) return propagateFloat64NaN( a, b );
 2849         if ( ( aExp | aSig ) == 0 ) {
float_raise( float_flag_invalid );
 2851             return float64_default_nan;
 2852         }
 2853         return packFloat64( zSign, 0x7FF, 0 );
 2854     }
 2855     if ( aExp == 0 ) {
 2856         if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
normalizeFloat64Subnormal( aSig, &aExp, &aSig );
 2858     }
 2859     if ( bExp == 0 ) {
 2860         if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );
normalizeFloat64Subnormal( bSig, &bExp, &bSig );
 2862     }
zExp = aExp + bExp - 0x3FF;
aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
mul64To128( aSig, bSig, &zSig0, &zSig1 );
zSig0 |= ( zSig1 != 0 );
 2868     if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
zSig0 <<= 1;
 2870         --zExp;
 2871     }
 2872     return roundAndPackFloat64( zSign, zExp, zSig0 );
 2874 }
 2876 /*
 2877 -------------------------------------------------------------------------------
 2878 Returns the result of dividing the double-precision floating-point value `a'
 2879 by the corresponding value `b'.  The operation is performed according to
 2880 the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 2881 -------------------------------------------------------------------------------
 2882 */
float64 float64_div( float64 a, float64 b )
 2884 {
flag aSign, bSign, zSign;
int16 aExp, bExp, zExp;
bits64 aSig, bSig, zSig;
bits64 rem0, rem1;
bits64 term0, term1;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
bSig = extractFloat64Frac( b );
bExp = extractFloat64Exp( b );
bSign = extractFloat64Sign( b );
zSign = aSign ^ bSign;
 2898     if ( aExp == 0x7FF ) {
 2899         if ( aSig ) return propagateFloat64NaN( a, b );
 2900         if ( bExp == 0x7FF ) {
 2901             if ( bSig ) return propagateFloat64NaN( a, b );
float_raise( float_flag_invalid );
 2903             return float64_default_nan;
 2904         }
 2905         return packFloat64( zSign, 0x7FF, 0 );
 2906     }
 2907     if ( bExp == 0x7FF ) {
 2908         if ( bSig ) return propagateFloat64NaN( a, b );
 2909         return packFloat64( zSign, 0, 0 );
 2910     }
 2911     if ( bExp == 0 ) {
 2912         if ( bSig == 0 ) {
 2913             if ( ( aExp | aSig ) == 0 ) {
float_raise( float_flag_invalid );
 2915                 return float64_default_nan;
 2916             }
float_raise( float_flag_divbyzero );
 2918             return packFloat64( zSign, 0x7FF, 0 );
 2919         }
normalizeFloat64Subnormal( bSig, &bExp, &bSig );
 2921     }
 2922     if ( aExp == 0 ) {
 2923         if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
normalizeFloat64Subnormal( aSig, &aExp, &aSig );
 2925     }
zExp = aExp - bExp + 0x3FD;
aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
 2929     if ( bSig <= ( aSig + aSig ) ) {
aSig >>= 1;
 2931         ++zExp;
 2932     }
zSig = estimateDiv128To64( aSig, 0, bSig );
 2934     if ( ( zSig & 0x1FF ) <= 2 ) {
mul64To128( bSig, zSig, &term0, &term1 );
sub128( aSig, 0, term0, term1, &rem0, &rem1 );
 2937         while ( (sbits64) rem0 < 0 ) {
 2938             --zSig;
add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
 2940         }
zSig |= ( rem1 != 0 );
 2942     }
 2943     return roundAndPackFloat64( zSign, zExp, zSig );
 2945 }
 2947 #ifndef SOFTFLOAT_FOR_GCC
 2948 /*
 2949 -------------------------------------------------------------------------------
 2950 Returns the remainder of the double-precision floating-point value `a'
 2951 with respect to the corresponding value `b'.  The operation is performed
 2952 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 2953 -------------------------------------------------------------------------------
 2954 */
float64 float64_rem( float64 a, float64 b )
 2956 {
flag aSign, bSign, zSign;
int16 aExp, bExp, expDiff;
bits64 aSig, bSig;
bits64 q, alternateASig;
sbits64 sigMean;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
bSig = extractFloat64Frac( b );
bExp = extractFloat64Exp( b );
bSign = extractFloat64Sign( b );
 2969     if ( aExp == 0x7FF ) {
 2970         if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
 2971             return propagateFloat64NaN( a, b );
 2972         }
float_raise( float_flag_invalid );
 2974         return float64_default_nan;
 2975     }
 2976     if ( bExp == 0x7FF ) {
 2977         if ( bSig ) return propagateFloat64NaN( a, b );
 2978         return a;
 2979     }
 2980     if ( bExp == 0 ) {
 2981         if ( bSig == 0 ) {
float_raise( float_flag_invalid );
 2983             return float64_default_nan;
 2984         }
normalizeFloat64Subnormal( bSig, &bExp, &bSig );
 2986     }
 2987     if ( aExp == 0 ) {
 2988         if ( aSig == 0 ) return a;
normalizeFloat64Subnormal( aSig, &aExp, &aSig );
 2990     }
expDiff = aExp - bExp;
aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;
bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
 2994     if ( expDiff < 0 ) {
 2995         if ( expDiff < -1 ) return a;
aSig >>= 1;
 2997     }
q = ( bSig <= aSig );
 2999     if ( q ) aSig -= bSig;
expDiff -= 64;
 3001     while ( 0 < expDiff ) {
q = estimateDiv128To64( aSig, 0, bSig );
q = ( 2 < q ) ? q - 2 : 0;
aSig = - ( ( bSig>>2 ) * q );
expDiff -= 62;
 3006     }
expDiff += 64;
 3008     if ( 0 < expDiff ) {
q = estimateDiv128To64( aSig, 0, bSig );
q = ( 2 < q ) ? q - 2 : 0;
q >>= 64 - expDiff;
bSig >>= 2;
aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
 3014     }
 3015     else {
aSig >>= 2;
bSig >>= 2;
 3018     }
 3019     do {
alternateASig = aSig;
 3021         ++q;
aSig -= bSig;
 3023     } while ( 0 <= (sbits64) aSig );
sigMean = aSig + alternateASig;
 3025     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
aSig = alternateASig;
 3027     }
zSign = ( (sbits64) aSig < 0 );
 3029     if ( zSign ) aSig = - aSig;
 3030     return normalizeRoundAndPackFloat64( aSign ^ zSign, bExp, aSig );
 3032 }
 3034 /*
 3035 -------------------------------------------------------------------------------
 3036 Returns the square root of the double-precision floating-point value `a'.
 3037 The operation is performed according to the IEC/IEEE Standard for Binary
 3038 Floating-Point Arithmetic.
 3039 -------------------------------------------------------------------------------
 3040 */
float64 float64_sqrt( float64 a )
 3042 {
flag aSign;
int16 aExp, zExp;
bits64 aSig, zSig, doubleZSig;
bits64 rem0, rem1, term0, term1;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 3051     if ( aExp == 0x7FF ) {
 3052         if ( aSig ) return propagateFloat64NaN( a, a );
 3053         if ( ! aSign ) return a;
float_raise( float_flag_invalid );
 3055         return float64_default_nan;
 3056     }
 3057     if ( aSign ) {
 3058         if ( ( aExp | aSig ) == 0 ) return a;
float_raise( float_flag_invalid );
 3060         return float64_default_nan;
 3061     }
 3062     if ( aExp == 0 ) {
 3063         if ( aSig == 0 ) return 0;
normalizeFloat64Subnormal( aSig, &aExp, &aSig );
 3065     }
zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
aSig |= LIT64( 0x0010000000000000 );
zSig = estimateSqrt32( aExp, aSig>>21 );
aSig <<= 9 - ( aExp & 1 );
zSig = estimateDiv128To64( aSig, 0, zSig<<32 ) + ( zSig<<30 );
 3071     if ( ( zSig & 0x1FF ) <= 5 ) {
doubleZSig = zSig<<1;
mul64To128( zSig, zSig, &term0, &term1 );
sub128( aSig, 0, term0, term1, &rem0, &rem1 );
 3075         while ( (sbits64) rem0 < 0 ) {
 3076             --zSig;
doubleZSig -= 2;
add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
 3079         }
zSig |= ( ( rem0 | rem1 ) != 0 );
 3081     }
 3082     return roundAndPackFloat64( 0, zExp, zSig );
 3084 }
 3085 #endif
 3087 /*
 3088 -------------------------------------------------------------------------------
 3089 Returns 1 if the double-precision floating-point value `a' is equal to the
 3090 corresponding value `b', and 0 otherwise.  The comparison is performed
 3091 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3092 -------------------------------------------------------------------------------
 3093 */
flag float64_eq( float64 a, float64 b )
 3095 {
 3097     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
 3098          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
 3099        ) {
 3100         if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 3102         }
 3103         return 0;
 3104     }
 3105     return ( a == b ) ||
 3106         ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) == 0 );
 3108 }
 3110 /*
 3111 -------------------------------------------------------------------------------
 3112 Returns 1 if the double-precision floating-point value `a' is less than or
 3113 equal to the corresponding value `b', and 0 otherwise.  The comparison is
 3114 performed according to the IEC/IEEE Standard for Binary Floating-Point
 3115 Arithmetic.
 3116 -------------------------------------------------------------------------------
 3117 */
flag float64_le( float64 a, float64 b )
 3119 {
flag aSign, bSign;
 3122     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
 3123          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
 3124        ) {
float_raise( float_flag_invalid );
 3126         return 0;
 3127     }
aSign = extractFloat64Sign( a );
bSign = extractFloat64Sign( b );
 3130     if ( aSign != bSign )
 3131         return aSign ||
 3132             ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) ==
 3133               0 );
 3134     return ( a == b ) ||
 3135         ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) );
 3137 }
 3139 /*
 3140 -------------------------------------------------------------------------------
 3141 Returns 1 if the double-precision floating-point value `a' is less than
 3142 the corresponding value `b', and 0 otherwise.  The comparison is performed
 3143 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3144 -------------------------------------------------------------------------------
 3145 */
flag float64_lt( float64 a, float64 b )
 3147 {
flag aSign, bSign;
 3150     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
 3151          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
 3152        ) {
float_raise( float_flag_invalid );
 3154         return 0;
 3155     }
aSign = extractFloat64Sign( a );
bSign = extractFloat64Sign( b );
 3158     if ( aSign != bSign )
 3159         return aSign &&
 3160             ( (bits64) ( ( FLOAT64_DEMANGLE(a) | FLOAT64_DEMANGLE(b) )<<1 ) !=
 3161               0 );
 3162     return ( a != b ) &&
 3163         ( aSign ^ ( FLOAT64_DEMANGLE(a) < FLOAT64_DEMANGLE(b) ) );
 3165 }
 3167 #ifndef SOFTFLOAT_FOR_GCC
 3168 /*
 3169 -------------------------------------------------------------------------------
 3170 Returns 1 if the double-precision floating-point value `a' is equal to the
 3171 corresponding value `b', and 0 otherwise.  The invalid exception is raised
 3172 if either operand is a NaN.  Otherwise, the comparison is performed
 3173 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3174 -------------------------------------------------------------------------------
 3175 */
flag float64_eq_signaling( float64 a, float64 b )
 3177 {
 3179     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
 3180          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
 3181        ) {
float_raise( float_flag_invalid );
 3183         return 0;
 3184     }
 3185     return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
 3187 }
 3189 /*
 3190 -------------------------------------------------------------------------------
 3191 Returns 1 if the double-precision floating-point value `a' is less than or
 3192 equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
 3193 cause an exception.  Otherwise, the comparison is performed according to the
 3194 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3195 -------------------------------------------------------------------------------
 3196 */
flag float64_le_quiet( float64 a, float64 b )
 3198 {
flag aSign, bSign;
 3201     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
 3202          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
 3203        ) {
 3204         if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 3206         }
 3207         return 0;
 3208     }
aSign = extractFloat64Sign( a );
bSign = extractFloat64Sign( b );
 3211     if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
 3212     return ( a == b ) || ( aSign ^ ( a < b ) );
 3214 }
 3216 /*
 3217 -------------------------------------------------------------------------------
 3218 Returns 1 if the double-precision floating-point value `a' is less than
 3219 the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
 3220 exception.  Otherwise, the comparison is performed according to the IEC/IEEE
 3221 Standard for Binary Floating-Point Arithmetic.
 3222 -------------------------------------------------------------------------------
 3223 */
flag float64_lt_quiet( float64 a, float64 b )
 3225 {
flag aSign, bSign;
 3228     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
 3229          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
 3230        ) {
 3231         if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 3233         }
 3234         return 0;
 3235     }
aSign = extractFloat64Sign( a );
bSign = extractFloat64Sign( b );
 3238     if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
 3239     return ( a != b ) && ( aSign ^ ( a < b ) );
 3241 }
 3242 #endif
 3244 #ifdef FLOATX80
 3246 /*
 3247 -------------------------------------------------------------------------------
 3248 Returns the result of converting the extended double-precision floating-
 3249 point value `a' to the 32-bit two's complement integer format.  The
 3250 conversion is performed according to the IEC/IEEE Standard for Binary
 3251 Floating-Point Arithmetic---which means in particular that the conversion
 3252 is rounded according to the current rounding mode.  If `a' is a NaN, the
 3253 largest positive integer is returned.  Otherwise, if the conversion
 3254 overflows, the largest integer with the same sign as `a' is returned.
 3255 -------------------------------------------------------------------------------
 3256 */
int32 floatx80_to_int32( floatx80 a )
 3258 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
 3266     if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
shiftCount = 0x4037 - aExp;
 3268     if ( shiftCount <= 0 ) shiftCount = 1;
shift64RightJamming( aSig, shiftCount, &aSig );
 3270     return roundAndPackInt32( aSign, aSig );
 3272 }
 3274 /*
 3275 -------------------------------------------------------------------------------
 3276 Returns the result of converting the extended double-precision floating-
 3277 point value `a' to the 32-bit two's complement integer format.  The
 3278 conversion is performed according to the IEC/IEEE Standard for Binary
 3279 Floating-Point Arithmetic, except that the conversion is always rounded
 3280 toward zero.  If `a' is a NaN, the largest positive integer is returned.
 3281 Otherwise, if the conversion overflows, the largest integer with the same
 3282 sign as `a' is returned.
 3283 -------------------------------------------------------------------------------
 3284 */
int32 floatx80_to_int32_round_to_zero( floatx80 a )
 3286 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig, savedASig;
int32 z;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
 3295     if ( 0x401E < aExp ) {
 3296         if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
 3297         goto invalid;
 3298     }
 3299     else if ( aExp < 0x3FFF ) {
 3300         if ( aExp || aSig ) float_set_inexact();
 3301         return 0;
 3302     }
shiftCount = 0x403E - aExp;
savedASig = aSig;
aSig >>= shiftCount;
z = aSig;
 3307     if ( aSign ) z = - z;
 3308     if ( ( z < 0 ) ^ aSign ) {
invalid:
float_raise( float_flag_invalid );
 3311         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
 3312     }
 3313     if ( ( aSig<<shiftCount ) != savedASig ) {
float_set_inexact();
 3315     }
 3316     return z;
 3318 }
 3320 /*
 3321 -------------------------------------------------------------------------------
 3322 Returns the result of converting the extended double-precision floating-
 3323 point value `a' to the 64-bit two's complement integer format.  The
 3324 conversion is performed according to the IEC/IEEE Standard for Binary
 3325 Floating-Point Arithmetic---which means in particular that the conversion
 3326 is rounded according to the current rounding mode.  If `a' is a NaN,
 3327 the largest positive integer is returned.  Otherwise, if the conversion
 3328 overflows, the largest integer with the same sign as `a' is returned.
 3329 -------------------------------------------------------------------------------
 3330 */
int64 floatx80_to_int64( floatx80 a )
 3332 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig, aSigExtra;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
shiftCount = 0x403E - aExp;
 3341     if ( shiftCount <= 0 ) {
 3342         if ( shiftCount ) {
float_raise( float_flag_invalid );
 3344             if (    ! aSign
 3345                  || (    ( aExp == 0x7FFF )
 3346                       && ( aSig != LIT64( 0x8000000000000000 ) ) )
 3347                ) {
 3348                 return LIT64( 0x7FFFFFFFFFFFFFFF );
 3349             }
 3350             return (sbits64) LIT64( 0x8000000000000000 );
 3351         }
aSigExtra = 0;
 3353     }
 3354     else {
shift64ExtraRightJamming( aSig, 0, shiftCount, &aSig, &aSigExtra );
 3356     }
 3357     return roundAndPackInt64( aSign, aSig, aSigExtra );
 3359 }
 3361 /*
 3362 -------------------------------------------------------------------------------
 3363 Returns the result of converting the extended double-precision floating-
 3364 point value `a' to the 64-bit two's complement integer format.  The
 3365 conversion is performed according to the IEC/IEEE Standard for Binary
 3366 Floating-Point Arithmetic, except that the conversion is always rounded
 3367 toward zero.  If `a' is a NaN, the largest positive integer is returned.
 3368 Otherwise, if the conversion overflows, the largest integer with the same
 3369 sign as `a' is returned.
 3370 -------------------------------------------------------------------------------
 3371 */
int64 floatx80_to_int64_round_to_zero( floatx80 a )
 3373 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig;
int64 z;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
shiftCount = aExp - 0x403E;
 3383     if ( 0 <= shiftCount ) {
aSig &= LIT64( 0x7FFFFFFFFFFFFFFF );
 3385         if ( ( a.high != 0xC03E ) || aSig ) {
float_raise( float_flag_invalid );
 3387             if ( ! aSign || ( ( aExp == 0x7FFF ) && aSig ) ) {
 3388                 return LIT64( 0x7FFFFFFFFFFFFFFF );
 3389             }
 3390         }
 3391         return (sbits64) LIT64( 0x8000000000000000 );
 3392     }
 3393     else if ( aExp < 0x3FFF ) {
 3394         if ( aExp | aSig ) float_set_inexact();
 3395         return 0;
 3396     }
z = aSig>>( - shiftCount );
 3398     if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
float_set_inexact();
 3400     }
 3401     if ( aSign ) z = - z;
 3402     return z;
 3404 }
 3406 /*
 3407 -------------------------------------------------------------------------------
 3408 Returns the result of converting the extended double-precision floating-
 3409 point value `a' to the single-precision floating-point format.  The
 3410 conversion is performed according to the IEC/IEEE Standard for Binary
 3411 Floating-Point Arithmetic.
 3412 -------------------------------------------------------------------------------
 3413 */
float32 floatx80_to_float32( floatx80 a )
 3415 {
flag aSign;
int32 aExp;
bits64 aSig;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
 3423     if ( aExp == 0x7FFF ) {
 3424         if ( (bits64) ( aSig<<1 ) ) {
 3425             return commonNaNToFloat32( floatx80ToCommonNaN( a ) );
 3426         }
 3427         return packFloat32( aSign, 0xFF, 0 );
 3428     }
shift64RightJamming( aSig, 33, &aSig );
 3430     if ( aExp || aSig ) aExp -= 0x3F81;
 3431     return roundAndPackFloat32( aSign, aExp, aSig );
 3433 }
 3435 /*
 3436 -------------------------------------------------------------------------------
 3437 Returns the result of converting the extended double-precision floating-
 3438 point value `a' to the double-precision floating-point format.  The
 3439 conversion is performed according to the IEC/IEEE Standard for Binary
 3440 Floating-Point Arithmetic.
 3441 -------------------------------------------------------------------------------
 3442 */
float64 floatx80_to_float64( floatx80 a )
 3444 {
flag aSign;
int32 aExp;
bits64 aSig, zSig;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
 3452     if ( aExp == 0x7FFF ) {
 3453         if ( (bits64) ( aSig<<1 ) ) {
 3454             return commonNaNToFloat64( floatx80ToCommonNaN( a ) );
 3455         }
 3456         return packFloat64( aSign, 0x7FF, 0 );
 3457     }
shift64RightJamming( aSig, 1, &zSig );
 3459     if ( aExp || aSig ) aExp -= 0x3C01;
 3460     return roundAndPackFloat64( aSign, aExp, zSig );
 3462 }
 3464 #ifdef FLOAT128
 3466 /*
 3467 -------------------------------------------------------------------------------
 3468 Returns the result of converting the extended double-precision floating-
 3469 point value `a' to the quadruple-precision floating-point format.  The
 3470 conversion is performed according to the IEC/IEEE Standard for Binary
 3471 Floating-Point Arithmetic.
 3472 -------------------------------------------------------------------------------
 3473 */
float128 floatx80_to_float128( floatx80 a )
 3475 {
flag aSign;
int16 aExp;
bits64 aSig, zSig0, zSig1;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
 3483     if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
 3484         return commonNaNToFloat128( floatx80ToCommonNaN( a ) );
 3485     }
shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
 3487     return packFloat128( aSign, aExp, zSig0, zSig1 );
 3489 }
 3491 #endif
 3493 /*
 3494 -------------------------------------------------------------------------------
 3495 Rounds the extended double-precision floating-point value `a' to an integer,
 3496 and returns the result as an extended quadruple-precision floating-point
 3497 value.  The operation is performed according to the IEC/IEEE Standard for
 3498 Binary Floating-Point Arithmetic.
 3499 -------------------------------------------------------------------------------
 3500 */
floatx80 floatx80_round_to_int( floatx80 a )
 3502 {
flag aSign;
int32 aExp;
bits64 lastBitMask, roundBitsMask;
int8 roundingMode;
floatx80 z;
aExp = extractFloatx80Exp( a );
 3510     if ( 0x403E <= aExp ) {
 3511         if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
 3512             return propagateFloatx80NaN( a, a );
 3513         }
 3514         return a;
 3515     }
 3516     if ( aExp < 0x3FFF ) {
 3517         if (    ( aExp == 0 )
 3518              && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
 3519             return a;
 3520         }
float_set_inexact();
aSign = extractFloatx80Sign( a );
 3523         switch ( float_rounding_mode() ) {
 3524          case float_round_nearest_even:
 3525             if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
 3526                ) {
 3527                 return
packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
 3529             }
 3530             break;
 3531          case float_round_down:
 3532             return
aSign ?
packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
 3535                 : packFloatx80( 0, 0, 0 );
 3536          case float_round_up:
 3537             return
aSign ? packFloatx80( 1, 0, 0 )
 3539                 : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
 3540         }
 3541         return packFloatx80( aSign, 0, 0 );
 3542     }
lastBitMask = 1;
lastBitMask <<= 0x403E - aExp;
roundBitsMask = lastBitMask - 1;
z = a;
roundingMode = float_rounding_mode();
 3548     if ( roundingMode == float_round_nearest_even ) {
z.low += lastBitMask>>1;
 3550         if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
 3551     }
 3552     else if ( roundingMode != float_round_to_zero ) {
 3553         if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) {
z.low += roundBitsMask;
 3555         }
 3556     }
z.low &= ~ roundBitsMask;
 3558     if ( z.low == 0 ) {
 3559         ++z.high;
z.low = LIT64( 0x8000000000000000 );
 3561     }
 3562     if ( z.low != a.low ) float_set_inexact();
 3563     return z;
 3565 }
 3567 /*
 3568 -------------------------------------------------------------------------------
 3569 Returns the result of adding the absolute values of the extended double-
 3570 precision floating-point values `a' and `b'.  If `zSign' is 1, the sum is
 3571 negated before being returned.  `zSign' is ignored if the result is a NaN.
 3572 The addition is performed according to the IEC/IEEE Standard for Binary
 3573 Floating-Point Arithmetic.
 3574 -------------------------------------------------------------------------------
 3575 */
 3576 static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign )
 3577 {
int32 aExp, bExp, zExp;
bits64 aSig, bSig, zSig0, zSig1;
int32 expDiff;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
bSig = extractFloatx80Frac( b );
bExp = extractFloatx80Exp( b );
expDiff = aExp - bExp;
 3587     if ( 0 < expDiff ) {
 3588         if ( aExp == 0x7FFF ) {
 3589             if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3590             return a;
 3591         }
 3592         if ( bExp == 0 ) --expDiff;
shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
zExp = aExp;
 3595     }
 3596     else if ( expDiff < 0 ) {
 3597         if ( bExp == 0x7FFF ) {
 3598             if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3599             return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 3600         }
 3601         if ( aExp == 0 ) ++expDiff;
shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
zExp = bExp;
 3604     }
 3605     else {
 3606         if ( aExp == 0x7FFF ) {
 3607             if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
 3608                 return propagateFloatx80NaN( a, b );
 3609             }
 3610             return a;
 3611         }
zSig1 = 0;
zSig0 = aSig + bSig;
 3614         if ( aExp == 0 ) {
normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
 3616             goto roundAndPack;
 3617         }
zExp = aExp;
 3619         goto shiftRight1;
 3620     }
zSig0 = aSig + bSig;
 3622     if ( (sbits64) zSig0 < 0 ) goto roundAndPack;
shiftRight1:
shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
zSig0 |= LIT64( 0x8000000000000000 );
 3626     ++zExp;
roundAndPack:
 3628     return
roundAndPackFloatx80(
floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
 3632 }
 3634 /*
 3635 -------------------------------------------------------------------------------
 3636 Returns the result of subtracting the absolute values of the extended
 3637 double-precision floating-point values `a' and `b'.  If `zSign' is 1, the
 3638 difference is negated before being returned.  `zSign' is ignored if the
 3639 result is a NaN.  The subtraction is performed according to the IEC/IEEE
 3640 Standard for Binary Floating-Point Arithmetic.
 3641 -------------------------------------------------------------------------------
 3642 */
 3643 static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign )
 3644 {
int32 aExp, bExp, zExp;
bits64 aSig, bSig, zSig0, zSig1;
int32 expDiff;
floatx80 z;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
bSig = extractFloatx80Frac( b );
bExp = extractFloatx80Exp( b );
expDiff = aExp - bExp;
 3655     if ( 0 < expDiff ) goto aExpBigger;
 3656     if ( expDiff < 0 ) goto bExpBigger;
 3657     if ( aExp == 0x7FFF ) {
 3658         if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
 3659             return propagateFloatx80NaN( a, b );
 3660         }
float_raise( float_flag_invalid );
z.low = floatx80_default_nan_low;
z.high = floatx80_default_nan_high;
 3664         return z;
 3665     }
 3666     if ( aExp == 0 ) {
aExp = 1;
bExp = 1;
 3669     }
zSig1 = 0;
 3671     if ( bSig < aSig ) goto aBigger;
 3672     if ( aSig < bSig ) goto bBigger;
 3673     return packFloatx80( float_rounding_mode() == float_round_down, 0, 0 );
bExpBigger:
 3675     if ( bExp == 0x7FFF ) {
 3676         if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3677         return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
 3678     }
 3679     if ( aExp == 0 ) ++expDiff;
shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
bBigger:
sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
zExp = bExp;
zSign ^= 1;
 3685     goto normalizeRoundAndPack;
aExpBigger:
 3687     if ( aExp == 0x7FFF ) {
 3688         if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3689         return a;
 3690     }
 3691     if ( bExp == 0 ) --expDiff;
shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
aBigger:
sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
zExp = aExp;
normalizeRoundAndPack:
 3697     return
normalizeRoundAndPackFloatx80(
floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
 3701 }
 3703 /*
 3704 -------------------------------------------------------------------------------
 3705 Returns the result of adding the extended double-precision floating-point
 3706 values `a' and `b'.  The operation is performed according to the IEC/IEEE
 3707 Standard for Binary Floating-Point Arithmetic.
 3708 -------------------------------------------------------------------------------
 3709 */
floatx80 floatx80_add( floatx80 a, floatx80 b )
 3711 {
flag aSign, bSign;
aSign = extractFloatx80Sign( a );
bSign = extractFloatx80Sign( b );
 3716     if ( aSign == bSign ) {
 3717         return addFloatx80Sigs( a, b, aSign );
 3718     }
 3719     else {
 3720         return subFloatx80Sigs( a, b, aSign );
 3721     }
 3723 }
 3725 /*
 3726 -------------------------------------------------------------------------------
 3727 Returns the result of subtracting the extended double-precision floating-
 3728 point values `a' and `b'.  The operation is performed according to the
 3729 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3730 -------------------------------------------------------------------------------
 3731 */
floatx80 floatx80_sub( floatx80 a, floatx80 b )
 3733 {
flag aSign, bSign;
aSign = extractFloatx80Sign( a );
bSign = extractFloatx80Sign( b );
 3738     if ( aSign == bSign ) {
 3739         return subFloatx80Sigs( a, b, aSign );
 3740     }
 3741     else {
 3742         return addFloatx80Sigs( a, b, aSign );
 3743     }
 3745 }
 3747 /*
 3748 -------------------------------------------------------------------------------
 3749 Returns the result of multiplying the extended double-precision floating-
 3750 point values `a' and `b'.  The operation is performed according to the
 3751 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3752 -------------------------------------------------------------------------------
 3753 */
floatx80 floatx80_mul( floatx80 a, floatx80 b )
 3755 {
flag aSign, bSign, zSign;
int32 aExp, bExp, zExp;
bits64 aSig, bSig, zSig0, zSig1;
floatx80 z;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
bSig = extractFloatx80Frac( b );
bExp = extractFloatx80Exp( b );
bSign = extractFloatx80Sign( b );
zSign = aSign ^ bSign;
 3768     if ( aExp == 0x7FFF ) {
 3769         if (    (bits64) ( aSig<<1 )
 3770              || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
 3771             return propagateFloatx80NaN( a, b );
 3772         }
 3773         if ( ( bExp | bSig ) == 0 ) goto invalid;
 3774         return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 3775     }
 3776     if ( bExp == 0x7FFF ) {
 3777         if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3778         if ( ( aExp | aSig ) == 0 ) {
invalid:
float_raise( float_flag_invalid );
z.low = floatx80_default_nan_low;
z.high = floatx80_default_nan_high;
 3783             return z;
 3784         }
 3785         return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 3786     }
 3787     if ( aExp == 0 ) {
 3788         if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
 3790     }
 3791     if ( bExp == 0 ) {
 3792         if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
 3794     }
zExp = aExp + bExp - 0x3FFE;
mul64To128( aSig, bSig, &zSig0, &zSig1 );
 3797     if ( 0 < (sbits64) zSig0 ) {
shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
 3799         --zExp;
 3800     }
 3801     return
roundAndPackFloatx80(
floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
 3805 }
 3807 /*
 3808 -------------------------------------------------------------------------------
 3809 Returns the result of dividing the extended double-precision floating-point
 3810 value `a' by the corresponding value `b'.  The operation is performed
 3811 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3812 -------------------------------------------------------------------------------
 3813 */
floatx80 floatx80_div( floatx80 a, floatx80 b )
 3815 {
flag aSign, bSign, zSign;
int32 aExp, bExp, zExp;
bits64 aSig, bSig, zSig0, zSig1;
bits64 rem0, rem1, rem2, term0, term1, term2;
floatx80 z;
aSig = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
bSig = extractFloatx80Frac( b );
bExp = extractFloatx80Exp( b );
bSign = extractFloatx80Sign( b );
zSign = aSign ^ bSign;
 3829     if ( aExp == 0x7FFF ) {
 3830         if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3831         if ( bExp == 0x7FFF ) {
 3832             if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3833             goto invalid;
 3834         }
 3835         return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 3836     }
 3837     if ( bExp == 0x7FFF ) {
 3838         if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3839         return packFloatx80( zSign, 0, 0 );
 3840     }
 3841     if ( bExp == 0 ) {
 3842         if ( bSig == 0 ) {
 3843             if ( ( aExp | aSig ) == 0 ) {
invalid:
float_raise( float_flag_invalid );
z.low = floatx80_default_nan_low;
z.high = floatx80_default_nan_high;
 3848                 return z;
 3849             }
float_raise( float_flag_divbyzero );
 3851             return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 3852         }
normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
 3854     }
 3855     if ( aExp == 0 ) {
 3856         if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
 3858     }
zExp = aExp - bExp + 0x3FFE;
rem1 = 0;
 3861     if ( bSig <= aSig ) {
shift128Right( aSig, 0, 1, &aSig, &rem1 );
 3863         ++zExp;
 3864     }
zSig0 = estimateDiv128To64( aSig, rem1, bSig );
mul64To128( bSig, zSig0, &term0, &term1 );
sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
 3868     while ( (sbits64) rem0 < 0 ) {
 3869         --zSig0;
add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
 3871     }
zSig1 = estimateDiv128To64( rem1, 0, bSig );
 3873     if ( (bits64) ( zSig1<<1 ) <= 8 ) {
mul64To128( bSig, zSig1, &term1, &term2 );
sub128( rem1, 0, term1, term2, &rem1, &rem2 );
 3876         while ( (sbits64) rem1 < 0 ) {
 3877             --zSig1;
add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
 3879         }
zSig1 |= ( ( rem1 | rem2 ) != 0 );
 3881     }
 3882     return
roundAndPackFloatx80(
floatx80_rounding_precision, zSign, zExp, zSig0, zSig1 );
 3886 }
 3888 /*
 3889 -------------------------------------------------------------------------------
 3890 Returns the remainder of the extended double-precision floating-point value
 3891 `a' with respect to the corresponding value `b'.  The operation is performed
 3892 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 3893 -------------------------------------------------------------------------------
 3894 */
floatx80 floatx80_rem( floatx80 a, floatx80 b )
 3896 {
flag aSign, bSign, zSign;
int32 aExp, bExp, expDiff;
bits64 aSig0, aSig1, bSig;
bits64 q, term0, term1, alternateASig0, alternateASig1;
floatx80 z;
aSig0 = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
bSig = extractFloatx80Frac( b );
bExp = extractFloatx80Exp( b );
bSign = extractFloatx80Sign( b );
 3909     if ( aExp == 0x7FFF ) {
 3910         if (    (bits64) ( aSig0<<1 )
 3911              || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
 3912             return propagateFloatx80NaN( a, b );
 3913         }
 3914         goto invalid;
 3915     }
 3916     if ( bExp == 0x7FFF ) {
 3917         if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
 3918         return a;
 3919     }
 3920     if ( bExp == 0 ) {
 3921         if ( bSig == 0 ) {
invalid:
float_raise( float_flag_invalid );
z.low = floatx80_default_nan_low;
z.high = floatx80_default_nan_high;
 3926             return z;
 3927         }
normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
 3929     }
 3930     if ( aExp == 0 ) {
 3931         if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
 3933     }
bSig |= LIT64( 0x8000000000000000 );
zSign = aSign;
expDiff = aExp - bExp;
aSig1 = 0;
 3938     if ( expDiff < 0 ) {
 3939         if ( expDiff < -1 ) return a;
shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
expDiff = 0;
 3942     }
q = ( bSig <= aSig0 );
 3944     if ( q ) aSig0 -= bSig;
expDiff -= 64;
 3946     while ( 0 < expDiff ) {
q = estimateDiv128To64( aSig0, aSig1, bSig );
q = ( 2 < q ) ? q - 2 : 0;
mul64To128( bSig, q, &term0, &term1 );
sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
expDiff -= 62;
 3953     }
expDiff += 64;
 3955     if ( 0 < expDiff ) {
q = estimateDiv128To64( aSig0, aSig1, bSig );
q = ( 2 < q ) ? q - 2 : 0;
q >>= 64 - expDiff;
mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
 3962         while ( le128( term0, term1, aSig0, aSig1 ) ) {
 3963             ++q;
sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
 3965         }
 3966     }
 3967     else {
term1 = 0;
term0 = bSig;
 3970     }
sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
 3972     if (    lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
 3973          || (    eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
 3974               && ( q & 1 ) )
 3975        ) {
aSig0 = alternateASig0;
aSig1 = alternateASig1;
zSign = ! zSign;
 3979     }
 3980     return
normalizeRoundAndPackFloatx80(
 3982             80, zSign, bExp + expDiff, aSig0, aSig1 );
 3984 }
 3986 /*
 3987 -------------------------------------------------------------------------------
 3988 Returns the square root of the extended double-precision floating-point
 3989 value `a'.  The operation is performed according to the IEC/IEEE Standard
 3990 for Binary Floating-Point Arithmetic.
 3991 -------------------------------------------------------------------------------
 3992 */
floatx80 floatx80_sqrt( floatx80 a )
 3994 {
flag aSign;
int32 aExp, zExp;
bits64 aSig0, aSig1, zSig0, zSig1, doubleZSig0;
bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
floatx80 z;
aSig0 = extractFloatx80Frac( a );
aExp = extractFloatx80Exp( a );
aSign = extractFloatx80Sign( a );
 4004     if ( aExp == 0x7FFF ) {
 4005         if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a );
 4006         if ( ! aSign ) return a;
 4007         goto invalid;
 4008     }
 4009     if ( aSign ) {
 4010         if ( ( aExp | aSig0 ) == 0 ) return a;
invalid:
float_raise( float_flag_invalid );
z.low = floatx80_default_nan_low;
z.high = floatx80_default_nan_high;
 4015         return z;
 4016     }
 4017     if ( aExp == 0 ) {
 4018         if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
 4020     }
zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
zSig0 = estimateSqrt32( aExp, aSig0>>32 );
shift128Right( aSig0, 0, 2 + ( aExp & 1 ), &aSig0, &aSig1 );
zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
doubleZSig0 = zSig0<<1;
mul64To128( zSig0, zSig0, &term0, &term1 );
sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
 4028     while ( (sbits64) rem0 < 0 ) {
 4029         --zSig0;
doubleZSig0 -= 2;
add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
 4032     }
zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
 4034     if ( ( zSig1 & LIT64( 0x3FFFFFFFFFFFFFFF ) ) <= 5 ) {
 4035         if ( zSig1 == 0 ) zSig1 = 1;
mul64To128( doubleZSig0, zSig1, &term1, &term2 );
sub128( rem1, 0, term1, term2, &rem1, &rem2 );
mul64To128( zSig1, zSig1, &term2, &term3 );
sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
 4040         while ( (sbits64) rem1 < 0 ) {
 4041             --zSig1;
shortShift128Left( 0, zSig1, 1, &term2, &term3 );
term3 |= 1;
term2 |= doubleZSig0;
add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
 4046         }
zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
 4048     }
shortShift128Left( 0, zSig1, 1, &zSig0, &zSig1 );
zSig0 |= doubleZSig0;
 4051     return
roundAndPackFloatx80(
floatx80_rounding_precision, 0, zExp, zSig0, zSig1 );
 4055 }
 4057 /*
 4058 -------------------------------------------------------------------------------
 4059 Returns 1 if the extended double-precision floating-point value `a' is
 4060 equal to the corresponding value `b', and 0 otherwise.  The comparison is
 4061 performed according to the IEC/IEEE Standard for Binary Floating-Point
 4062 Arithmetic.
 4063 -------------------------------------------------------------------------------
 4064 */
flag floatx80_eq( floatx80 a, floatx80 b )
 4066 {
 4068     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
 4069               && (bits64) ( extractFloatx80Frac( a )<<1 ) )
 4070          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
 4071               && (bits64) ( extractFloatx80Frac( b )<<1 ) )
 4072        ) {
 4073         if (    floatx80_is_signaling_nan( a )
 4074              || floatx80_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 4076         }
 4077         return 0;
 4078     }
 4079     return
 4080            ( a.low == b.low )
 4081         && (    ( a.high == b.high )
 4082              || (    ( a.low == 0 )
 4083                   && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
 4084            );
 4086 }
 4088 /*
 4089 -------------------------------------------------------------------------------
 4090 Returns 1 if the extended double-precision floating-point value `a' is
 4091 less than or equal to the corresponding value `b', and 0 otherwise.  The
 4092 comparison is performed according to the IEC/IEEE Standard for Binary
 4093 Floating-Point Arithmetic.
 4094 -------------------------------------------------------------------------------
 4095 */
flag floatx80_le( floatx80 a, floatx80 b )
 4097 {
flag aSign, bSign;
 4100     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
 4101               && (bits64) ( extractFloatx80Frac( a )<<1 ) )
 4102          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
 4103               && (bits64) ( extractFloatx80Frac( b )<<1 ) )
 4104        ) {
float_raise( float_flag_invalid );
 4106         return 0;
 4107     }
aSign = extractFloatx80Sign( a );
bSign = extractFloatx80Sign( b );
 4110     if ( aSign != bSign ) {
 4111         return
aSign
 4113             || (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 4114                  == 0 );
 4115     }
 4116     return
aSign ? le128( b.high, b.low, a.high, a.low )
 4118         : le128( a.high, a.low, b.high, b.low );
 4120 }
 4122 /*
 4123 -------------------------------------------------------------------------------
 4124 Returns 1 if the extended double-precision floating-point value `a' is
 4125 less than the corresponding value `b', and 0 otherwise.  The comparison
 4126 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4127 Arithmetic.
 4128 -------------------------------------------------------------------------------
 4129 */
flag floatx80_lt( floatx80 a, floatx80 b )
 4131 {
flag aSign, bSign;
 4134     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
 4135               && (bits64) ( extractFloatx80Frac( a )<<1 ) )
 4136          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
 4137               && (bits64) ( extractFloatx80Frac( b )<<1 ) )
 4138        ) {
float_raise( float_flag_invalid );
 4140         return 0;
 4141     }
aSign = extractFloatx80Sign( a );
bSign = extractFloatx80Sign( b );
 4144     if ( aSign != bSign ) {
 4145         return
aSign
 4147             && (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 4148                  != 0 );
 4149     }
 4150     return
aSign ? lt128( b.high, b.low, a.high, a.low )
 4152         : lt128( a.high, a.low, b.high, b.low );
 4154 }
 4156 /*
 4157 -------------------------------------------------------------------------------
 4158 Returns 1 if the extended double-precision floating-point value `a' is equal
 4159 to the corresponding value `b', and 0 otherwise.  The invalid exception is
 4160 raised if either operand is a NaN.  Otherwise, the comparison is performed
 4161 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 4162 -------------------------------------------------------------------------------
 4163 */
flag floatx80_eq_signaling( floatx80 a, floatx80 b )
 4165 {
 4167     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
 4168               && (bits64) ( extractFloatx80Frac( a )<<1 ) )
 4169          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
 4170               && (bits64) ( extractFloatx80Frac( b )<<1 ) )
 4171        ) {
float_raise( float_flag_invalid );
 4173         return 0;
 4174     }
 4175     return
 4176            ( a.low == b.low )
 4177         && (    ( a.high == b.high )
 4178              || (    ( a.low == 0 )
 4179                   && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
 4180            );
 4182 }
 4184 /*
 4185 -------------------------------------------------------------------------------
 4186 Returns 1 if the extended double-precision floating-point value `a' is less
 4187 than or equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs
 4188 do not cause an exception.  Otherwise, the comparison is performed according
 4189 to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 4190 -------------------------------------------------------------------------------
 4191 */
flag floatx80_le_quiet( floatx80 a, floatx80 b )
 4193 {
flag aSign, bSign;
 4196     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
 4197               && (bits64) ( extractFloatx80Frac( a )<<1 ) )
 4198          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
 4199               && (bits64) ( extractFloatx80Frac( b )<<1 ) )
 4200        ) {
 4201         if (    floatx80_is_signaling_nan( a )
 4202              || floatx80_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 4204         }
 4205         return 0;
 4206     }
aSign = extractFloatx80Sign( a );
bSign = extractFloatx80Sign( b );
 4209     if ( aSign != bSign ) {
 4210         return
aSign
 4212             || (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 4213                  == 0 );
 4214     }
 4215     return
aSign ? le128( b.high, b.low, a.high, a.low )
 4217         : le128( a.high, a.low, b.high, b.low );
 4219 }
 4221 /*
 4222 -------------------------------------------------------------------------------
 4223 Returns 1 if the extended double-precision floating-point value `a' is less
 4224 than the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause
 4225 an exception.  Otherwise, the comparison is performed according to the
 4226 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 4227 -------------------------------------------------------------------------------
 4228 */
flag floatx80_lt_quiet( floatx80 a, floatx80 b )
 4230 {
flag aSign, bSign;
 4233     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
 4234               && (bits64) ( extractFloatx80Frac( a )<<1 ) )
 4235          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
 4236               && (bits64) ( extractFloatx80Frac( b )<<1 ) )
 4237        ) {
 4238         if (    floatx80_is_signaling_nan( a )
 4239              || floatx80_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 4241         }
 4242         return 0;
 4243     }
aSign = extractFloatx80Sign( a );
bSign = extractFloatx80Sign( b );
 4246     if ( aSign != bSign ) {
 4247         return
aSign
 4249             && (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 4250                  != 0 );
 4251     }
 4252     return
aSign ? lt128( b.high, b.low, a.high, a.low )
 4254         : lt128( a.high, a.low, b.high, b.low );
 4256 }
 4258 #endif
 4260 #ifdef FLOAT128
 4262 /*
 4263 -------------------------------------------------------------------------------
 4264 Returns the result of converting the quadruple-precision floating-point
 4265 value `a' to the 32-bit two's complement integer format.  The conversion
 4266 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4267 Arithmetic---which means in particular that the conversion is rounded
 4268 according to the current rounding mode.  If `a' is a NaN, the largest
 4269 positive integer is returned.  Otherwise, if the conversion overflows, the
 4270 largest integer with the same sign as `a' is returned.
 4271 -------------------------------------------------------------------------------
 4272 */
int32 float128_to_int32( float128 a )
 4274 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig0, aSig1;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 4283     if ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) aSign = 0;
 4284     if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
aSig0 |= ( aSig1 != 0 );
shiftCount = 0x4028 - aExp;
 4287     if ( 0 < shiftCount ) shift64RightJamming( aSig0, shiftCount, &aSig0 );
 4288     return roundAndPackInt32( aSign, aSig0 );
 4290 }
 4292 /*
 4293 -------------------------------------------------------------------------------
 4294 Returns the result of converting the quadruple-precision floating-point
 4295 value `a' to the 32-bit two's complement integer format.  The conversion
 4296 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4297 Arithmetic, except that the conversion is always rounded toward zero.  If
 4298 `a' is a NaN, the largest positive integer is returned.  Otherwise, if the
 4299 conversion overflows, the largest integer with the same sign as `a' is
 4300 returned.
 4301 -------------------------------------------------------------------------------
 4302 */
int32 float128_to_int32_round_to_zero( float128 a )
 4304 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig0, aSig1, savedASig;
int32 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
aSig0 |= ( aSig1 != 0 );
 4315     if ( 0x401E < aExp ) {
 4316         if ( ( aExp == 0x7FFF ) && aSig0 ) aSign = 0;
 4317         goto invalid;
 4318     }
 4319     else if ( aExp < 0x3FFF ) {
 4320         if ( aExp || aSig0 ) float_set_inexact();
 4321         return 0;
 4322     }
aSig0 |= LIT64( 0x0001000000000000 );
shiftCount = 0x402F - aExp;
savedASig = aSig0;
aSig0 >>= shiftCount;
z = aSig0;
 4328     if ( aSign ) z = - z;
 4329     if ( ( z < 0 ) ^ aSign ) {
invalid:
float_raise( float_flag_invalid );
 4332         return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
 4333     }
 4334     if ( ( aSig0<<shiftCount ) != savedASig ) {
float_set_inexact();
 4336     }
 4337     return z;
 4339 }
 4341 /*
 4342 -------------------------------------------------------------------------------
 4343 Returns the result of converting the quadruple-precision floating-point
 4344 value `a' to the 64-bit two's complement integer format.  The conversion
 4345 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4346 Arithmetic---which means in particular that the conversion is rounded
 4347 according to the current rounding mode.  If `a' is a NaN, the largest
 4348 positive integer is returned.  Otherwise, if the conversion overflows, the
 4349 largest integer with the same sign as `a' is returned.
 4350 -------------------------------------------------------------------------------
 4351 */
int64 float128_to_int64( float128 a )
 4353 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig0, aSig1;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 4362     if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
shiftCount = 0x402F - aExp;
 4364     if ( shiftCount <= 0 ) {
 4365         if ( 0x403E < aExp ) {
float_raise( float_flag_invalid );
 4367             if (    ! aSign
 4368                  || (    ( aExp == 0x7FFF )
 4369                       && ( aSig1 || ( aSig0 != LIT64( 0x0001000000000000 ) ) )
 4370                     )
 4371                ) {
 4372                 return LIT64( 0x7FFFFFFFFFFFFFFF );
 4373             }
 4374             return (sbits64) LIT64( 0x8000000000000000 );
 4375         }
shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
 4377     }
 4378     else {
shift64ExtraRightJamming( aSig0, aSig1, shiftCount, &aSig0, &aSig1 );
 4380     }
 4381     return roundAndPackInt64( aSign, aSig0, aSig1 );
 4383 }
 4385 /*
 4386 -------------------------------------------------------------------------------
 4387 Returns the result of converting the quadruple-precision floating-point
 4388 value `a' to the 64-bit two's complement integer format.  The conversion
 4389 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4390 Arithmetic, except that the conversion is always rounded toward zero.
 4391 If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
 4392 the conversion overflows, the largest integer with the same sign as `a' is
 4393 returned.
 4394 -------------------------------------------------------------------------------
 4395 */
int64 float128_to_int64_round_to_zero( float128 a )
 4397 {
flag aSign;
int32 aExp, shiftCount;
bits64 aSig0, aSig1;
int64 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 4407     if ( aExp ) aSig0 |= LIT64( 0x0001000000000000 );
shiftCount = aExp - 0x402F;
 4409     if ( 0 < shiftCount ) {
 4410         if ( 0x403E <= aExp ) {
aSig0 &= LIT64( 0x0000FFFFFFFFFFFF );
 4412             if (    ( a.high == LIT64( 0xC03E000000000000 ) )
 4413                  && ( aSig1 < LIT64( 0x0002000000000000 ) ) ) {
 4414                 if ( aSig1 ) float_set_inexact();
 4415             }
 4416             else {
float_raise( float_flag_invalid );
 4418                 if ( ! aSign || ( ( aExp == 0x7FFF ) && ( aSig0 | aSig1 ) ) ) {
 4419                     return LIT64( 0x7FFFFFFFFFFFFFFF );
 4420                 }
 4421             }
 4422             return (sbits64) LIT64( 0x8000000000000000 );
 4423         }
z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
 4425         if ( (bits64) ( aSig1<<shiftCount ) ) {
float_set_inexact();
 4427         }
 4428     }
 4429     else {
 4430         if ( aExp < 0x3FFF ) {
 4431             if ( aExp | aSig0 | aSig1 ) {
float_set_inexact();
 4433             }
 4434             return 0;
 4435         }
z = aSig0>>( - shiftCount );
 4437         if (    aSig1
 4438              || ( shiftCount && (bits64) ( aSig0<<( shiftCount & 63 ) ) ) ) {
float_set_inexact();
 4440         }
 4441     }
 4442     if ( aSign ) z = - z;
 4443     return z;
 4445 }
 4447 /*
 4448 -------------------------------------------------------------------------------
 4449 Returns the result of converting the quadruple-precision floating-point
 4450 value `a' to the single-precision floating-point format.  The conversion
 4451 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4452 Arithmetic.
 4453 -------------------------------------------------------------------------------
 4454 */
float32 float128_to_float32( float128 a )
 4456 {
flag aSign;
int32 aExp;
bits64 aSig0, aSig1;
bits32 zSig;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 4466     if ( aExp == 0x7FFF ) {
 4467         if ( aSig0 | aSig1 ) {
 4468             return commonNaNToFloat32( float128ToCommonNaN( a ) );
 4469         }
 4470         return packFloat32( aSign, 0xFF, 0 );
 4471     }
aSig0 |= ( aSig1 != 0 );
shift64RightJamming( aSig0, 18, &aSig0 );
zSig = aSig0;
 4475     if ( aExp || zSig ) {
zSig |= 0x40000000;
aExp -= 0x3F81;
 4478     }
 4479     return roundAndPackFloat32( aSign, aExp, zSig );
 4481 }
 4483 /*
 4484 -------------------------------------------------------------------------------
 4485 Returns the result of converting the quadruple-precision floating-point
 4486 value `a' to the double-precision floating-point format.  The conversion
 4487 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 4488 Arithmetic.
 4489 -------------------------------------------------------------------------------
 4490 */
float64 float128_to_float64( float128 a )
 4492 {
flag aSign;
int32 aExp;
bits64 aSig0, aSig1;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 4501     if ( aExp == 0x7FFF ) {
 4502         if ( aSig0 | aSig1 ) {
 4503             return commonNaNToFloat64( float128ToCommonNaN( a ) );
 4504         }
 4505         return packFloat64( aSign, 0x7FF, 0 );
 4506     }
shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
aSig0 |= ( aSig1 != 0 );
 4509     if ( aExp || aSig0 ) {
aSig0 |= LIT64( 0x4000000000000000 );
aExp -= 0x3C01;
 4512     }
 4513     return roundAndPackFloat64( aSign, aExp, aSig0 );
 4515 }
 4517 #ifdef FLOATX80
 4519 /*
 4520 -------------------------------------------------------------------------------
 4521 Returns the result of converting the quadruple-precision floating-point
 4522 value `a' to the extended double-precision floating-point format.  The
 4523 conversion is performed according to the IEC/IEEE Standard for Binary
 4524 Floating-Point Arithmetic.
 4525 -------------------------------------------------------------------------------
 4526 */
floatx80 float128_to_floatx80( float128 a )
 4528 {
flag aSign;
int32 aExp;
bits64 aSig0, aSig1;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 4537     if ( aExp == 0x7FFF ) {
 4538         if ( aSig0 | aSig1 ) {
 4539             return commonNaNToFloatx80( float128ToCommonNaN( a ) );
 4540         }
 4541         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
 4542     }
 4543     if ( aExp == 0 ) {
 4544         if ( ( aSig0 | aSig1 ) == 0 ) return packFloatx80( aSign, 0, 0 );
normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
 4546     }
 4547     else {
aSig0 |= LIT64( 0x0001000000000000 );
 4549     }
shortShift128Left( aSig0, aSig1, 15, &aSig0, &aSig1 );
 4551     return roundAndPackFloatx80( 80, aSign, aExp, aSig0, aSig1 );
 4553 }
 4555 #endif
 4557 /*
 4558 -------------------------------------------------------------------------------
 4559 Rounds the quadruple-precision floating-point value `a' to an integer, and
 4560 returns the result as a quadruple-precision floating-point value.  The
 4561 operation is performed according to the IEC/IEEE Standard for Binary
 4562 Floating-Point Arithmetic.
 4563 -------------------------------------------------------------------------------
 4564 */
float128 float128_round_to_int( float128 a )
 4566 {
flag aSign;
int32 aExp;
bits64 lastBitMask, roundBitsMask;
int8 roundingMode;
float128 z;
aExp = extractFloat128Exp( a );
 4574     if ( 0x402F <= aExp ) {
 4575         if ( 0x406F <= aExp ) {
 4576             if (    ( aExp == 0x7FFF )
 4577                  && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) )
 4578                ) {
 4579                 return propagateFloat128NaN( a, a );
 4580             }
 4581             return a;
 4582         }
lastBitMask = 1;
lastBitMask = ( lastBitMask<<( 0x406E - aExp ) )<<1;
roundBitsMask = lastBitMask - 1;
z = a;
roundingMode = float_rounding_mode();
 4588         if ( roundingMode == float_round_nearest_even ) {
 4589             if ( lastBitMask ) {
add128( z.high, z.low, 0, lastBitMask>>1, &z.high, &z.low );
 4591                 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
 4592             }
 4593             else {
 4594                 if ( (sbits64) z.low < 0 ) {
 4595                     ++z.high;
 4596                     if ( (bits64) ( z.low<<1 ) == 0 ) z.high &= ~1;
 4597                 }
 4598             }
 4599         }
 4600         else if ( roundingMode != float_round_to_zero ) {
 4601             if (   extractFloat128Sign( z )
 4602                  ^ ( roundingMode == float_round_up ) ) {
add128( z.high, z.low, 0, roundBitsMask, &z.high, &z.low );
 4604             }
 4605         }
z.low &= ~ roundBitsMask;
 4607     }
 4608     else {
 4609         if ( aExp < 0x3FFF ) {
 4610             if ( ( ( (bits64) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
float_set_inexact();
aSign = extractFloat128Sign( a );
 4613             switch ( float_rounding_mode() ) {
 4614              case float_round_nearest_even:
 4615                 if (    ( aExp == 0x3FFE )
 4616                      && (   extractFloat128Frac0( a )
 4617                           | extractFloat128Frac1( a ) )
 4618                    ) {
 4619                     return packFloat128( aSign, 0x3FFF, 0, 0 );
 4620                 }
 4621                 break;
 4622              case float_round_down:
 4623                 return
aSign ? packFloat128( 1, 0x3FFF, 0, 0 )
 4625                     : packFloat128( 0, 0, 0, 0 );
 4626              case float_round_up:
 4627                 return
aSign ? packFloat128( 1, 0, 0, 0 )
 4629                     : packFloat128( 0, 0x3FFF, 0, 0 );
 4630             }
 4631             return packFloat128( aSign, 0, 0, 0 );
 4632         }
lastBitMask = 1;
lastBitMask <<= 0x402F - aExp;
roundBitsMask = lastBitMask - 1;
z.low = 0;
z.high = a.high;
roundingMode = float_rounding_mode();
 4639         if ( roundingMode == float_round_nearest_even ) {
z.high += lastBitMask>>1;
 4641             if ( ( ( z.high & roundBitsMask ) | a.low ) == 0 ) {
z.high &= ~ lastBitMask;
 4643             }
 4644         }
 4645         else if ( roundingMode != float_round_to_zero ) {
 4646             if (   extractFloat128Sign( z )
 4647                  ^ ( roundingMode == float_round_up ) ) {
z.high |= ( a.low != 0 );
z.high += roundBitsMask;
 4650             }
 4651         }
z.high &= ~ roundBitsMask;
 4653     }
 4654     if ( ( z.low != a.low ) || ( z.high != a.high ) ) {
float_set_inexact();
 4656     }
 4657     return z;
 4659 }
 4661 /*
 4662 -------------------------------------------------------------------------------
 4663 Returns the result of adding the absolute values of the quadruple-precision
 4664 floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
 4665 before being returned.  `zSign' is ignored if the result is a NaN.
 4666 The addition is performed according to the IEC/IEEE Standard for Binary
 4667 Floating-Point Arithmetic.
 4668 -------------------------------------------------------------------------------
 4669 */
 4670 static float128 addFloat128Sigs( float128 a, float128 b, flag zSign )
 4671 {
int32 aExp, bExp, zExp;
bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
int32 expDiff;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
bSig1 = extractFloat128Frac1( b );
bSig0 = extractFloat128Frac0( b );
bExp = extractFloat128Exp( b );
expDiff = aExp - bExp;
 4683     if ( 0 < expDiff ) {
 4684         if ( aExp == 0x7FFF ) {
 4685             if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b );
 4686             return a;
 4687         }
 4688         if ( bExp == 0 ) {
 4689             --expDiff;
 4690         }
 4691         else {
bSig0 |= LIT64( 0x0001000000000000 );
 4693         }
shift128ExtraRightJamming(
bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
zExp = aExp;
 4697     }
 4698     else if ( expDiff < 0 ) {
 4699         if ( bExp == 0x7FFF ) {
 4700             if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
 4701             return packFloat128( zSign, 0x7FFF, 0, 0 );
 4702         }
 4703         if ( aExp == 0 ) {
 4704             ++expDiff;
 4705         }
 4706         else {
aSig0 |= LIT64( 0x0001000000000000 );
 4708         }
shift128ExtraRightJamming(
aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
zExp = bExp;
 4712     }
 4713     else {
 4714         if ( aExp == 0x7FFF ) {
 4715             if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
 4716                 return propagateFloat128NaN( a, b );
 4717             }
 4718             return a;
 4719         }
add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
 4721         if ( aExp == 0 ) return packFloat128( zSign, 0, zSig0, zSig1 );
zSig2 = 0;
zSig0 |= LIT64( 0x0002000000000000 );
zExp = aExp;
 4725         goto shiftRight1;
 4726     }
aSig0 |= LIT64( 0x0001000000000000 );
add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
 4729     --zExp;
 4730     if ( zSig0 < LIT64( 0x0002000000000000 ) ) goto roundAndPack;
 4731     ++zExp;
shiftRight1:
shift128ExtraRightJamming(
zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
roundAndPack:
 4736     return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
 4738 }
 4740 /*
 4741 -------------------------------------------------------------------------------
 4742 Returns the result of subtracting the absolute values of the quadruple-
 4743 precision floating-point values `a' and `b'.  If `zSign' is 1, the
 4744 difference is negated before being returned.  `zSign' is ignored if the
 4745 result is a NaN.  The subtraction is performed according to the IEC/IEEE
 4746 Standard for Binary Floating-Point Arithmetic.
 4747 -------------------------------------------------------------------------------
 4748 */
 4749 static float128 subFloat128Sigs( float128 a, float128 b, flag zSign )
 4750 {
int32 aExp, bExp, zExp;
bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
int32 expDiff;
float128 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
bSig1 = extractFloat128Frac1( b );
bSig0 = extractFloat128Frac0( b );
bExp = extractFloat128Exp( b );
expDiff = aExp - bExp;
shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
 4765     if ( 0 < expDiff ) goto aExpBigger;
 4766     if ( expDiff < 0 ) goto bExpBigger;
 4767     if ( aExp == 0x7FFF ) {
 4768         if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
 4769             return propagateFloat128NaN( a, b );
 4770         }
float_raise( float_flag_invalid );
z.low = float128_default_nan_low;
z.high = float128_default_nan_high;
 4774         return z;
 4775     }
 4776     if ( aExp == 0 ) {
aExp = 1;
bExp = 1;
 4779     }
 4780     if ( bSig0 < aSig0 ) goto aBigger;
 4781     if ( aSig0 < bSig0 ) goto bBigger;
 4782     if ( bSig1 < aSig1 ) goto aBigger;
 4783     if ( aSig1 < bSig1 ) goto bBigger;
 4784     return packFloat128( float_rounding_mode() == float_round_down, 0, 0, 0 );
bExpBigger:
 4786     if ( bExp == 0x7FFF ) {
 4787         if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
 4788         return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
 4789     }
 4790     if ( aExp == 0 ) {
 4791         ++expDiff;
 4792     }
 4793     else {
aSig0 |= LIT64( 0x4000000000000000 );
 4795     }
shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
bSig0 |= LIT64( 0x4000000000000000 );
bBigger:
sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
zExp = bExp;
zSign ^= 1;
 4802     goto normalizeRoundAndPack;
aExpBigger:
 4804     if ( aExp == 0x7FFF ) {
 4805         if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b );
 4806         return a;
 4807     }
 4808     if ( bExp == 0 ) {
 4809         --expDiff;
 4810     }
 4811     else {
bSig0 |= LIT64( 0x4000000000000000 );
 4813     }
shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
aSig0 |= LIT64( 0x4000000000000000 );
aBigger:
sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
zExp = aExp;
normalizeRoundAndPack:
 4820     --zExp;
 4821     return normalizeRoundAndPackFloat128( zSign, zExp - 14, zSig0, zSig1 );
 4823 }
 4825 /*
 4826 -------------------------------------------------------------------------------
 4827 Returns the result of adding the quadruple-precision floating-point values
 4828 `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
 4829 for Binary Floating-Point Arithmetic.
 4830 -------------------------------------------------------------------------------
 4831 */
float128 float128_add( float128 a, float128 b )
 4833 {
flag aSign, bSign;
aSign = extractFloat128Sign( a );
bSign = extractFloat128Sign( b );
 4838     if ( aSign == bSign ) {
 4839         return addFloat128Sigs( a, b, aSign );
 4840     }
 4841     else {
 4842         return subFloat128Sigs( a, b, aSign );
 4843     }
 4845 }
 4847 /*
 4848 -------------------------------------------------------------------------------
 4849 Returns the result of subtracting the quadruple-precision floating-point
 4850 values `a' and `b'.  The operation is performed according to the IEC/IEEE
 4851 Standard for Binary Floating-Point Arithmetic.
 4852 -------------------------------------------------------------------------------
 4853 */
float128 float128_sub( float128 a, float128 b )
 4855 {
flag aSign, bSign;
aSign = extractFloat128Sign( a );
bSign = extractFloat128Sign( b );
 4860     if ( aSign == bSign ) {
 4861         return subFloat128Sigs( a, b, aSign );
 4862     }
 4863     else {
 4864         return addFloat128Sigs( a, b, aSign );
 4865     }
 4867 }
 4869 /*
 4870 -------------------------------------------------------------------------------
 4871 Returns the result of multiplying the quadruple-precision floating-point
 4872 values `a' and `b'.  The operation is performed according to the IEC/IEEE
 4873 Standard for Binary Floating-Point Arithmetic.
 4874 -------------------------------------------------------------------------------
 4875 */
float128 float128_mul( float128 a, float128 b )
 4877 {
flag aSign, bSign, zSign;
int32 aExp, bExp, zExp;
bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
float128 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
bSig1 = extractFloat128Frac1( b );
bSig0 = extractFloat128Frac0( b );
bExp = extractFloat128Exp( b );
bSign = extractFloat128Sign( b );
zSign = aSign ^ bSign;
 4892     if ( aExp == 0x7FFF ) {
 4893         if (    ( aSig0 | aSig1 )
 4894              || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
 4895             return propagateFloat128NaN( a, b );
 4896         }
 4897         if ( ( bExp | bSig0 | bSig1 ) == 0 ) goto invalid;
 4898         return packFloat128( zSign, 0x7FFF, 0, 0 );
 4899     }
 4900     if ( bExp == 0x7FFF ) {
 4901         if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
 4902         if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
invalid:
float_raise( float_flag_invalid );
z.low = float128_default_nan_low;
z.high = float128_default_nan_high;
 4907             return z;
 4908         }
 4909         return packFloat128( zSign, 0x7FFF, 0, 0 );
 4910     }
 4911     if ( aExp == 0 ) {
 4912         if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
 4914     }
 4915     if ( bExp == 0 ) {
 4916         if ( ( bSig0 | bSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
 4918     }
zExp = aExp + bExp - 0x4000;
aSig0 |= LIT64( 0x0001000000000000 );
shortShift128Left( bSig0, bSig1, 16, &bSig0, &bSig1 );
mul128To256( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1, &zSig2, &zSig3 );
add128( zSig0, zSig1, aSig0, aSig1, &zSig0, &zSig1 );
zSig2 |= ( zSig3 != 0 );
 4925     if ( LIT64( 0x0002000000000000 ) <= zSig0 ) {
shift128ExtraRightJamming(
zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
 4928         ++zExp;
 4929     }
 4930     return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
 4932 }
 4934 /*
 4935 -------------------------------------------------------------------------------
 4936 Returns the result of dividing the quadruple-precision floating-point value
 4937 `a' by the corresponding value `b'.  The operation is performed according to
 4938 the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 4939 -------------------------------------------------------------------------------
 4940 */
float128 float128_div( float128 a, float128 b )
 4942 {
flag aSign, bSign, zSign;
int32 aExp, bExp, zExp;
bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
float128 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
bSig1 = extractFloat128Frac1( b );
bSig0 = extractFloat128Frac0( b );
bExp = extractFloat128Exp( b );
bSign = extractFloat128Sign( b );
zSign = aSign ^ bSign;
 4958     if ( aExp == 0x7FFF ) {
 4959         if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, b );
 4960         if ( bExp == 0x7FFF ) {
 4961             if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
 4962             goto invalid;
 4963         }
 4964         return packFloat128( zSign, 0x7FFF, 0, 0 );
 4965     }
 4966     if ( bExp == 0x7FFF ) {
 4967         if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
 4968         return packFloat128( zSign, 0, 0, 0 );
 4969     }
 4970     if ( bExp == 0 ) {
 4971         if ( ( bSig0 | bSig1 ) == 0 ) {
 4972             if ( ( aExp | aSig0 | aSig1 ) == 0 ) {
invalid:
float_raise( float_flag_invalid );
z.low = float128_default_nan_low;
z.high = float128_default_nan_high;
 4977                 return z;
 4978             }
float_raise( float_flag_divbyzero );
 4980             return packFloat128( zSign, 0x7FFF, 0, 0 );
 4981         }
normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
 4983     }
 4984     if ( aExp == 0 ) {
 4985         if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( zSign, 0, 0, 0 );
normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
 4987     }
zExp = aExp - bExp + 0x3FFD;
shortShift128Left(
aSig0 | LIT64( 0x0001000000000000 ), aSig1, 15, &aSig0, &aSig1 );
shortShift128Left(
bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
 4993     if ( le128( bSig0, bSig1, aSig0, aSig1 ) ) {
shift128Right( aSig0, aSig1, 1, &aSig0, &aSig1 );
 4995         ++zExp;
 4996     }
zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
 5000     while ( (sbits64) rem0 < 0 ) {
 5001         --zSig0;
add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
 5003     }
zSig1 = estimateDiv128To64( rem1, rem2, bSig0 );
 5005     if ( ( zSig1 & 0x3FFF ) <= 4 ) {
mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
 5008         while ( (sbits64) rem1 < 0 ) {
 5009             --zSig1;
add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
 5011         }
zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
 5013     }
shift128ExtraRightJamming( zSig0, zSig1, 0, 15, &zSig0, &zSig1, &zSig2 );
 5015     return roundAndPackFloat128( zSign, zExp, zSig0, zSig1, zSig2 );
 5017 }
 5019 /*
 5020 -------------------------------------------------------------------------------
 5021 Returns the remainder of the quadruple-precision floating-point value `a'
 5022 with respect to the corresponding value `b'.  The operation is performed
 5023 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 5024 -------------------------------------------------------------------------------
 5025 */
float128 float128_rem( float128 a, float128 b )
 5027 {
flag aSign, bSign, zSign;
int32 aExp, bExp, expDiff;
bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
bits64 allZero, alternateASig0, alternateASig1, sigMean1;
sbits64 sigMean0;
float128 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
bSig1 = extractFloat128Frac1( b );
bSig0 = extractFloat128Frac0( b );
bExp = extractFloat128Exp( b );
bSign = extractFloat128Sign( b );
 5043     if ( aExp == 0x7FFF ) {
 5044         if (    ( aSig0 | aSig1 )
 5045              || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
 5046             return propagateFloat128NaN( a, b );
 5047         }
 5048         goto invalid;
 5049     }
 5050     if ( bExp == 0x7FFF ) {
 5051         if ( bSig0 | bSig1 ) return propagateFloat128NaN( a, b );
 5052         return a;
 5053     }
 5054     if ( bExp == 0 ) {
 5055         if ( ( bSig0 | bSig1 ) == 0 ) {
invalid:
float_raise( float_flag_invalid );
z.low = float128_default_nan_low;
z.high = float128_default_nan_high;
 5060             return z;
 5061         }
normalizeFloat128Subnormal( bSig0, bSig1, &bExp, &bSig0, &bSig1 );
 5063     }
 5064     if ( aExp == 0 ) {
 5065         if ( ( aSig0 | aSig1 ) == 0 ) return a;
normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
 5067     }
expDiff = aExp - bExp;
 5069     if ( expDiff < -1 ) return a;
shortShift128Left(
aSig0 | LIT64( 0x0001000000000000 ),
aSig1,
 5073         15 - ( expDiff < 0 ),
 5074         &aSig0,
 5075         &aSig1
 5076     );
shortShift128Left(
bSig0 | LIT64( 0x0001000000000000 ), bSig1, 15, &bSig0, &bSig1 );
q = le128( bSig0, bSig1, aSig0, aSig1 );
 5080     if ( q ) sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
expDiff -= 64;
 5082     while ( 0 < expDiff ) {
q = estimateDiv128To64( aSig0, aSig1, bSig0 );
q = ( 4 < q ) ? q - 4 : 0;
mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
shortShift192Left( term0, term1, term2, 61, &term1, &term2, &allZero );
shortShift128Left( aSig0, aSig1, 61, &aSig0, &allZero );
sub128( aSig0, 0, term1, term2, &aSig0, &aSig1 );
expDiff -= 61;
 5090     }
 5091     if ( -64 < expDiff ) {
q = estimateDiv128To64( aSig0, aSig1, bSig0 );
q = ( 4 < q ) ? q - 4 : 0;
q >>= - expDiff;
shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
expDiff += 52;
 5097         if ( expDiff < 0 ) {
shift128Right( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
 5099         }
 5100         else {
shortShift128Left( aSig0, aSig1, expDiff, &aSig0, &aSig1 );
 5102         }
mul128By64To192( bSig0, bSig1, q, &term0, &term1, &term2 );
sub128( aSig0, aSig1, term1, term2, &aSig0, &aSig1 );
 5105     }
 5106     else {
shift128Right( aSig0, aSig1, 12, &aSig0, &aSig1 );
shift128Right( bSig0, bSig1, 12, &bSig0, &bSig1 );
 5109     }
 5110     do {
alternateASig0 = aSig0;
alternateASig1 = aSig1;
 5113         ++q;
sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
 5115     } while ( 0 <= (sbits64) aSig0 );
add128(
aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 );
 5118     if (    ( sigMean0 < 0 )
 5119          || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
aSig0 = alternateASig0;
aSig1 = alternateASig1;
 5122     }
zSign = ( (sbits64) aSig0 < 0 );
 5124     if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
 5125     return
normalizeRoundAndPackFloat128( aSign ^ zSign, bExp - 4, aSig0, aSig1 );
 5128 }
 5130 /*
 5131 -------------------------------------------------------------------------------
 5132 Returns the square root of the quadruple-precision floating-point value `a'.
 5133 The operation is performed according to the IEC/IEEE Standard for Binary
 5134 Floating-Point Arithmetic.
 5135 -------------------------------------------------------------------------------
 5136 */
float128 float128_sqrt( float128 a )
 5138 {
flag aSign;
int32 aExp, zExp;
bits64 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
float128 z;
aSig1 = extractFloat128Frac1( a );
aSig0 = extractFloat128Frac0( a );
aExp = extractFloat128Exp( a );
aSign = extractFloat128Sign( a );
 5149     if ( aExp == 0x7FFF ) {
 5150         if ( aSig0 | aSig1 ) return propagateFloat128NaN( a, a );
 5151         if ( ! aSign ) return a;
 5152         goto invalid;
 5153     }
 5154     if ( aSign ) {
 5155         if ( ( aExp | aSig0 | aSig1 ) == 0 ) return a;
invalid:
float_raise( float_flag_invalid );
z.low = float128_default_nan_low;
z.high = float128_default_nan_high;
 5160         return z;
 5161     }
 5162     if ( aExp == 0 ) {
 5163         if ( ( aSig0 | aSig1 ) == 0 ) return packFloat128( 0, 0, 0, 0 );
normalizeFloat128Subnormal( aSig0, aSig1, &aExp, &aSig0, &aSig1 );
 5165     }
zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFE;
aSig0 |= LIT64( 0x0001000000000000 );
zSig0 = estimateSqrt32( aExp, aSig0>>17 );
shortShift128Left( aSig0, aSig1, 13 - ( aExp & 1 ), &aSig0, &aSig1 );
zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0<<32 ) + ( zSig0<<30 );
doubleZSig0 = zSig0<<1;
mul64To128( zSig0, zSig0, &term0, &term1 );
sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
 5174     while ( (sbits64) rem0 < 0 ) {
 5175         --zSig0;
doubleZSig0 -= 2;
add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
 5178     }
zSig1 = estimateDiv128To64( rem1, 0, doubleZSig0 );
 5180     if ( ( zSig1 & 0x1FFF ) <= 5 ) {
 5181         if ( zSig1 == 0 ) zSig1 = 1;
mul64To128( doubleZSig0, zSig1, &term1, &term2 );
sub128( rem1, 0, term1, term2, &rem1, &rem2 );
mul64To128( zSig1, zSig1, &term2, &term3 );
sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
 5186         while ( (sbits64) rem1 < 0 ) {
 5187             --zSig1;
shortShift128Left( 0, zSig1, 1, &term2, &term3 );
term3 |= 1;
term2 |= doubleZSig0;
add192( rem1, rem2, rem3, 0, term2, term3, &rem1, &rem2, &rem3 );
 5192         }
zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
 5194     }
shift128ExtraRightJamming( zSig0, zSig1, 0, 14, &zSig0, &zSig1, &zSig2 );
 5196     return roundAndPackFloat128( 0, zExp, zSig0, zSig1, zSig2 );
 5198 }
 5200 /*
 5201 -------------------------------------------------------------------------------
 5202 Returns 1 if the quadruple-precision floating-point value `a' is equal to
 5203 the corresponding value `b', and 0 otherwise.  The comparison is performed
 5204 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 5205 -------------------------------------------------------------------------------
 5206 */
flag float128_eq( float128 a, float128 b )
 5208 {
 5210     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
 5211               && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
 5212          || (    ( extractFloat128Exp( b ) == 0x7FFF )
 5213               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
 5214        ) {
 5215         if (    float128_is_signaling_nan( a )
 5216              || float128_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 5218         }
 5219         return 0;
 5220     }
 5221     return
 5222            ( a.low == b.low )
 5223         && (    ( a.high == b.high )
 5224              || (    ( a.low == 0 )
 5225                   && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
 5226            );
 5228 }
 5230 /*
 5231 -------------------------------------------------------------------------------
 5232 Returns 1 if the quadruple-precision floating-point value `a' is less than
 5233 or equal to the corresponding value `b', and 0 otherwise.  The comparison
 5234 is performed according to the IEC/IEEE Standard for Binary Floating-Point
 5235 Arithmetic.
 5236 -------------------------------------------------------------------------------
 5237 */
flag float128_le( float128 a, float128 b )
 5239 {
flag aSign, bSign;
 5242     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
 5243               && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
 5244          || (    ( extractFloat128Exp( b ) == 0x7FFF )
 5245               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
 5246        ) {
float_raise( float_flag_invalid );
 5248         return 0;
 5249     }
aSign = extractFloat128Sign( a );
bSign = extractFloat128Sign( b );
 5252     if ( aSign != bSign ) {
 5253         return
aSign
 5255             || (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 5256                  == 0 );
 5257     }
 5258     return
aSign ? le128( b.high, b.low, a.high, a.low )
 5260         : le128( a.high, a.low, b.high, b.low );
 5262 }
 5264 /*
 5265 -------------------------------------------------------------------------------
 5266 Returns 1 if the quadruple-precision floating-point value `a' is less than
 5267 the corresponding value `b', and 0 otherwise.  The comparison is performed
 5268 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 5269 -------------------------------------------------------------------------------
 5270 */
flag float128_lt( float128 a, float128 b )
 5272 {
flag aSign, bSign;
 5275     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
 5276               && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
 5277          || (    ( extractFloat128Exp( b ) == 0x7FFF )
 5278               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
 5279        ) {
float_raise( float_flag_invalid );
 5281         return 0;
 5282     }
aSign = extractFloat128Sign( a );
bSign = extractFloat128Sign( b );
 5285     if ( aSign != bSign ) {
 5286         return
aSign
 5288             && (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 5289                  != 0 );
 5290     }
 5291     return
aSign ? lt128( b.high, b.low, a.high, a.low )
 5293         : lt128( a.high, a.low, b.high, b.low );
 5295 }
 5297 /*
 5298 -------------------------------------------------------------------------------
 5299 Returns 1 if the quadruple-precision floating-point value `a' is equal to
 5300 the corresponding value `b', and 0 otherwise.  The invalid exception is
 5301 raised if either operand is a NaN.  Otherwise, the comparison is performed
 5302 according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 5303 -------------------------------------------------------------------------------
 5304 */
flag float128_eq_signaling( float128 a, float128 b )
 5306 {
 5308     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
 5309               && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
 5310          || (    ( extractFloat128Exp( b ) == 0x7FFF )
 5311               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
 5312        ) {
float_raise( float_flag_invalid );
 5314         return 0;
 5315     }
 5316     return
 5317            ( a.low == b.low )
 5318         && (    ( a.high == b.high )
 5319              || (    ( a.low == 0 )
 5320                   && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
 5321            );
 5323 }
 5325 /*
 5326 -------------------------------------------------------------------------------
 5327 Returns 1 if the quadruple-precision floating-point value `a' is less than
 5328 or equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
 5329 cause an exception.  Otherwise, the comparison is performed according to the
 5330 IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 5331 -------------------------------------------------------------------------------
 5332 */
flag float128_le_quiet( float128 a, float128 b )
 5334 {
flag aSign, bSign;
 5337     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
 5338               && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
 5339          || (    ( extractFloat128Exp( b ) == 0x7FFF )
 5340               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
 5341        ) {
 5342         if (    float128_is_signaling_nan( a )
 5343              || float128_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 5345         }
 5346         return 0;
 5347     }
aSign = extractFloat128Sign( a );
bSign = extractFloat128Sign( b );
 5350     if ( aSign != bSign ) {
 5351         return
aSign
 5353             || (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 5354                  == 0 );
 5355     }
 5356     return
aSign ? le128( b.high, b.low, a.high, a.low )
 5358         : le128( a.high, a.low, b.high, b.low );
 5360 }
 5362 /*
 5363 -------------------------------------------------------------------------------
 5364 Returns 1 if the quadruple-precision floating-point value `a' is less than
 5365 the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
 5366 exception.  Otherwise, the comparison is performed according to the IEC/IEEE
 5367 Standard for Binary Floating-Point Arithmetic.
 5368 -------------------------------------------------------------------------------
 5369 */
flag float128_lt_quiet( float128 a, float128 b )
 5371 {
flag aSign, bSign;
 5374     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
 5375               && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
 5376          || (    ( extractFloat128Exp( b ) == 0x7FFF )
 5377               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
 5378        ) {
 5379         if (    float128_is_signaling_nan( a )
 5380              || float128_is_signaling_nan( b ) ) {
float_raise( float_flag_invalid );
 5382         }
 5383         return 0;
 5384     }
aSign = extractFloat128Sign( a );
bSign = extractFloat128Sign( b );
 5387     if ( aSign != bSign ) {
 5388         return
aSign
 5390             && (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
 5391                  != 0 );
 5392     }
 5393     return
aSign ? lt128( b.high, b.low, a.high, a.low )
 5395         : lt128( a.high, a.low, b.high, b.low );
 5397 }
 5399 #endif
 5402 #if defined(SOFTFLOAT_FOR_GCC) && defined(SOFTFLOAT_NEED_FIXUNS)
 5404 /*
 5405  * These two routines are not part of the original softfloat distribution.
 5406  *
 5407  * They are based on the corresponding conversions to integer but return
 5408  * unsigned numbers instead since these functions are required by GCC.
 5409  *
 5410  * Added by Mark Brinicombe <mark@netbsd.org>   27/09/97
 5411  *
 5412  * float64 version overhauled for SoftFloat 2a [bjh21 2000-07-15]
 5413  */
 5415 /*
 5416 -------------------------------------------------------------------------------
 5417 Returns the result of converting the double-precision floating-point value
 5418 `a' to the 32-bit unsigned integer format.  The conversion is
 5419 performed according to the IEC/IEEE Standard for Binary Floating-point
 5420 Arithmetic, except that the conversion is always rounded toward zero.  If
 5421 `a' is a NaN, the largest positive integer is returned.  If the conversion
 5422 overflows, the largest integer positive is returned.
 5423 -------------------------------------------------------------------------------
 5424 */
uint32 float64_to_uint32_round_to_zero( float64 a )
 5426 {
flag aSign;
int16 aExp, shiftCount;
bits64 aSig, savedASig;
uint32 z;
aSig = extractFloat64Frac( a );
aExp = extractFloat64Exp( a );
aSign = extractFloat64Sign( a );
 5436     if (aSign) {
float_raise( float_flag_invalid );
 5438         return(0);
 5439     }
 5441     if ( 0x41E < aExp ) {
float_raise( float_flag_invalid );
 5443         return 0xffffffff;
 5444     }
 5445     else if ( aExp < 0x3FF ) {
 5446         if ( aExp || aSig ) float_set_inexact();
 5447         return 0;
 5448     }
aSig |= LIT64( 0x0010000000000000 );
shiftCount = 0x433 - aExp;
savedASig = aSig;
aSig >>= shiftCount;
z = aSig;
 5454     if ( ( aSig<<shiftCount ) != savedASig ) {
float_set_inexact();
 5456     }
 5457     return z;
 5459 }
 5461 /*
 5462 -------------------------------------------------------------------------------
 5463 Returns the result of converting the single-precision floating-point value
 5464 `a' to the 32-bit unsigned integer format.  The conversion is
 5465 performed according to the IEC/IEEE Standard for Binary Floating-point
 5466 Arithmetic, except that the conversion is always rounded toward zero.  If
 5467 `a' is a NaN, the largest positive integer is returned.  If the conversion
 5468 overflows, the largest positive integer is returned.
 5469 -------------------------------------------------------------------------------
 5470 */
uint32 float32_to_uint32_round_to_zero( float32 a )
 5472 {
flag aSign;
int16 aExp, shiftCount;
bits32 aSig;
uint32 z;
aSig = extractFloat32Frac( a );
aExp = extractFloat32Exp( a );
aSign = extractFloat32Sign( a );
shiftCount = aExp - 0x9E;
 5483     if (aSign) {
float_raise( float_flag_invalid );
 5485         return(0);
 5486     }
 5487     if ( 0 < shiftCount ) {
float_raise( float_flag_invalid );
 5489         return 0xFFFFFFFF;
 5490     }
 5491     else if ( aExp <= 0x7E ) {
 5492         if ( aExp | aSig ) float_set_inexact();
 5493         return 0;
 5494     }
aSig = ( aSig | 0x800000 )<<8;
z = aSig>>( - shiftCount );
 5497     if ( aSig<<( shiftCount & 31 ) ) {
float_set_inexact();
 5499     }
 5500     return z;
 5502 }
 5504 #endif
 5506 #endif /* !NO_IEEE */