]> git.proxmox.com Git - mirror_qemu.git/blobdiff - fpu/softfloat.c
softfloat: Replace uint16 type with uint_fast16_t
[mirror_qemu.git] / fpu / softfloat.c
index 0b8279798cb1eeee70c70501e9b0d418cac858a0..50fce01c3a3b4973ff1fa54a9f351083ee71dd37 100644 (file)
@@ -1,3 +1,8 @@
+/*
+ * QEMU float support
+ *
+ * Derived from SoftFloat.
+ */
 
 /*============================================================================
 
@@ -30,8 +35,11 @@ these four paragraphs for those parts of this code that are retained.
 
 =============================================================================*/
 
-/* FIXME: Flush-To-Zero only effects results.  Denormal inputs should also
-   be flushed to zero.  */
+/* softfloat (and in particular the code in softfloat-specialize.h) is
+ * target-dependent and needs the TARGET_* macros.
+ */
+#include "config.h"
+
 #include "softfloat.h"
 
 /*----------------------------------------------------------------------------
@@ -61,12 +69,37 @@ void set_float_exception_flags(int val STATUS_PARAM)
     STATUS(float_exception_flags) = val;
 }
 
-#ifdef FLOATX80
 void set_floatx80_rounding_precision(int val STATUS_PARAM)
 {
     STATUS(floatx80_rounding_precision) = val;
 }
-#endif
+
+/*----------------------------------------------------------------------------
+| Returns the fraction bits of the half-precision floating-point value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE uint32_t extractFloat16Frac(float16 a)
+{
+    return float16_val(a) & 0x3ff;
+}
+
+/*----------------------------------------------------------------------------
+| Returns the exponent bits of the half-precision floating-point value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE int16 extractFloat16Exp(float16 a)
+{
+    return (float16_val(a) >> 10) & 0x1f;
+}
+
+/*----------------------------------------------------------------------------
+| Returns the sign bit of the single-precision floating-point value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE flag extractFloat16Sign(float16 a)
+{
+    return float16_val(a)>>15;
+}
 
 /*----------------------------------------------------------------------------
 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
@@ -79,12 +112,12 @@ void set_floatx80_rounding_precision(int val STATUS_PARAM)
 | positive or negative integer is returned.
 *----------------------------------------------------------------------------*/
 
-static int32 roundAndPackInt32( flag zSign, bits64 absZ STATUS_PARAM)
+static int32 roundAndPackInt32( flag zSign, uint64_t absZ STATUS_PARAM)
 {
     int8 roundingMode;
     flag roundNearestEven;
     int8 roundIncrement, roundBits;
-    int32 z;
+    int32_t z;
 
     roundingMode = STATUS(float_rounding_mode);
     roundNearestEven = ( roundingMode == float_round_nearest_even );
@@ -110,7 +143,7 @@ static int32 roundAndPackInt32( flag zSign, bits64 absZ STATUS_PARAM)
     if ( zSign ) z = - z;
     if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
         float_raise( float_flag_invalid STATUS_VAR);
-        return zSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
+        return zSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
     }
     if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
     return z;
@@ -129,15 +162,15 @@ static int32 roundAndPackInt32( flag zSign, bits64 absZ STATUS_PARAM)
 | returned.
 *----------------------------------------------------------------------------*/
 
-static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PARAM)
+static int64 roundAndPackInt64( flag zSign, uint64_t absZ0, uint64_t absZ1 STATUS_PARAM)
 {
     int8 roundingMode;
     flag roundNearestEven, increment;
-    int64 z;
+    int64_t z;
 
     roundingMode = STATUS(float_rounding_mode);
     roundNearestEven = ( roundingMode == float_round_nearest_even );
-    increment = ( (sbits64) absZ1 < 0 );
+    increment = ( (int64_t) absZ1 < 0 );
     if ( ! roundNearestEven ) {
         if ( roundingMode == float_round_to_zero ) {
             increment = 0;
@@ -154,7 +187,7 @@ static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PA
     if ( increment ) {
         ++absZ0;
         if ( absZ0 == 0 ) goto overflow;
-        absZ0 &= ~ ( ( (bits64) ( absZ1<<1 ) == 0 ) & roundNearestEven );
+        absZ0 &= ~ ( ( (uint64_t) ( absZ1<<1 ) == 0 ) & roundNearestEven );
     }
     z = absZ0;
     if ( zSign ) z = - z;
@@ -162,7 +195,7 @@ static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PA
  overflow:
         float_raise( float_flag_invalid STATUS_VAR);
         return
-              zSign ? (sbits64) LIT64( 0x8000000000000000 )
+              zSign ? (int64_t) LIT64( 0x8000000000000000 )
             : LIT64( 0x7FFFFFFFFFFFFFFF );
     }
     if ( absZ1 ) STATUS(float_exception_flags) |= float_flag_inexact;
@@ -174,7 +207,7 @@ static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PA
 | Returns the fraction bits of the single-precision floating-point value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE bits32 extractFloat32Frac( float32 a )
+INLINE uint32_t extractFloat32Frac( float32 a )
 {
 
     return float32_val(a) & 0x007FFFFF;
@@ -203,6 +236,21 @@ INLINE flag extractFloat32Sign( float32 a )
 
 }
 
+/*----------------------------------------------------------------------------
+| If `a' is denormal and we are in flush-to-zero mode then set the
+| input-denormal exception and return zero. Otherwise just return the value.
+*----------------------------------------------------------------------------*/
+static float32 float32_squash_input_denormal(float32 a STATUS_PARAM)
+{
+    if (STATUS(flush_inputs_to_zero)) {
+        if (extractFloat32Exp(a) == 0 && extractFloat32Frac(a) != 0) {
+            float_raise(float_flag_input_denormal STATUS_VAR);
+            return make_float32(float32_val(a) & 0x80000000);
+        }
+    }
+    return a;
+}
+
 /*----------------------------------------------------------------------------
 | Normalizes the subnormal single-precision floating-point value represented
 | by the denormalized significand `aSig'.  The normalized exponent and
@@ -211,7 +259,7 @@ INLINE flag extractFloat32Sign( float32 a )
 *----------------------------------------------------------------------------*/
 
 static void
- normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
+ normalizeFloat32Subnormal( uint32_t aSig, int16 *zExpPtr, uint32_t *zSigPtr )
 {
     int8 shiftCount;
 
@@ -232,11 +280,11 @@ static void
 | significand.
 *----------------------------------------------------------------------------*/
 
-INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
+INLINE float32 packFloat32( flag zSign, int16 zExp, uint32_t zSig )
 {
 
     return make_float32(
-          ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig);
+          ( ( (uint32_t) zSign )<<31 ) + ( ( (uint32_t) zExp )<<23 ) + zSig);
 
 }
 
@@ -262,7 +310,7 @@ INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
 | Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
+static float32 roundAndPackFloat32( flag zSign, int16 zExp, uint32_t zSig STATUS_PARAM)
 {
     int8 roundingMode;
     flag roundNearestEven;
@@ -287,16 +335,19 @@ static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_P
         }
     }
     roundBits = zSig & 0x7F;
-    if ( 0xFD <= (bits16) zExp ) {
+    if ( 0xFD <= (uint16_t) zExp ) {
         if (    ( 0xFD < zExp )
              || (    ( zExp == 0xFD )
-                  && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
+                  && ( (int32_t) ( zSig + roundIncrement ) < 0 ) )
            ) {
             float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
             return packFloat32( zSign, 0xFF, - ( roundIncrement == 0 ));
         }
         if ( zExp < 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloat32( zSign, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                float_raise(float_flag_output_denormal STATUS_VAR);
+                return packFloat32(zSign, 0, 0);
+            }
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < -1 )
@@ -325,7 +376,7 @@ static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_P
 *----------------------------------------------------------------------------*/
 
 static float32
- normalizeRoundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_PARAM)
+ normalizeRoundAndPackFloat32( flag zSign, int16 zExp, uint32_t zSig STATUS_PARAM)
 {
     int8 shiftCount;
 
@@ -338,7 +389,7 @@ static float32
 | Returns the fraction bits of the double-precision floating-point value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE bits64 extractFloat64Frac( float64 a )
+INLINE uint64_t extractFloat64Frac( float64 a )
 {
 
     return float64_val(a) & LIT64( 0x000FFFFFFFFFFFFF );
@@ -367,6 +418,21 @@ INLINE flag extractFloat64Sign( float64 a )
 
 }
 
+/*----------------------------------------------------------------------------
+| If `a' is denormal and we are in flush-to-zero mode then set the
+| input-denormal exception and return zero. Otherwise just return the value.
+*----------------------------------------------------------------------------*/
+static float64 float64_squash_input_denormal(float64 a STATUS_PARAM)
+{
+    if (STATUS(flush_inputs_to_zero)) {
+        if (extractFloat64Exp(a) == 0 && extractFloat64Frac(a) != 0) {
+            float_raise(float_flag_input_denormal STATUS_VAR);
+            return make_float64(float64_val(a) & (1ULL << 63));
+        }
+    }
+    return a;
+}
+
 /*----------------------------------------------------------------------------
 | Normalizes the subnormal double-precision floating-point value represented
 | by the denormalized significand `aSig'.  The normalized exponent and
@@ -375,7 +441,7 @@ INLINE flag extractFloat64Sign( float64 a )
 *----------------------------------------------------------------------------*/
 
 static void
- normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
+ normalizeFloat64Subnormal( uint64_t aSig, int16 *zExpPtr, uint64_t *zSigPtr )
 {
     int8 shiftCount;
 
@@ -396,11 +462,11 @@ static void
 | significand.
 *----------------------------------------------------------------------------*/
 
-INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
+INLINE float64 packFloat64( flag zSign, int16 zExp, uint64_t zSig )
 {
 
     return make_float64(
-        ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig);
+        ( ( (uint64_t) zSign )<<63 ) + ( ( (uint64_t) zExp )<<52 ) + zSig);
 
 }
 
@@ -426,7 +492,7 @@ INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
 | Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
+static float64 roundAndPackFloat64( flag zSign, int16 zExp, uint64_t zSig STATUS_PARAM)
 {
     int8 roundingMode;
     flag roundNearestEven;
@@ -451,16 +517,19 @@ static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_P
         }
     }
     roundBits = zSig & 0x3FF;
-    if ( 0x7FD <= (bits16) zExp ) {
+    if ( 0x7FD <= (uint16_t) zExp ) {
         if (    ( 0x7FD < zExp )
              || (    ( zExp == 0x7FD )
-                  && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
+                  && ( (int64_t) ( zSig + roundIncrement ) < 0 ) )
            ) {
             float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
             return packFloat64( zSign, 0x7FF, - ( roundIncrement == 0 ));
         }
         if ( zExp < 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloat64( zSign, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                float_raise(float_flag_output_denormal STATUS_VAR);
+                return packFloat64(zSign, 0, 0);
+            }
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < -1 )
@@ -489,7 +558,7 @@ static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_P
 *----------------------------------------------------------------------------*/
 
 static float64
- normalizeRoundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_PARAM)
+ normalizeRoundAndPackFloat64( flag zSign, int16 zExp, uint64_t zSig STATUS_PARAM)
 {
     int8 shiftCount;
 
@@ -498,14 +567,12 @@ static float64
 
 }
 
-#ifdef FLOATX80
-
 /*----------------------------------------------------------------------------
 | Returns the fraction bits of the extended double-precision floating-point
 | value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE bits64 extractFloatx80Frac( floatx80 a )
+INLINE uint64_t extractFloatx80Frac( floatx80 a )
 {
 
     return a.low;
@@ -544,7 +611,7 @@ INLINE flag extractFloatx80Sign( floatx80 a )
 *----------------------------------------------------------------------------*/
 
 static void
- normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
+ normalizeFloatx80Subnormal( uint64_t aSig, int32 *zExpPtr, uint64_t *zSigPtr )
 {
     int8 shiftCount;
 
@@ -559,12 +626,12 @@ static void
 | extended double-precision floating-point value, returning the result.
 *----------------------------------------------------------------------------*/
 
-INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
+INLINE floatx80 packFloatx80( flag zSign, int32 zExp, uint64_t zSig )
 {
     floatx80 z;
 
     z.low = zSig;
-    z.high = ( ( (bits16) zSign )<<15 ) + zExp;
+    z.high = ( ( (uint16_t) zSign )<<15 ) + zExp;
     return z;
 
 }
@@ -595,7 +662,7 @@ INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
 
 static floatx80
  roundAndPackFloatx80(
-     int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
+     int8 roundingPrecision, flag zSign, int32 zExp, uint64_t zSig0, uint64_t zSig1
  STATUS_PARAM)
 {
     int8 roundingMode;
@@ -632,14 +699,17 @@ static floatx80
         }
     }
     roundBits = zSig0 & roundMask;
-    if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
+    if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
         if (    ( 0x7FFE < zExp )
              || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
            ) {
             goto overflow;
         }
         if ( zExp <= 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloatx80( zSign, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                float_raise(float_flag_output_denormal STATUS_VAR);
+                return packFloatx80(zSign, 0, 0);
+            }
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < 0 )
@@ -650,7 +720,7 @@ static floatx80
             if ( isTiny && roundBits ) float_raise( float_flag_underflow STATUS_VAR);
             if ( roundBits ) STATUS(float_exception_flags) |= float_flag_inexact;
             zSig0 += roundIncrement;
-            if ( (sbits64) zSig0 < 0 ) zExp = 1;
+            if ( (int64_t) zSig0 < 0 ) zExp = 1;
             roundIncrement = roundMask + 1;
             if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
                 roundMask |= roundIncrement;
@@ -673,7 +743,7 @@ static floatx80
     if ( zSig0 == 0 ) zExp = 0;
     return packFloatx80( zSign, zExp, zSig0 );
  precision80:
-    increment = ( (sbits64) zSig1 < 0 );
+    increment = ( (int64_t) zSig1 < 0 );
     if ( ! roundNearestEven ) {
         if ( roundingMode == float_round_to_zero ) {
             increment = 0;
@@ -687,7 +757,7 @@ static floatx80
             }
         }
     }
-    if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
+    if ( 0x7FFD <= (uint32_t) ( zExp - 1 ) ) {
         if (    ( 0x7FFE < zExp )
              || (    ( zExp == 0x7FFE )
                   && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
@@ -716,7 +786,7 @@ static floatx80
             if ( isTiny && zSig1 ) float_raise( float_flag_underflow STATUS_VAR);
             if ( zSig1 ) STATUS(float_exception_flags) |= float_flag_inexact;
             if ( roundNearestEven ) {
-                increment = ( (sbits64) zSig1 < 0 );
+                increment = ( (int64_t) zSig1 < 0 );
             }
             else {
                 if ( zSign ) {
@@ -729,8 +799,8 @@ static floatx80
             if ( increment ) {
                 ++zSig0;
                 zSig0 &=
-                    ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
-                if ( (sbits64) zSig0 < 0 ) zExp = 1;
+                    ~ ( ( (uint64_t) ( zSig1<<1 ) == 0 ) & roundNearestEven );
+                if ( (int64_t) zSig0 < 0 ) zExp = 1;
             }
             return packFloatx80( zSign, zExp, zSig0 );
         }
@@ -743,7 +813,7 @@ static floatx80
             zSig0 = LIT64( 0x8000000000000000 );
         }
         else {
-            zSig0 &= ~ ( ( (bits64) ( zSig1<<1 ) == 0 ) & roundNearestEven );
+            zSig0 &= ~ ( ( (uint64_t) ( zSig1<<1 ) == 0 ) & roundNearestEven );
         }
     }
     else {
@@ -764,7 +834,7 @@ static floatx80
 
 static floatx80
  normalizeRoundAndPackFloatx80(
-     int8 roundingPrecision, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
+     int8 roundingPrecision, flag zSign, int32 zExp, uint64_t zSig0, uint64_t zSig1
  STATUS_PARAM)
 {
     int8 shiftCount;
@@ -782,16 +852,12 @@ static floatx80
 
 }
 
-#endif
-
-#ifdef FLOAT128
-
 /*----------------------------------------------------------------------------
 | Returns the least-significant 64 fraction bits of the quadruple-precision
 | floating-point value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE bits64 extractFloat128Frac1( float128 a )
+INLINE uint64_t extractFloat128Frac1( float128 a )
 {
 
     return a.low;
@@ -803,7 +869,7 @@ INLINE bits64 extractFloat128Frac1( float128 a )
 | floating-point value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE bits64 extractFloat128Frac0( float128 a )
+INLINE uint64_t extractFloat128Frac0( float128 a )
 {
 
     return a.high & LIT64( 0x0000FFFFFFFFFFFF );
@@ -845,11 +911,11 @@ INLINE flag extractFloat128Sign( float128 a )
 
 static void
  normalizeFloat128Subnormal(
-     bits64 aSig0,
-     bits64 aSig1,
+     uint64_t aSig0,
+     uint64_t aSig1,
      int32 *zExpPtr,
-     bits64 *zSig0Ptr,
-     bits64 *zSig1Ptr
+     uint64_t *zSig0Ptr,
+     uint64_t *zSig1Ptr
  )
 {
     int8 shiftCount;
@@ -888,12 +954,12 @@ static void
 *----------------------------------------------------------------------------*/
 
 INLINE float128
- packFloat128( flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 )
+ packFloat128( flag zSign, int32 zExp, uint64_t zSig0, uint64_t zSig1 )
 {
     float128 z;
 
     z.low = zSig1;
-    z.high = ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<48 ) + zSig0;
+    z.high = ( ( (uint64_t) zSign )<<63 ) + ( ( (uint64_t) zExp )<<48 ) + zSig0;
     return z;
 
 }
@@ -921,14 +987,14 @@ INLINE float128
 
 static float128
  roundAndPackFloat128(
-     flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1, bits64 zSig2 STATUS_PARAM)
+     flag zSign, int32 zExp, uint64_t zSig0, uint64_t zSig1, uint64_t zSig2 STATUS_PARAM)
 {
     int8 roundingMode;
     flag roundNearestEven, increment, isTiny;
 
     roundingMode = STATUS(float_rounding_mode);
     roundNearestEven = ( roundingMode == float_round_nearest_even );
-    increment = ( (sbits64) zSig2 < 0 );
+    increment = ( (int64_t) zSig2 < 0 );
     if ( ! roundNearestEven ) {
         if ( roundingMode == float_round_to_zero ) {
             increment = 0;
@@ -942,7 +1008,7 @@ static float128
             }
         }
     }
-    if ( 0x7FFD <= (bits32) zExp ) {
+    if ( 0x7FFD <= (uint32_t) zExp ) {
         if (    ( 0x7FFD < zExp )
              || (    ( zExp == 0x7FFD )
                   && eq128(
@@ -970,7 +1036,10 @@ static float128
             return packFloat128( zSign, 0x7FFF, 0, 0 );
         }
         if ( zExp < 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloat128( zSign, 0, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                float_raise(float_flag_output_denormal STATUS_VAR);
+                return packFloat128(zSign, 0, 0, 0);
+            }
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < -1 )
@@ -986,7 +1055,7 @@ static float128
             zExp = 0;
             if ( isTiny && zSig2 ) float_raise( float_flag_underflow STATUS_VAR);
             if ( roundNearestEven ) {
-                increment = ( (sbits64) zSig2 < 0 );
+                increment = ( (int64_t) zSig2 < 0 );
             }
             else {
                 if ( zSign ) {
@@ -1022,10 +1091,10 @@ static float128
 
 static float128
  normalizeRoundAndPackFloat128(
-     flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1 STATUS_PARAM)
+     flag zSign, int32 zExp, uint64_t zSig0, uint64_t zSig1 STATUS_PARAM)
 {
     int8 shiftCount;
-    bits64 zSig2;
+    uint64_t zSig2;
 
     if ( zSig0 == 0 ) {
         zSig0 = zSig1;
@@ -1046,8 +1115,6 @@ static float128
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 32-bit two's complement integer `a'
 | to the single-precision floating-point format.  The conversion is performed
@@ -1059,7 +1126,7 @@ float32 int32_to_float32( int32 a STATUS_PARAM )
     flag zSign;
 
     if ( a == 0 ) return float32_zero;
-    if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
+    if ( a == (int32_t) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
     zSign = ( a < 0 );
     return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a STATUS_VAR );
 
@@ -1076,7 +1143,7 @@ float64 int32_to_float64( int32 a STATUS_PARAM )
     flag zSign;
     uint32 absA;
     int8 shiftCount;
-    bits64 zSig;
+    uint64_t zSig;
 
     if ( a == 0 ) return float64_zero;
     zSign = ( a < 0 );
@@ -1087,8 +1154,6 @@ float64 int32_to_float64( int32 a STATUS_PARAM )
 
 }
 
-#ifdef FLOATX80
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 32-bit two's complement integer `a'
 | to the extended double-precision floating-point format.  The conversion
@@ -1101,7 +1166,7 @@ floatx80 int32_to_floatx80( int32 a STATUS_PARAM )
     flag zSign;
     uint32 absA;
     int8 shiftCount;
-    bits64 zSig;
+    uint64_t zSig;
 
     if ( a == 0 ) return packFloatx80( 0, 0, 0 );
     zSign = ( a < 0 );
@@ -1112,10 +1177,6 @@ floatx80 int32_to_floatx80( int32 a STATUS_PARAM )
 
 }
 
-#endif
-
-#ifdef FLOAT128
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 32-bit two's complement integer `a' to
 | the quadruple-precision floating-point format.  The conversion is performed
@@ -1127,7 +1188,7 @@ float128 int32_to_float128( int32 a STATUS_PARAM )
     flag zSign;
     uint32 absA;
     int8 shiftCount;
-    bits64 zSig0;
+    uint64_t zSig0;
 
     if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
     zSign = ( a < 0 );
@@ -1138,8 +1199,6 @@ float128 int32_to_float128( int32 a STATUS_PARAM )
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 64-bit two's complement integer `a'
 | to the single-precision floating-point format.  The conversion is performed
@@ -1204,7 +1263,7 @@ float64 int64_to_float64( int64 a STATUS_PARAM )
     flag zSign;
 
     if ( a == 0 ) return float64_zero;
-    if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
+    if ( a == (int64_t) LIT64( 0x8000000000000000 ) ) {
         return packFloat64( 1, 0x43E, 0 );
     }
     zSign = ( a < 0 );
@@ -1219,8 +1278,6 @@ float64 uint64_to_float64( uint64 a STATUS_PARAM )
 
 }
 
-#ifdef FLOATX80
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 64-bit two's complement integer `a'
 | to the extended double-precision floating-point format.  The conversion
@@ -1242,10 +1299,6 @@ floatx80 int64_to_floatx80( int64 a STATUS_PARAM )
 
 }
 
-#endif
-
-#ifdef FLOAT128
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 64-bit two's complement integer `a' to
 | the quadruple-precision floating-point format.  The conversion is performed
@@ -1258,7 +1311,7 @@ float128 int64_to_float128( int64 a STATUS_PARAM )
     uint64 absA;
     int8 shiftCount;
     int32 zExp;
-    bits64 zSig0, zSig1;
+    uint64_t zSig0, zSig1;
 
     if ( a == 0 ) return packFloat128( 0, 0, 0, 0 );
     zSign = ( a < 0 );
@@ -1279,8 +1332,6 @@ float128 int64_to_float128( int64 a STATUS_PARAM )
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the single-precision floating-point value
 | `a' to the 32-bit two's complement integer format.  The conversion is
@@ -1295,9 +1346,10 @@ int32 float32_to_int32( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits32 aSig;
-    bits64 aSig64;
+    uint32_t aSig;
+    uint64_t aSig64;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
@@ -1325,8 +1377,9 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits32 aSig;
-    int32 z;
+    uint32_t aSig;
+    int32_t z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1337,7 +1390,7 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
             float_raise( float_flag_invalid STATUS_VAR);
             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
         }
-        return (sbits32) 0x80000000;
+        return (int32_t) 0x80000000;
     }
     else if ( aExp <= 0x7E ) {
         if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
@@ -1345,7 +1398,7 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
     }
     aSig = ( aSig | 0x00800000 )<<8;
     z = aSig>>( - shiftCount );
-    if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
+    if ( (uint32_t) ( aSig<<( shiftCount & 31 ) ) ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
     }
     if ( aSign ) z = - z;
@@ -1353,6 +1406,55 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the result of converting the single-precision floating-point value
+| `a' to the 16-bit two's complement integer format.  The conversion is
+| performed according to the IEC/IEEE Standard for Binary Floating-Point
+| Arithmetic, except that the conversion is always rounded toward zero.
+| If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
+| the conversion overflows, the largest integer with the same sign as `a' is
+| returned.
+*----------------------------------------------------------------------------*/
+
+int16 float32_to_int16_round_to_zero( float32 a STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    uint32_t aSig;
+    int32 z;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    shiftCount = aExp - 0x8E;
+    if ( 0 <= shiftCount ) {
+        if ( float32_val(a) != 0xC7000000 ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+            if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
+                return 0x7FFF;
+            }
+        }
+        return (int32_t) 0xffff8000;
+    }
+    else if ( aExp <= 0x7E ) {
+        if ( aExp | aSig ) {
+            STATUS(float_exception_flags) |= float_flag_inexact;
+        }
+        return 0;
+    }
+    shiftCount -= 0x10;
+    aSig = ( aSig | 0x00800000 )<<8;
+    z = aSig>>( - shiftCount );
+    if ( (uint32_t) ( aSig<<( shiftCount & 31 ) ) ) {
+        STATUS(float_exception_flags) |= float_flag_inexact;
+    }
+    if ( aSign ) {
+        z = - z;
+    }
+    return z;
+
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of converting the single-precision floating-point value
 | `a' to the 64-bit two's complement integer format.  The conversion is
@@ -1367,8 +1469,9 @@ int64 float32_to_int64( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits32 aSig;
-    bits64 aSig64, aSigExtra;
+    uint32_t aSig;
+    uint64_t aSig64, aSigExtra;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1379,7 +1482,7 @@ int64 float32_to_int64( float32 a STATUS_PARAM )
         if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
             return LIT64( 0x7FFFFFFFFFFFFFFF );
         }
-        return (sbits64) LIT64( 0x8000000000000000 );
+        return (int64_t) LIT64( 0x8000000000000000 );
     }
     if ( aExp ) aSig |= 0x00800000;
     aSig64 = aSig;
@@ -1403,9 +1506,10 @@ int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits32 aSig;
-    bits64 aSig64;
+    uint32_t aSig;
+    uint64_t aSig64;
     int64 z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1418,7 +1522,7 @@ int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
             }
         }
-        return (sbits64) LIT64( 0x8000000000000000 );
+        return (int64_t) LIT64( 0x8000000000000000 );
     }
     else if ( aExp <= 0x7E ) {
         if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
@@ -1427,7 +1531,7 @@ int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
     aSig64 = aSig | 0x00800000;
     aSig64 <<= 40;
     z = aSig64>>( - shiftCount );
-    if ( (bits64) ( aSig64<<( shiftCount & 63 ) ) ) {
+    if ( (uint64_t) ( aSig64<<( shiftCount & 63 ) ) ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
     }
     if ( aSign ) z = - z;
@@ -1446,13 +1550,14 @@ float64 float32_to_float64( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits32 aSig;
+    uint32_t aSig;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     if ( aExp == 0xFF ) {
-        if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ));
+        if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloat64( aSign, 0x7FF, 0 );
     }
     if ( aExp == 0 ) {
@@ -1460,12 +1565,10 @@ float64 float32_to_float64( float32 a STATUS_PARAM )
         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
         --aExp;
     }
-    return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
+    return packFloat64( aSign, aExp + 0x380, ( (uint64_t) aSig )<<29 );
 
 }
 
-#ifdef FLOATX80
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the single-precision floating-point value
 | `a' to the extended double-precision floating-point format.  The conversion
@@ -1477,13 +1580,14 @@ floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits32 aSig;
+    uint32_t aSig;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     if ( aExp == 0xFF ) {
-        if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( aExp == 0 ) {
@@ -1491,14 +1595,10 @@ floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
     }
     aSig |= 0x00800000;
-    return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
+    return packFloatx80( aSign, aExp + 0x3F80, ( (uint64_t) aSig )<<40 );
 
 }
 
-#endif
-
-#ifdef FLOAT128
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the single-precision floating-point value
 | `a' to the double-precision floating-point format.  The conversion is
@@ -1510,13 +1610,14 @@ float128 float32_to_float128( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits32 aSig;
+    uint32_t aSig;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     if ( aExp == 0xFF ) {
-        if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloat128( aSign, 0x7FFF, 0, 0 );
     }
     if ( aExp == 0 ) {
@@ -1524,12 +1625,10 @@ float128 float32_to_float128( float32 a STATUS_PARAM )
         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
         --aExp;
     }
-    return packFloat128( aSign, aExp + 0x3F80, ( (bits64) aSig )<<25, 0 );
+    return packFloat128( aSign, aExp + 0x3F80, ( (uint64_t) aSig )<<25, 0 );
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Rounds the single-precision floating-point value `a' to an integer, and
 | returns the result as a single-precision floating-point value.  The
@@ -1541,9 +1640,10 @@ float32 float32_round_to_int( float32 a STATUS_PARAM)
 {
     flag aSign;
     int16 aExp;
-    bits32 lastBitMask, roundBitsMask;
+    uint32_t lastBitMask, roundBitsMask;
     int8 roundingMode;
-    bits32 z;
+    uint32_t z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aExp = extractFloat32Exp( a );
     if ( 0x96 <= aExp ) {
@@ -1553,7 +1653,7 @@ float32 float32_round_to_int( float32 a STATUS_PARAM)
         return a;
     }
     if ( aExp <= 0x7E ) {
-        if ( (bits32) ( float32_val(a)<<1 ) == 0 ) return a;
+        if ( (uint32_t) ( float32_val(a)<<1 ) == 0 ) return a;
         STATUS(float_exception_flags) |= float_flag_inexact;
         aSign = extractFloat32Sign( a );
         switch ( STATUS(float_rounding_mode) ) {
@@ -1600,7 +1700,7 @@ float32 float32_round_to_int( float32 a STATUS_PARAM)
 static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
 {
     int16 aExp, bExp, zExp;
-    bits32 aSig, bSig, zSig;
+    uint32_t aSig, bSig, zSig;
     int16 expDiff;
 
     aSig = extractFloat32Frac( a );
@@ -1644,7 +1744,12 @@ static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
             return a;
         }
         if ( aExp == 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloat32( zSign, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                if (aSig | bSig) {
+                    float_raise(float_flag_output_denormal STATUS_VAR);
+                }
+                return packFloat32(zSign, 0, 0);
+            }
             return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
         }
         zSig = 0x40000000 + aSig + bSig;
@@ -1654,7 +1759,7 @@ static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
     aSig |= 0x20000000;
     zSig = ( aSig + bSig )<<1;
     --zExp;
-    if ( (sbits32) zSig < 0 ) {
+    if ( (int32_t) zSig < 0 ) {
         zSig = aSig + bSig;
         ++zExp;
     }
@@ -1674,7 +1779,7 @@ static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
 static float32 subFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
 {
     int16 aExp, bExp, zExp;
-    bits32 aSig, bSig, zSig;
+    uint32_t aSig, bSig, zSig;
     int16 expDiff;
 
     aSig = extractFloat32Frac( a );
@@ -1747,6 +1852,8 @@ static float32 subFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
 float32 float32_add( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     aSign = extractFloat32Sign( a );
     bSign = extractFloat32Sign( b );
@@ -1768,6 +1875,8 @@ float32 float32_add( float32 a, float32 b STATUS_PARAM )
 float32 float32_sub( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     aSign = extractFloat32Sign( a );
     bSign = extractFloat32Sign( b );
@@ -1790,9 +1899,12 @@ float32 float32_mul( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int16 aExp, bExp, zExp;
-    bits32 aSig, bSig;
-    bits64 zSig64;
-    bits32 zSig;
+    uint32_t aSig, bSig;
+    uint64_t zSig64;
+    uint32_t zSig;
+
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1830,9 +1942,9 @@ float32 float32_mul( float32 a, float32 b STATUS_PARAM )
     zExp = aExp + bExp - 0x7F;
     aSig = ( aSig | 0x00800000 )<<7;
     bSig = ( bSig | 0x00800000 )<<8;
-    shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
+    shift64RightJamming( ( (uint64_t) aSig ) * bSig, 32, &zSig64 );
     zSig = zSig64;
-    if ( 0 <= (sbits32) ( zSig<<1 ) ) {
+    if ( 0 <= (int32_t) ( zSig<<1 ) ) {
         zSig <<= 1;
         --zExp;
     }
@@ -1850,7 +1962,9 @@ float32 float32_div( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int16 aExp, bExp, zExp;
-    bits32 aSig, bSig, zSig;
+    uint32_t aSig, bSig, zSig;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1894,9 +2008,9 @@ float32 float32_div( float32 a, float32 b STATUS_PARAM )
         aSig >>= 1;
         ++zExp;
     }
-    zSig = ( ( (bits64) aSig )<<32 ) / bSig;
+    zSig = ( ( (uint64_t) aSig )<<32 ) / bSig;
     if ( ( zSig & 0x3F ) == 0 ) {
-        zSig |= ( (bits64) bSig * zSig != ( (bits64) aSig )<<32 );
+        zSig |= ( (uint64_t) bSig * zSig != ( (uint64_t) aSig )<<32 );
     }
     return roundAndPackFloat32( zSign, zExp, zSig STATUS_VAR );
 
@@ -1912,11 +2026,13 @@ float32 float32_rem( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, zSign;
     int16 aExp, bExp, expDiff;
-    bits32 aSig, bSig;
-    bits32 q;
-    bits64 aSig64, bSig64, q64;
-    bits32 alternateASig;
-    sbits32 sigMean;
+    uint32_t aSig, bSig;
+    uint32_t q;
+    uint64_t aSig64, bSig64, q64;
+    uint32_t alternateASig;
+    int32_t sigMean;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1958,7 +2074,7 @@ float32 float32_rem( float32 a, float32 b STATUS_PARAM )
         q = ( bSig <= aSig );
         if ( q ) aSig -= bSig;
         if ( 0 < expDiff ) {
-            q = ( ( (bits64) aSig )<<32 ) / bSig;
+            q = ( ( (uint64_t) aSig )<<32 ) / bSig;
             q >>= 32 - expDiff;
             bSig >>= 2;
             aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
@@ -1970,8 +2086,8 @@ float32 float32_rem( float32 a, float32 b STATUS_PARAM )
     }
     else {
         if ( bSig <= aSig ) aSig -= bSig;
-        aSig64 = ( (bits64) aSig )<<40;
-        bSig64 = ( (bits64) bSig )<<40;
+        aSig64 = ( (uint64_t) aSig )<<40;
+        bSig64 = ( (uint64_t) bSig )<<40;
         expDiff -= 64;
         while ( 0 < expDiff ) {
             q64 = estimateDiv128To64( aSig64, 0, bSig64 );
@@ -1990,17 +2106,224 @@ float32 float32_rem( float32 a, float32 b STATUS_PARAM )
         alternateASig = aSig;
         ++q;
         aSig -= bSig;
-    } while ( 0 <= (sbits32) aSig );
+    } while ( 0 <= (int32_t) aSig );
     sigMean = aSig + alternateASig;
     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
         aSig = alternateASig;
     }
-    zSign = ( (sbits32) aSig < 0 );
+    zSign = ( (int32_t) aSig < 0 );
     if ( zSign ) aSig = - aSig;
     return normalizeRoundAndPackFloat32( aSign ^ zSign, bExp, aSig STATUS_VAR );
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the result of multiplying the single-precision floating-point values
+| `a' and `b' then adding 'c', with no intermediate rounding step after the
+| multiplication.  The operation is performed according to the IEC/IEEE
+| Standard for Binary Floating-Point Arithmetic 754-2008.
+| The flags argument allows the caller to select negation of the
+| addend, the intermediate product, or the final result. (The difference
+| between this and having the caller do a separate negation is that negating
+| externally will flip the sign bit on NaNs.)
+*----------------------------------------------------------------------------*/
+
+float32 float32_muladd(float32 a, float32 b, float32 c, int flags STATUS_PARAM)
+{
+    flag aSign, bSign, cSign, zSign;
+    int16 aExp, bExp, cExp, pExp, zExp, expDiff;
+    uint32_t aSig, bSig, cSig;
+    flag pInf, pZero, pSign;
+    uint64_t pSig64, cSig64, zSig64;
+    uint32_t pSig;
+    int shiftcount;
+    flag signflip, infzero;
+
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
+    c = float32_squash_input_denormal(c STATUS_VAR);
+    aSig = extractFloat32Frac(a);
+    aExp = extractFloat32Exp(a);
+    aSign = extractFloat32Sign(a);
+    bSig = extractFloat32Frac(b);
+    bExp = extractFloat32Exp(b);
+    bSign = extractFloat32Sign(b);
+    cSig = extractFloat32Frac(c);
+    cExp = extractFloat32Exp(c);
+    cSign = extractFloat32Sign(c);
+
+    infzero = ((aExp == 0 && aSig == 0 && bExp == 0xff && bSig == 0) ||
+               (aExp == 0xff && aSig == 0 && bExp == 0 && bSig == 0));
+
+    /* It is implementation-defined whether the cases of (0,inf,qnan)
+     * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
+     * they return if they do), so we have to hand this information
+     * off to the target-specific pick-a-NaN routine.
+     */
+    if (((aExp == 0xff) && aSig) ||
+        ((bExp == 0xff) && bSig) ||
+        ((cExp == 0xff) && cSig)) {
+        return propagateFloat32MulAddNaN(a, b, c, infzero STATUS_VAR);
+    }
+
+    if (infzero) {
+        float_raise(float_flag_invalid STATUS_VAR);
+        return float32_default_nan;
+    }
+
+    if (flags & float_muladd_negate_c) {
+        cSign ^= 1;
+    }
+
+    signflip = (flags & float_muladd_negate_result) ? 1 : 0;
+
+    /* Work out the sign and type of the product */
+    pSign = aSign ^ bSign;
+    if (flags & float_muladd_negate_product) {
+        pSign ^= 1;
+    }
+    pInf = (aExp == 0xff) || (bExp == 0xff);
+    pZero = ((aExp | aSig) == 0) || ((bExp | bSig) == 0);
+
+    if (cExp == 0xff) {
+        if (pInf && (pSign ^ cSign)) {
+            /* addition of opposite-signed infinities => InvalidOperation */
+            float_raise(float_flag_invalid STATUS_VAR);
+            return float32_default_nan;
+        }
+        /* Otherwise generate an infinity of the same sign */
+        return packFloat32(cSign ^ signflip, 0xff, 0);
+    }
+
+    if (pInf) {
+        return packFloat32(pSign ^ signflip, 0xff, 0);
+    }
+
+    if (pZero) {
+        if (cExp == 0) {
+            if (cSig == 0) {
+                /* Adding two exact zeroes */
+                if (pSign == cSign) {
+                    zSign = pSign;
+                } else if (STATUS(float_rounding_mode) == float_round_down) {
+                    zSign = 1;
+                } else {
+                    zSign = 0;
+                }
+                return packFloat32(zSign ^ signflip, 0, 0);
+            }
+            /* Exact zero plus a denorm */
+            if (STATUS(flush_to_zero)) {
+                float_raise(float_flag_output_denormal STATUS_VAR);
+                return packFloat32(cSign ^ signflip, 0, 0);
+            }
+        }
+        /* Zero plus something non-zero : just return the something */
+        return make_float32(float32_val(c) ^ (signflip << 31));
+    }
+
+    if (aExp == 0) {
+        normalizeFloat32Subnormal(aSig, &aExp, &aSig);
+    }
+    if (bExp == 0) {
+        normalizeFloat32Subnormal(bSig, &bExp, &bSig);
+    }
+
+    /* Calculate the actual result a * b + c */
+
+    /* Multiply first; this is easy. */
+    /* NB: we subtract 0x7e where float32_mul() subtracts 0x7f
+     * because we want the true exponent, not the "one-less-than"
+     * flavour that roundAndPackFloat32() takes.
+     */
+    pExp = aExp + bExp - 0x7e;
+    aSig = (aSig | 0x00800000) << 7;
+    bSig = (bSig | 0x00800000) << 8;
+    pSig64 = (uint64_t)aSig * bSig;
+    if ((int64_t)(pSig64 << 1) >= 0) {
+        pSig64 <<= 1;
+        pExp--;
+    }
+
+    zSign = pSign ^ signflip;
+
+    /* Now pSig64 is the significand of the multiply, with the explicit bit in
+     * position 62.
+     */
+    if (cExp == 0) {
+        if (!cSig) {
+            /* Throw out the special case of c being an exact zero now */
+            shift64RightJamming(pSig64, 32, &pSig64);
+            pSig = pSig64;
+            return roundAndPackFloat32(zSign, pExp - 1,
+                                       pSig STATUS_VAR);
+        }
+        normalizeFloat32Subnormal(cSig, &cExp, &cSig);
+    }
+
+    cSig64 = (uint64_t)cSig << (62 - 23);
+    cSig64 |= LIT64(0x4000000000000000);
+    expDiff = pExp - cExp;
+
+    if (pSign == cSign) {
+        /* Addition */
+        if (expDiff > 0) {
+            /* scale c to match p */
+            shift64RightJamming(cSig64, expDiff, &cSig64);
+            zExp = pExp;
+        } else if (expDiff < 0) {
+            /* scale p to match c */
+            shift64RightJamming(pSig64, -expDiff, &pSig64);
+            zExp = cExp;
+        } else {
+            /* no scaling needed */
+            zExp = cExp;
+        }
+        /* Add significands and make sure explicit bit ends up in posn 62 */
+        zSig64 = pSig64 + cSig64;
+        if ((int64_t)zSig64 < 0) {
+            shift64RightJamming(zSig64, 1, &zSig64);
+        } else {
+            zExp--;
+        }
+    } else {
+        /* Subtraction */
+        if (expDiff > 0) {
+            shift64RightJamming(cSig64, expDiff, &cSig64);
+            zSig64 = pSig64 - cSig64;
+            zExp = pExp;
+        } else if (expDiff < 0) {
+            shift64RightJamming(pSig64, -expDiff, &pSig64);
+            zSig64 = cSig64 - pSig64;
+            zExp = cExp;
+            zSign ^= 1;
+        } else {
+            zExp = pExp;
+            if (cSig64 < pSig64) {
+                zSig64 = pSig64 - cSig64;
+            } else if (pSig64 < cSig64) {
+                zSig64 = cSig64 - pSig64;
+                zSign ^= 1;
+            } else {
+                /* Exact zero */
+                zSign = signflip;
+                if (STATUS(float_rounding_mode) == float_round_down) {
+                    zSign ^= 1;
+                }
+                return packFloat32(zSign, 0, 0);
+            }
+        }
+        --zExp;
+        /* Normalize to put the explicit bit back into bit 62. */
+        shiftcount = countLeadingZeros64(zSig64) - 1;
+        zSig64 <<= shiftcount;
+        zExp -= shiftcount;
+    }
+    shift64RightJamming(zSig64, 32, &zSig64);
+    return roundAndPackFloat32(zSign, zExp, zSig64 STATUS_VAR);
+}
+
+
 /*----------------------------------------------------------------------------
 | Returns the square root of the single-precision floating-point value `a'.
 | The operation is performed according to the IEC/IEEE Standard for Binary
@@ -2011,8 +2334,9 @@ float32 float32_sqrt( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, zExp;
-    bits32 aSig, zSig;
-    bits64 rem, term;
+    uint32_t aSig, zSig;
+    uint64_t rem, term;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -2041,11 +2365,11 @@ float32 float32_sqrt( float32 a STATUS_PARAM )
             goto roundAndPack;
         }
         aSig >>= aExp & 1;
-        term = ( (bits64) zSig ) * zSig;
-        rem = ( ( (bits64) aSig )<<32 ) - term;
-        while ( (sbits64) rem < 0 ) {
+        term = ( (uint64_t) zSig ) * zSig;
+        rem = ( ( (uint64_t) aSig )<<32 ) - term;
+        while ( (int64_t) rem < 0 ) {
             --zSig;
-            rem += ( ( (bits64) zSig )<<1 ) | 1;
+            rem += ( ( (uint64_t) zSig )<<1 ) | 1;
         }
         zSig |= ( rem != 0 );
     }
@@ -2075,30 +2399,31 @@ float32 float32_sqrt( float32 a STATUS_PARAM )
 
 static const float64 float32_exp2_coefficients[15] =
 {
-    make_float64( 0x3ff0000000000000ll ), /*  1 */
-    make_float64( 0x3fe0000000000000ll ), /*  2 */
-    make_float64( 0x3fc5555555555555ll ), /*  3 */
-    make_float64( 0x3fa5555555555555ll ), /*  4 */
-    make_float64( 0x3f81111111111111ll ), /*  5 */
-    make_float64( 0x3f56c16c16c16c17ll ), /*  6 */
-    make_float64( 0x3f2a01a01a01a01all ), /*  7 */
-    make_float64( 0x3efa01a01a01a01all ), /*  8 */
-    make_float64( 0x3ec71de3a556c734ll ), /*  9 */
-    make_float64( 0x3e927e4fb7789f5cll ), /* 10 */
-    make_float64( 0x3e5ae64567f544e4ll ), /* 11 */
-    make_float64( 0x3e21eed8eff8d898ll ), /* 12 */
-    make_float64( 0x3de6124613a86d09ll ), /* 13 */
-    make_float64( 0x3da93974a8c07c9dll ), /* 14 */
-    make_float64( 0x3d6ae7f3e733b81fll ), /* 15 */
+    const_float64( 0x3ff0000000000000ll ), /*  1 */
+    const_float64( 0x3fe0000000000000ll ), /*  2 */
+    const_float64( 0x3fc5555555555555ll ), /*  3 */
+    const_float64( 0x3fa5555555555555ll ), /*  4 */
+    const_float64( 0x3f81111111111111ll ), /*  5 */
+    const_float64( 0x3f56c16c16c16c17ll ), /*  6 */
+    const_float64( 0x3f2a01a01a01a01all ), /*  7 */
+    const_float64( 0x3efa01a01a01a01all ), /*  8 */
+    const_float64( 0x3ec71de3a556c734ll ), /*  9 */
+    const_float64( 0x3e927e4fb7789f5cll ), /* 10 */
+    const_float64( 0x3e5ae64567f544e4ll ), /* 11 */
+    const_float64( 0x3e21eed8eff8d898ll ), /* 12 */
+    const_float64( 0x3de6124613a86d09ll ), /* 13 */
+    const_float64( 0x3da93974a8c07c9dll ), /* 14 */
+    const_float64( 0x3d6ae7f3e733b81fll ), /* 15 */
 };
 
 float32 float32_exp2( float32 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits32 aSig;
+    uint32_t aSig;
     float64 r, x, xn;
     int i;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -2143,8 +2468,9 @@ float32 float32_log2( float32 a STATUS_PARAM )
 {
     flag aSign, zSign;
     int16 aExp;
-    bits32 aSig, zSig, i;
+    uint32_t aSig, zSig, i;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
@@ -2168,7 +2494,7 @@ float32 float32_log2( float32 a STATUS_PARAM )
     zSig = aExp << 23;
 
     for (i = 1 << 22; i > 0; i >>= 1) {
-        aSig = ( (bits64)aSig * aSig ) >> 23;
+        aSig = ( (uint64_t)aSig * aSig ) >> 23;
         if ( aSig & 0x01000000 ) {
             aSig >>= 1;
             zSig |= i;
@@ -2183,37 +2509,41 @@ float32 float32_log2( float32 a STATUS_PARAM )
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the single-precision floating-point value `a' is equal to
-| the corresponding value `b', and 0 otherwise.  The comparison is performed
+| the corresponding value `b', and 0 otherwise.  The invalid exception is
+| raised if either operand is a NaN.  Otherwise, the comparison is performed
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float32_eq( float32 a, float32 b STATUS_PARAM )
 {
+    uint32_t av, bv;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
        ) {
-        if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
-            float_raise( float_flag_invalid STATUS_VAR);
-        }
+        float_raise( float_flag_invalid STATUS_VAR);
         return 0;
     }
-    return ( float32_val(a) == float32_val(b) ) ||
-            ( (bits32) ( ( float32_val(a) | float32_val(b) )<<1 ) == 0 );
-
+    av = float32_val(a);
+    bv = float32_val(b);
+    return ( av == bv ) || ( (uint32_t) ( ( av | bv )<<1 ) == 0 );
 }
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the single-precision floating-point value `a' is less than
-| or equal to the corresponding value `b', and 0 otherwise.  The comparison
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
+| or equal to the corresponding value `b', and 0 otherwise.  The invalid
+| exception is raised if either operand is a NaN.  The comparison is performed
+| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float32_le( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits32 av, bv;
+    uint32_t av, bv;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
@@ -2225,21 +2555,24 @@ int float32_le( float32 a, float32 b STATUS_PARAM )
     bSign = extractFloat32Sign( b );
     av = float32_val(a);
     bv = float32_val(b);
-    if ( aSign != bSign ) return aSign || ( (bits32) ( ( av | bv )<<1 ) == 0 );
+    if ( aSign != bSign ) return aSign || ( (uint32_t) ( ( av | bv )<<1 ) == 0 );
     return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the single-precision floating-point value `a' is less than
-| the corresponding value `b', and 0 otherwise.  The comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| the corresponding value `b', and 0 otherwise.  The invalid exception is
+| raised if either operand is a NaN.  The comparison is performed according
+| to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float32_lt( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits32 av, bv;
+    uint32_t av, bv;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
@@ -2251,32 +2584,54 @@ int float32_lt( float32 a, float32 b STATUS_PARAM )
     bSign = extractFloat32Sign( b );
     av = float32_val(a);
     bv = float32_val(b);
-    if ( aSign != bSign ) return aSign && ( (bits32) ( ( av | bv )<<1 ) != 0 );
+    if ( aSign != bSign ) return aSign && ( (uint32_t) ( ( av | bv )<<1 ) != 0 );
     return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
 /*----------------------------------------------------------------------------
-| Returns 1 if the single-precision floating-point value `a' is equal to
-| the corresponding value `b', and 0 otherwise.  The invalid exception is
-| raised if either operand is a NaN.  Otherwise, the comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| Returns 1 if the single-precision floating-point values `a' and `b' cannot
+| be compared, and 0 otherwise.  The invalid exception is raised if either
+| operand is a NaN.  The comparison is performed according to the IEC/IEEE
+| Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-int float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
+int float32_unordered( float32 a, float32 b STATUS_PARAM )
 {
-    bits32 av, bv;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
        ) {
         float_raise( float_flag_invalid STATUS_VAR);
-        return 0;
+        return 1;
     }
-    av = float32_val(a);
-    bv = float32_val(b);
-    return ( av == bv ) || ( (bits32) ( ( av | bv )<<1 ) == 0 );
+    return 0;
+}
 
+/*----------------------------------------------------------------------------
+| Returns 1 if the single-precision floating-point value `a' is equal to
+| the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
+| exception.  The comparison is performed according to the IEC/IEEE Standard
+| for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+int float32_eq_quiet( float32 a, float32 b STATUS_PARAM )
+{
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return 0;
+    }
+    return ( float32_val(a) == float32_val(b) ) ||
+            ( (uint32_t) ( ( float32_val(a) | float32_val(b) )<<1 ) == 0 );
 }
 
 /*----------------------------------------------------------------------------
@@ -2289,7 +2644,9 @@ int float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
 int float32_le_quiet( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits32 av, bv;
+    uint32_t av, bv;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
@@ -2303,7 +2660,7 @@ int float32_le_quiet( float32 a, float32 b STATUS_PARAM )
     bSign = extractFloat32Sign( b );
     av = float32_val(a);
     bv = float32_val(b);
-    if ( aSign != bSign ) return aSign || ( (bits32) ( ( av | bv )<<1 ) == 0 );
+    if ( aSign != bSign ) return aSign || ( (uint32_t) ( ( av | bv )<<1 ) == 0 );
     return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
@@ -2318,7 +2675,9 @@ int float32_le_quiet( float32 a, float32 b STATUS_PARAM )
 int float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits32 av, bv;
+    uint32_t av, bv;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
          || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
@@ -2332,11 +2691,34 @@ int float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
     bSign = extractFloat32Sign( b );
     av = float32_val(a);
     bv = float32_val(b);
-    if ( aSign != bSign ) return aSign && ( (bits32) ( ( av | bv )<<1 ) != 0 );
+    if ( aSign != bSign ) return aSign && ( (uint32_t) ( ( av | bv )<<1 ) != 0 );
     return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns 1 if the single-precision floating-point values `a' and `b' cannot
+| be compared, and 0 otherwise.  Quiet NaNs do not cause an exception.  The
+| comparison is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+int float32_unordered_quiet( float32 a, float32 b STATUS_PARAM )
+{
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return 1;
+    }
+    return 0;
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of converting the double-precision floating-point value
 | `a' to the 32-bit two's complement integer format.  The conversion is
@@ -2351,7 +2733,8 @@ int32 float64_to_int32( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits64 aSig;
+    uint64_t aSig;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2378,8 +2761,9 @@ int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits64 aSig, savedASig;
-    int32 z;
+    uint64_t aSig, savedASig;
+    int32_t z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2401,7 +2785,7 @@ int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
     if ( ( z < 0 ) ^ aSign ) {
  invalid:
         float_raise( float_flag_invalid STATUS_VAR);
-        return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
+        return aSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
     }
     if ( ( aSig<<shiftCount ) != savedASig ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
@@ -2410,6 +2794,57 @@ int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the result of converting the double-precision floating-point value
+| `a' to the 16-bit two's complement integer format.  The conversion is
+| performed according to the IEC/IEEE Standard for Binary Floating-Point
+| Arithmetic, except that the conversion is always rounded toward zero.
+| If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
+| the conversion overflows, the largest integer with the same sign as `a' is
+| returned.
+*----------------------------------------------------------------------------*/
+
+int16 float64_to_int16_round_to_zero( float64 a STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    uint64_t aSig, savedASig;
+    int32 z;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    if ( 0x40E < aExp ) {
+        if ( ( aExp == 0x7FF ) && aSig ) {
+            aSign = 0;
+        }
+        goto invalid;
+    }
+    else if ( aExp < 0x3FF ) {
+        if ( aExp || aSig ) {
+            STATUS(float_exception_flags) |= float_flag_inexact;
+        }
+        return 0;
+    }
+    aSig |= LIT64( 0x0010000000000000 );
+    shiftCount = 0x433 - aExp;
+    savedASig = aSig;
+    aSig >>= shiftCount;
+    z = aSig;
+    if ( aSign ) {
+        z = - z;
+    }
+    if ( ( (int16_t)z < 0 ) ^ aSign ) {
+ invalid:
+        float_raise( float_flag_invalid STATUS_VAR);
+        return aSign ? (int32_t) 0xffff8000 : 0x7FFF;
+    }
+    if ( ( aSig<<shiftCount ) != savedASig ) {
+        STATUS(float_exception_flags) |= float_flag_inexact;
+    }
+    return z;
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of converting the double-precision floating-point value
 | `a' to the 64-bit two's complement integer format.  The conversion is
@@ -2424,7 +2859,8 @@ int64 float64_to_int64( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits64 aSig, aSigExtra;
+    uint64_t aSig, aSigExtra;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2440,7 +2876,7 @@ int64 float64_to_int64( float64 a STATUS_PARAM )
                ) {
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
             }
-            return (sbits64) LIT64( 0x8000000000000000 );
+            return (int64_t) LIT64( 0x8000000000000000 );
         }
         aSigExtra = 0;
         aSig <<= - shiftCount;
@@ -2466,8 +2902,9 @@ int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, shiftCount;
-    bits64 aSig;
+    uint64_t aSig;
     int64 z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2485,7 +2922,7 @@ int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
                     return LIT64( 0x7FFFFFFFFFFFFFFF );
                 }
             }
-            return (sbits64) LIT64( 0x8000000000000000 );
+            return (int64_t) LIT64( 0x8000000000000000 );
         }
         z = aSig<<shiftCount;
     }
@@ -2495,7 +2932,7 @@ int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
             return 0;
         }
         z = aSig>>( - shiftCount );
-        if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
+        if ( (uint64_t) ( aSig<<( shiftCount & 63 ) ) ) {
             STATUS(float_exception_flags) |= float_flag_inexact;
         }
     }
@@ -2515,14 +2952,15 @@ float32 float64_to_float32( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits64 aSig;
-    bits32 zSig;
+    uint64_t aSig;
+    uint32_t zSig;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
     aSign = extractFloat64Sign( a );
     if ( aExp == 0x7FF ) {
-        if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloat32( aSign, 0xFF, 0 );
     }
     shift64RightJamming( aSig, 22, &aSig );
@@ -2546,29 +2984,28 @@ float32 float64_to_float32( float64 a STATUS_PARAM )
 | than the desired result exponent whenever `zSig' is a complete, normalized
 | significand.
 *----------------------------------------------------------------------------*/
-static bits16 packFloat16(flag zSign, int16 zExp, bits16 zSig)
+static float16 packFloat16(flag zSign, int16 zExp, uint16_t zSig)
 {
-    return (((bits32)zSign) << 15) + (((bits32)zExp) << 10) + zSig;
+    return make_float16(
+        (((uint32_t)zSign) << 15) + (((uint32_t)zExp) << 10) + zSig);
 }
 
 /* Half precision floats come in two formats: standard IEEE and "ARM" format.
    The latter gains extra exponent range by omitting the NaN/Inf encodings.  */
-  
-float32 float16_to_float32( bits16 a, flag ieee STATUS_PARAM )
+
+float32 float16_to_float32(float16 a, flag ieee STATUS_PARAM)
 {
     flag aSign;
     int16 aExp;
-    bits32 aSig;
+    uint32_t aSig;
 
-    aSign = a >> 15;
-    aExp = (a >> 10) & 0x1f;
-    aSig = a & 0x3ff;
+    aSign = extractFloat16Sign(a);
+    aExp = extractFloat16Exp(a);
+    aSig = extractFloat16Frac(a);
 
     if (aExp == 0x1f && ieee) {
         if (aSig) {
-            /* Make sure correct exceptions are raised.  */
-            float32ToCommonNaN(a STATUS_VAR);
-            aSig |= 0x200;
+            return commonNaNToFloat32(float16ToCommonNaN(a STATUS_VAR) STATUS_VAR);
         }
         return packFloat32(aSign, 0xff, aSig << 13);
     }
@@ -2586,38 +3023,45 @@ float32 float16_to_float32( bits16 a, flag ieee STATUS_PARAM )
     return packFloat32( aSign, aExp + 0x70, aSig << 13);
 }
 
-bits16 float32_to_float16( float32 a, flag ieee STATUS_PARAM)
+float16 float32_to_float16(float32 a, flag ieee STATUS_PARAM)
 {
     flag aSign;
     int16 aExp;
-    bits32 aSig;
-    bits32 mask;
-    bits32 increment;
+    uint32_t aSig;
+    uint32_t mask;
+    uint32_t increment;
     int8 roundingMode;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     if ( aExp == 0xFF ) {
         if (aSig) {
-            /* Make sure correct exceptions are raised.  */
-            float32ToCommonNaN(a STATUS_VAR);
-            aSig |= 0x00400000;
+            /* Input is a NaN */
+            float16 r = commonNaNToFloat16( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
+            if (!ieee) {
+                return packFloat16(aSign, 0, 0);
+            }
+            return r;
+        }
+        /* Infinity */
+        if (!ieee) {
+            float_raise(float_flag_invalid STATUS_VAR);
+            return packFloat16(aSign, 0x1f, 0x3ff);
         }
-        return packFloat16(aSign, 0x1f, aSig >> 13);
+        return packFloat16(aSign, 0x1f, 0);
     }
-    if (aExp == 0 && aSign == 0) {
+    if (aExp == 0 && aSig == 0) {
         return packFloat16(aSign, 0, 0);
     }
     /* Decimal point between bits 22 and 23.  */
     aSig |= 0x00800000;
     aExp -= 0x7f;
     if (aExp < -14) {
-        mask = 0x007fffff;
-        if (aExp < -24) {
-            aExp = -25;
-        } else {
-            mask >>= 24 + aExp;
+        mask = 0x00ffffff;
+        if (aExp >= -24) {
+            mask >>= 25 + aExp;
         }
     } else {
         mask = 0x00001fff;
@@ -2659,7 +3103,7 @@ bits16 float32_to_float16( float32 a, flag ieee STATUS_PARAM)
         }
     } else {
         if (aExp > 16) {
-            float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
+            float_raise(float_flag_invalid | float_flag_inexact STATUS_VAR);
             return packFloat16(aSign, 0x1f, 0x3ff);
         }
     }
@@ -2673,8 +3117,6 @@ bits16 float32_to_float16( float32 a, flag ieee STATUS_PARAM)
     return packFloat16(aSign, aExp + 14, aSig >> 13);
 }
 
-#ifdef FLOATX80
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the double-precision floating-point value
 | `a' to the extended double-precision floating-point format.  The conversion
@@ -2686,13 +3128,14 @@ floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits64 aSig;
+    uint64_t aSig;
 
+    a = float64_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
     aSign = extractFloat64Sign( a );
     if ( aExp == 0x7FF ) {
-        if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( aExp == 0 ) {
@@ -2705,10 +3148,6 @@ floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
 
 }
 
-#endif
-
-#ifdef FLOAT128
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the double-precision floating-point value
 | `a' to the quadruple-precision floating-point format.  The conversion is
@@ -2720,13 +3159,14 @@ float128 float64_to_float128( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits64 aSig, zSig0, zSig1;
+    uint64_t aSig, zSig0, zSig1;
 
+    a = float64_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
     aSign = extractFloat64Sign( a );
     if ( aExp == 0x7FF ) {
-        if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloat128( float64ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloat128( aSign, 0x7FFF, 0, 0 );
     }
     if ( aExp == 0 ) {
@@ -2739,8 +3179,6 @@ float128 float64_to_float128( float64 a STATUS_PARAM )
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Rounds the double-precision floating-point value `a' to an integer, and
 | returns the result as a double-precision floating-point value.  The
@@ -2752,9 +3190,10 @@ float64 float64_round_to_int( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits64 lastBitMask, roundBitsMask;
+    uint64_t lastBitMask, roundBitsMask;
     int8 roundingMode;
-    bits64 z;
+    uint64_t z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aExp = extractFloat64Exp( a );
     if ( 0x433 <= aExp ) {
@@ -2764,7 +3203,7 @@ float64 float64_round_to_int( float64 a STATUS_PARAM )
         return a;
     }
     if ( aExp < 0x3FF ) {
-        if ( (bits64) ( float64_val(a)<<1 ) == 0 ) return a;
+        if ( (uint64_t) ( float64_val(a)<<1 ) == 0 ) return a;
         STATUS(float_exception_flags) |= float_flag_inexact;
         aSign = extractFloat64Sign( a );
         switch ( STATUS(float_rounding_mode) ) {
@@ -2824,7 +3263,7 @@ float64 float64_trunc_to_int( float64 a STATUS_PARAM)
 static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
 {
     int16 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig;
+    uint64_t aSig, bSig, zSig;
     int16 expDiff;
 
     aSig = extractFloat64Frac( a );
@@ -2868,7 +3307,12 @@ static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
             return a;
         }
         if ( aExp == 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloat64( zSign, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                if (aSig | bSig) {
+                    float_raise(float_flag_output_denormal STATUS_VAR);
+                }
+                return packFloat64(zSign, 0, 0);
+            }
             return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
         }
         zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
@@ -2878,7 +3322,7 @@ static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
     aSig |= LIT64( 0x2000000000000000 );
     zSig = ( aSig + bSig )<<1;
     --zExp;
-    if ( (sbits64) zSig < 0 ) {
+    if ( (int64_t) zSig < 0 ) {
         zSig = aSig + bSig;
         ++zExp;
     }
@@ -2898,7 +3342,7 @@ static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
 static float64 subFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
 {
     int16 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig;
+    uint64_t aSig, bSig, zSig;
     int16 expDiff;
 
     aSig = extractFloat64Frac( a );
@@ -2971,6 +3415,8 @@ static float64 subFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
 float64 float64_add( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     aSign = extractFloat64Sign( a );
     bSign = extractFloat64Sign( b );
@@ -2992,6 +3438,8 @@ float64 float64_add( float64 a, float64 b STATUS_PARAM )
 float64 float64_sub( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     aSign = extractFloat64Sign( a );
     bSign = extractFloat64Sign( b );
@@ -3014,7 +3462,10 @@ float64 float64_mul( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int16 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig0, zSig1;
+    uint64_t aSig, bSig, zSig0, zSig1;
+
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -3054,7 +3505,7 @@ float64 float64_mul( float64 a, float64 b STATUS_PARAM )
     bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
     mul64To128( aSig, bSig, &zSig0, &zSig1 );
     zSig0 |= ( zSig1 != 0 );
-    if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
+    if ( 0 <= (int64_t) ( zSig0<<1 ) ) {
         zSig0 <<= 1;
         --zExp;
     }
@@ -3072,9 +3523,11 @@ float64 float64_div( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int16 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig;
-    bits64 rem0, rem1;
-    bits64 term0, term1;
+    uint64_t aSig, bSig, zSig;
+    uint64_t rem0, rem1;
+    uint64_t term0, term1;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -3122,7 +3575,7 @@ float64 float64_div( float64 a, float64 b STATUS_PARAM )
     if ( ( zSig & 0x1FF ) <= 2 ) {
         mul64To128( bSig, zSig, &term0, &term1 );
         sub128( aSig, 0, term0, term1, &rem0, &rem1 );
-        while ( (sbits64) rem0 < 0 ) {
+        while ( (int64_t) rem0 < 0 ) {
             --zSig;
             add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
         }
@@ -3142,10 +3595,12 @@ float64 float64_rem( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, zSign;
     int16 aExp, bExp, expDiff;
-    bits64 aSig, bSig;
-    bits64 q, alternateASig;
-    sbits64 sigMean;
+    uint64_t aSig, bSig;
+    uint64_t q, alternateASig;
+    int64_t sigMean;
 
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
     aSign = extractFloat64Sign( a );
@@ -3205,29 +3660,250 @@ float64 float64_rem( float64 a, float64 b STATUS_PARAM )
         alternateASig = aSig;
         ++q;
         aSig -= bSig;
-    } while ( 0 <= (sbits64) aSig );
+    } while ( 0 <= (int64_t) aSig );
     sigMean = aSig + alternateASig;
     if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
         aSig = alternateASig;
     }
-    zSign = ( (sbits64) aSig < 0 );
+    zSign = ( (int64_t) aSig < 0 );
     if ( zSign ) aSig = - aSig;
     return normalizeRoundAndPackFloat64( aSign ^ zSign, bExp, aSig STATUS_VAR );
 
 }
 
 /*----------------------------------------------------------------------------
-| Returns the square root of the double-precision floating-point value `a'.
-| The operation is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
+| Returns the result of multiplying the double-precision floating-point values
+| `a' and `b' then adding 'c', with no intermediate rounding step after the
+| multiplication.  The operation is performed according to the IEC/IEEE
+| Standard for Binary Floating-Point Arithmetic 754-2008.
+| The flags argument allows the caller to select negation of the
+| addend, the intermediate product, or the final result. (The difference
+| between this and having the caller do a separate negation is that negating
+| externally will flip the sign bit on NaNs.)
+*----------------------------------------------------------------------------*/
+
+float64 float64_muladd(float64 a, float64 b, float64 c, int flags STATUS_PARAM)
+{
+    flag aSign, bSign, cSign, zSign;
+    int16 aExp, bExp, cExp, pExp, zExp, expDiff;
+    uint64_t aSig, bSig, cSig;
+    flag pInf, pZero, pSign;
+    uint64_t pSig0, pSig1, cSig0, cSig1, zSig0, zSig1;
+    int shiftcount;
+    flag signflip, infzero;
+
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
+    c = float64_squash_input_denormal(c STATUS_VAR);
+    aSig = extractFloat64Frac(a);
+    aExp = extractFloat64Exp(a);
+    aSign = extractFloat64Sign(a);
+    bSig = extractFloat64Frac(b);
+    bExp = extractFloat64Exp(b);
+    bSign = extractFloat64Sign(b);
+    cSig = extractFloat64Frac(c);
+    cExp = extractFloat64Exp(c);
+    cSign = extractFloat64Sign(c);
+
+    infzero = ((aExp == 0 && aSig == 0 && bExp == 0x7ff && bSig == 0) ||
+               (aExp == 0x7ff && aSig == 0 && bExp == 0 && bSig == 0));
+
+    /* It is implementation-defined whether the cases of (0,inf,qnan)
+     * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
+     * they return if they do), so we have to hand this information
+     * off to the target-specific pick-a-NaN routine.
+     */
+    if (((aExp == 0x7ff) && aSig) ||
+        ((bExp == 0x7ff) && bSig) ||
+        ((cExp == 0x7ff) && cSig)) {
+        return propagateFloat64MulAddNaN(a, b, c, infzero STATUS_VAR);
+    }
+
+    if (infzero) {
+        float_raise(float_flag_invalid STATUS_VAR);
+        return float64_default_nan;
+    }
 
-float64 float64_sqrt( float64 a STATUS_PARAM )
+    if (flags & float_muladd_negate_c) {
+        cSign ^= 1;
+    }
+
+    signflip = (flags & float_muladd_negate_result) ? 1 : 0;
+
+    /* Work out the sign and type of the product */
+    pSign = aSign ^ bSign;
+    if (flags & float_muladd_negate_product) {
+        pSign ^= 1;
+    }
+    pInf = (aExp == 0x7ff) || (bExp == 0x7ff);
+    pZero = ((aExp | aSig) == 0) || ((bExp | bSig) == 0);
+
+    if (cExp == 0x7ff) {
+        if (pInf && (pSign ^ cSign)) {
+            /* addition of opposite-signed infinities => InvalidOperation */
+            float_raise(float_flag_invalid STATUS_VAR);
+            return float64_default_nan;
+        }
+        /* Otherwise generate an infinity of the same sign */
+        return packFloat64(cSign ^ signflip, 0x7ff, 0);
+    }
+
+    if (pInf) {
+        return packFloat64(pSign ^ signflip, 0x7ff, 0);
+    }
+
+    if (pZero) {
+        if (cExp == 0) {
+            if (cSig == 0) {
+                /* Adding two exact zeroes */
+                if (pSign == cSign) {
+                    zSign = pSign;
+                } else if (STATUS(float_rounding_mode) == float_round_down) {
+                    zSign = 1;
+                } else {
+                    zSign = 0;
+                }
+                return packFloat64(zSign ^ signflip, 0, 0);
+            }
+            /* Exact zero plus a denorm */
+            if (STATUS(flush_to_zero)) {
+                float_raise(float_flag_output_denormal STATUS_VAR);
+                return packFloat64(cSign ^ signflip, 0, 0);
+            }
+        }
+        /* Zero plus something non-zero : just return the something */
+        return make_float64(float64_val(c) ^ ((uint64_t)signflip << 63));
+    }
+
+    if (aExp == 0) {
+        normalizeFloat64Subnormal(aSig, &aExp, &aSig);
+    }
+    if (bExp == 0) {
+        normalizeFloat64Subnormal(bSig, &bExp, &bSig);
+    }
+
+    /* Calculate the actual result a * b + c */
+
+    /* Multiply first; this is easy. */
+    /* NB: we subtract 0x3fe where float64_mul() subtracts 0x3ff
+     * because we want the true exponent, not the "one-less-than"
+     * flavour that roundAndPackFloat64() takes.
+     */
+    pExp = aExp + bExp - 0x3fe;
+    aSig = (aSig | LIT64(0x0010000000000000))<<10;
+    bSig = (bSig | LIT64(0x0010000000000000))<<11;
+    mul64To128(aSig, bSig, &pSig0, &pSig1);
+    if ((int64_t)(pSig0 << 1) >= 0) {
+        shortShift128Left(pSig0, pSig1, 1, &pSig0, &pSig1);
+        pExp--;
+    }
+
+    zSign = pSign ^ signflip;
+
+    /* Now [pSig0:pSig1] is the significand of the multiply, with the explicit
+     * bit in position 126.
+     */
+    if (cExp == 0) {
+        if (!cSig) {
+            /* Throw out the special case of c being an exact zero now */
+            shift128RightJamming(pSig0, pSig1, 64, &pSig0, &pSig1);
+            return roundAndPackFloat64(zSign, pExp - 1,
+                                       pSig1 STATUS_VAR);
+        }
+        normalizeFloat64Subnormal(cSig, &cExp, &cSig);
+    }
+
+    /* Shift cSig and add the explicit bit so [cSig0:cSig1] is the
+     * significand of the addend, with the explicit bit in position 126.
+     */
+    cSig0 = cSig << (126 - 64 - 52);
+    cSig1 = 0;
+    cSig0 |= LIT64(0x4000000000000000);
+    expDiff = pExp - cExp;
+
+    if (pSign == cSign) {
+        /* Addition */
+        if (expDiff > 0) {
+            /* scale c to match p */
+            shift128RightJamming(cSig0, cSig1, expDiff, &cSig0, &cSig1);
+            zExp = pExp;
+        } else if (expDiff < 0) {
+            /* scale p to match c */
+            shift128RightJamming(pSig0, pSig1, -expDiff, &pSig0, &pSig1);
+            zExp = cExp;
+        } else {
+            /* no scaling needed */
+            zExp = cExp;
+        }
+        /* Add significands and make sure explicit bit ends up in posn 126 */
+        add128(pSig0, pSig1, cSig0, cSig1, &zSig0, &zSig1);
+        if ((int64_t)zSig0 < 0) {
+            shift128RightJamming(zSig0, zSig1, 1, &zSig0, &zSig1);
+        } else {
+            zExp--;
+        }
+        shift128RightJamming(zSig0, zSig1, 64, &zSig0, &zSig1);
+        return roundAndPackFloat64(zSign, zExp, zSig1 STATUS_VAR);
+    } else {
+        /* Subtraction */
+        if (expDiff > 0) {
+            shift128RightJamming(cSig0, cSig1, expDiff, &cSig0, &cSig1);
+            sub128(pSig0, pSig1, cSig0, cSig1, &zSig0, &zSig1);
+            zExp = pExp;
+        } else if (expDiff < 0) {
+            shift128RightJamming(pSig0, pSig1, -expDiff, &pSig0, &pSig1);
+            sub128(cSig0, cSig1, pSig0, pSig1, &zSig0, &zSig1);
+            zExp = cExp;
+            zSign ^= 1;
+        } else {
+            zExp = pExp;
+            if (lt128(cSig0, cSig1, pSig0, pSig1)) {
+                sub128(pSig0, pSig1, cSig0, cSig1, &zSig0, &zSig1);
+            } else if (lt128(pSig0, pSig1, cSig0, cSig1)) {
+                sub128(cSig0, cSig1, pSig0, pSig1, &zSig0, &zSig1);
+                zSign ^= 1;
+            } else {
+                /* Exact zero */
+                zSign = signflip;
+                if (STATUS(float_rounding_mode) == float_round_down) {
+                    zSign ^= 1;
+                }
+                return packFloat64(zSign, 0, 0);
+            }
+        }
+        --zExp;
+        /* Do the equivalent of normalizeRoundAndPackFloat64() but
+         * starting with the significand in a pair of uint64_t.
+         */
+        if (zSig0) {
+            shiftcount = countLeadingZeros64(zSig0) - 1;
+            shortShift128Left(zSig0, zSig1, shiftcount, &zSig0, &zSig1);
+            if (zSig1) {
+                zSig0 |= 1;
+            }
+            zExp -= shiftcount;
+        } else {
+            shiftcount = countLeadingZeros64(zSig1) - 1;
+            zSig0 = zSig1 << shiftcount;
+            zExp -= (shiftcount + 64);
+        }
+        return roundAndPackFloat64(zSign, zExp, zSig0 STATUS_VAR);
+    }
+}
+
+/*----------------------------------------------------------------------------
+| Returns the square root of the double-precision floating-point value `a'.
+| The operation is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+float64 float64_sqrt( float64 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp, zExp;
-    bits64 aSig, zSig, doubleZSig;
-    bits64 rem0, rem1, term0, term1;
+    uint64_t aSig, zSig, doubleZSig;
+    uint64_t rem0, rem1, term0, term1;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -3256,7 +3932,7 @@ float64 float64_sqrt( float64 a STATUS_PARAM )
         doubleZSig = zSig<<1;
         mul64To128( zSig, zSig, &term0, &term1 );
         sub128( aSig, 0, term0, term1, &rem0, &rem1 );
-        while ( (sbits64) rem0 < 0 ) {
+        while ( (int64_t) rem0 < 0 ) {
             --zSig;
             doubleZSig -= 2;
             add128( rem0, rem1, zSig>>63, doubleZSig | 1, &rem0, &rem1 );
@@ -3276,7 +3952,8 @@ float64 float64_log2( float64 a STATUS_PARAM )
 {
     flag aSign, zSign;
     int16 aExp;
-    bits64 aSig, aSig0, aSig1, zSig, i;
+    uint64_t aSig, aSig0, aSig1, zSig, i;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -3298,7 +3975,7 @@ float64 float64_log2( float64 a STATUS_PARAM )
     aExp -= 0x3FF;
     aSig |= LIT64( 0x0010000000000000 );
     zSign = aExp < 0;
-    zSig = (bits64)aExp << 52;
+    zSig = (uint64_t)aExp << 52;
     for (i = 1LL << 51; i > 0; i >>= 1) {
         mul64To128( aSig, aSig, &aSig0, &aSig1 );
         aSig = ( aSig0 << 12 ) | ( aSig1 >> 52 );
@@ -3315,39 +3992,42 @@ float64 float64_log2( float64 a STATUS_PARAM )
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the double-precision floating-point value `a' is equal to the
-| corresponding value `b', and 0 otherwise.  The comparison is performed
+| corresponding value `b', and 0 otherwise.  The invalid exception is raised
+| if either operand is a NaN.  Otherwise, the comparison is performed
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float64_eq( float64 a, float64 b STATUS_PARAM )
 {
-    bits64 av, bv;
+    uint64_t av, bv;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
        ) {
-        if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
-            float_raise( float_flag_invalid STATUS_VAR);
-        }
+        float_raise( float_flag_invalid STATUS_VAR);
         return 0;
     }
     av = float64_val(a);
     bv = float64_val(b);
-    return ( av == bv ) || ( (bits64) ( ( av | bv )<<1 ) == 0 );
+    return ( av == bv ) || ( (uint64_t) ( ( av | bv )<<1 ) == 0 );
 
 }
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the double-precision floating-point value `a' is less than or
-| equal to the corresponding value `b', and 0 otherwise.  The comparison is
-| performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
+| equal to the corresponding value `b', and 0 otherwise.  The invalid
+| exception is raised if either operand is a NaN.  The comparison is performed
+| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float64_le( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits64 av, bv;
+    uint64_t av, bv;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
@@ -3359,22 +4039,25 @@ int float64_le( float64 a, float64 b STATUS_PARAM )
     bSign = extractFloat64Sign( b );
     av = float64_val(a);
     bv = float64_val(b);
-    if ( aSign != bSign ) return aSign || ( (bits64) ( ( av | bv )<<1 ) == 0 );
+    if ( aSign != bSign ) return aSign || ( (uint64_t) ( ( av | bv )<<1 ) == 0 );
     return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the double-precision floating-point value `a' is less than
-| the corresponding value `b', and 0 otherwise.  The comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| the corresponding value `b', and 0 otherwise.  The invalid exception is
+| raised if either operand is a NaN.  The comparison is performed according
+| to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float64_lt( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits64 av, bv;
+    uint64_t av, bv;
 
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
        ) {
@@ -3385,31 +4068,56 @@ int float64_lt( float64 a, float64 b STATUS_PARAM )
     bSign = extractFloat64Sign( b );
     av = float64_val(a);
     bv = float64_val(b);
-    if ( aSign != bSign ) return aSign && ( (bits64) ( ( av | bv )<<1 ) != 0 );
+    if ( aSign != bSign ) return aSign && ( (uint64_t) ( ( av | bv )<<1 ) != 0 );
     return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
 /*----------------------------------------------------------------------------
-| Returns 1 if the double-precision floating-point value `a' is equal to the
-| corresponding value `b', and 0 otherwise.  The invalid exception is raised
-| if either operand is a NaN.  Otherwise, the comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| Returns 1 if the double-precision floating-point values `a' and `b' cannot
+| be compared, and 0 otherwise.  The invalid exception is raised if either
+| operand is a NaN.  The comparison is performed according to the IEC/IEEE
+| Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-int float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
+int float64_unordered( float64 a, float64 b STATUS_PARAM )
 {
-    bits64 av, bv;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
        ) {
         float_raise( float_flag_invalid STATUS_VAR);
+        return 1;
+    }
+    return 0;
+}
+
+/*----------------------------------------------------------------------------
+| Returns 1 if the double-precision floating-point value `a' is equal to the
+| corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
+| exception.The comparison is performed according to the IEC/IEEE Standard
+| for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+int float64_eq_quiet( float64 a, float64 b STATUS_PARAM )
+{
+    uint64_t av, bv;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
         return 0;
     }
     av = float64_val(a);
     bv = float64_val(b);
-    return ( av == bv ) || ( (bits64) ( ( av | bv )<<1 ) == 0 );
+    return ( av == bv ) || ( (uint64_t) ( ( av | bv )<<1 ) == 0 );
 
 }
 
@@ -3423,7 +4131,9 @@ int float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
 int float64_le_quiet( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits64 av, bv;
+    uint64_t av, bv;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
@@ -3437,7 +4147,7 @@ int float64_le_quiet( float64 a, float64 b STATUS_PARAM )
     bSign = extractFloat64Sign( b );
     av = float64_val(a);
     bv = float64_val(b);
-    if ( aSign != bSign ) return aSign || ( (bits64) ( ( av | bv )<<1 ) == 0 );
+    if ( aSign != bSign ) return aSign || ( (uint64_t) ( ( av | bv )<<1 ) == 0 );
     return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
@@ -3452,7 +4162,9 @@ int float64_le_quiet( float64 a, float64 b STATUS_PARAM )
 int float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
-    bits64 av, bv;
+    uint64_t av, bv;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
          || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
@@ -3466,12 +4178,33 @@ int float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
     bSign = extractFloat64Sign( b );
     av = float64_val(a);
     bv = float64_val(b);
-    if ( aSign != bSign ) return aSign && ( (bits64) ( ( av | bv )<<1 ) != 0 );
+    if ( aSign != bSign ) return aSign && ( (uint64_t) ( ( av | bv )<<1 ) != 0 );
     return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
-#ifdef FLOATX80
+/*----------------------------------------------------------------------------
+| Returns 1 if the double-precision floating-point values `a' and `b' cannot
+| be compared, and 0 otherwise.  Quiet NaNs do not cause an exception.  The
+| comparison is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+int float64_unordered_quiet( float64 a, float64 b STATUS_PARAM )
+{
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return 1;
+    }
+    return 0;
+}
 
 /*----------------------------------------------------------------------------
 | Returns the result of converting the extended double-precision floating-
@@ -3487,12 +4220,12 @@ int32 floatx80_to_int32( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig;
+    uint64_t aSig;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
-    if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
+    if ( ( aExp == 0x7FFF ) && (uint64_t) ( aSig<<1 ) ) aSign = 0;
     shiftCount = 0x4037 - aExp;
     if ( shiftCount <= 0 ) shiftCount = 1;
     shift64RightJamming( aSig, shiftCount, &aSig );
@@ -3514,14 +4247,14 @@ int32 floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig, savedASig;
-    int32 z;
+    uint64_t aSig, savedASig;
+    int32_t z;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
     if ( 0x401E < aExp ) {
-        if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
+        if ( ( aExp == 0x7FFF ) && (uint64_t) ( aSig<<1 ) ) aSign = 0;
         goto invalid;
     }
     else if ( aExp < 0x3FFF ) {
@@ -3536,7 +4269,7 @@ int32 floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM )
     if ( ( z < 0 ) ^ aSign ) {
  invalid:
         float_raise( float_flag_invalid STATUS_VAR);
-        return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
+        return aSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
     }
     if ( ( aSig<<shiftCount ) != savedASig ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
@@ -3559,7 +4292,7 @@ int64 floatx80_to_int64( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig, aSigExtra;
+    uint64_t aSig, aSigExtra;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
@@ -3574,7 +4307,7 @@ int64 floatx80_to_int64( floatx80 a STATUS_PARAM )
                ) {
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
             }
-            return (sbits64) LIT64( 0x8000000000000000 );
+            return (int64_t) LIT64( 0x8000000000000000 );
         }
         aSigExtra = 0;
     }
@@ -3599,7 +4332,7 @@ int64 floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig;
+    uint64_t aSig;
     int64 z;
 
     aSig = extractFloatx80Frac( a );
@@ -3614,14 +4347,14 @@ int64 floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM )
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
             }
         }
-        return (sbits64) LIT64( 0x8000000000000000 );
+        return (int64_t) LIT64( 0x8000000000000000 );
     }
     else if ( aExp < 0x3FFF ) {
         if ( aExp | aSig ) STATUS(float_exception_flags) |= float_flag_inexact;
         return 0;
     }
     z = aSig>>( - shiftCount );
-    if ( (bits64) ( aSig<<( shiftCount & 63 ) ) ) {
+    if ( (uint64_t) ( aSig<<( shiftCount & 63 ) ) ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
     }
     if ( aSign ) z = - z;
@@ -3640,14 +4373,14 @@ float32 floatx80_to_float32( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 aSig;
+    uint64_t aSig;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
     if ( aExp == 0x7FFF ) {
-        if ( (bits64) ( aSig<<1 ) ) {
-            return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) );
+        if ( (uint64_t) ( aSig<<1 ) ) {
+            return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloat32( aSign, 0xFF, 0 );
     }
@@ -3668,14 +4401,14 @@ float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 aSig, zSig;
+    uint64_t aSig, zSig;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
     if ( aExp == 0x7FFF ) {
-        if ( (bits64) ( aSig<<1 ) ) {
-            return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) );
+        if ( (uint64_t) ( aSig<<1 ) ) {
+            return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloat64( aSign, 0x7FF, 0 );
     }
@@ -3685,8 +4418,6 @@ float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
 
 }
 
-#ifdef FLOAT128
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the extended double-precision floating-
 | point value `a' to the quadruple-precision floating-point format.  The
@@ -3698,21 +4429,19 @@ float128 floatx80_to_float128( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int16 aExp;
-    bits64 aSig, zSig0, zSig1;
+    uint64_t aSig, zSig0, zSig1;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
-    if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
-        return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) );
+    if ( ( aExp == 0x7FFF ) && (uint64_t) ( aSig<<1 ) ) {
+        return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
     }
     shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
     return packFloat128( aSign, aExp, zSig0, zSig1 );
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Rounds the extended double-precision floating-point value `a' to an integer,
 | and returns the result as an extended quadruple-precision floating-point
@@ -3724,27 +4453,27 @@ floatx80 floatx80_round_to_int( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 lastBitMask, roundBitsMask;
+    uint64_t lastBitMask, roundBitsMask;
     int8 roundingMode;
     floatx80 z;
 
     aExp = extractFloatx80Exp( a );
     if ( 0x403E <= aExp ) {
-        if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
+        if ( ( aExp == 0x7FFF ) && (uint64_t) ( extractFloatx80Frac( a )<<1 ) ) {
             return propagateFloatx80NaN( a, a STATUS_VAR );
         }
         return a;
     }
     if ( aExp < 0x3FFF ) {
         if (    ( aExp == 0 )
-             && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
+             && ( (uint64_t) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
             return a;
         }
         STATUS(float_exception_flags) |= float_flag_inexact;
         aSign = extractFloatx80Sign( a );
         switch ( STATUS(float_rounding_mode) ) {
          case float_round_nearest_even:
-            if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
+            if ( ( aExp == 0x3FFE ) && (uint64_t) ( extractFloatx80Frac( a )<<1 )
                ) {
                 return
                     packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
@@ -3797,7 +4526,7 @@ floatx80 floatx80_round_to_int( floatx80 a STATUS_PARAM )
 static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM)
 {
     int32 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig0, zSig1;
+    uint64_t aSig, bSig, zSig0, zSig1;
     int32 expDiff;
 
     aSig = extractFloatx80Frac( a );
@@ -3807,7 +4536,7 @@ static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
     expDiff = aExp - bExp;
     if ( 0 < expDiff ) {
         if ( aExp == 0x7FFF ) {
-            if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+            if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
             return a;
         }
         if ( bExp == 0 ) --expDiff;
@@ -3816,7 +4545,7 @@ static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
     }
     else if ( expDiff < 0 ) {
         if ( bExp == 0x7FFF ) {
-            if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+            if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
             return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
         }
         if ( aExp == 0 ) ++expDiff;
@@ -3825,7 +4554,7 @@ static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
     }
     else {
         if ( aExp == 0x7FFF ) {
-            if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
+            if ( (uint64_t) ( ( aSig | bSig )<<1 ) ) {
                 return propagateFloatx80NaN( a, b STATUS_VAR );
             }
             return a;
@@ -3840,7 +4569,7 @@ static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
         goto shiftRight1;
     }
     zSig0 = aSig + bSig;
-    if ( (sbits64) zSig0 < 0 ) goto roundAndPack;
+    if ( (int64_t) zSig0 < 0 ) goto roundAndPack;
  shiftRight1:
     shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
     zSig0 |= LIT64( 0x8000000000000000 );
@@ -3863,7 +4592,7 @@ static floatx80 addFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
 static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM )
 {
     int32 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig0, zSig1;
+    uint64_t aSig, bSig, zSig0, zSig1;
     int32 expDiff;
     floatx80 z;
 
@@ -3875,7 +4604,7 @@ static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
     if ( 0 < expDiff ) goto aExpBigger;
     if ( expDiff < 0 ) goto bExpBigger;
     if ( aExp == 0x7FFF ) {
-        if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
+        if ( (uint64_t) ( ( aSig | bSig )<<1 ) ) {
             return propagateFloatx80NaN( a, b STATUS_VAR );
         }
         float_raise( float_flag_invalid STATUS_VAR);
@@ -3893,7 +4622,7 @@ static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
     return packFloatx80( STATUS(float_rounding_mode) == float_round_down, 0, 0 );
  bExpBigger:
     if ( bExp == 0x7FFF ) {
-        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+        if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
         return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( aExp == 0 ) ++expDiff;
@@ -3905,7 +4634,7 @@ static floatx80 subFloatx80Sigs( floatx80 a, floatx80 b, flag zSign STATUS_PARAM
     goto normalizeRoundAndPack;
  aExpBigger:
     if ( aExp == 0x7FFF ) {
-        if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+        if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
         return a;
     }
     if ( bExp == 0 ) --expDiff;
@@ -3972,7 +4701,7 @@ floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int32 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig0, zSig1;
+    uint64_t aSig, bSig, zSig0, zSig1;
     floatx80 z;
 
     aSig = extractFloatx80Frac( a );
@@ -3983,15 +4712,15 @@ floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM )
     bSign = extractFloatx80Sign( b );
     zSign = aSign ^ bSign;
     if ( aExp == 0x7FFF ) {
-        if (    (bits64) ( aSig<<1 )
-             || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
+        if (    (uint64_t) ( aSig<<1 )
+             || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
             return propagateFloatx80NaN( a, b STATUS_VAR );
         }
         if ( ( bExp | bSig ) == 0 ) goto invalid;
         return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( bExp == 0x7FFF ) {
-        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+        if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
         if ( ( aExp | aSig ) == 0 ) {
  invalid:
             float_raise( float_flag_invalid STATUS_VAR);
@@ -4011,7 +4740,7 @@ floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM )
     }
     zExp = aExp + bExp - 0x3FFE;
     mul64To128( aSig, bSig, &zSig0, &zSig1 );
-    if ( 0 < (sbits64) zSig0 ) {
+    if ( 0 < (int64_t) zSig0 ) {
         shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
         --zExp;
     }
@@ -4031,8 +4760,8 @@ floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int32 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig0, zSig1;
-    bits64 rem0, rem1, rem2, term0, term1, term2;
+    uint64_t aSig, bSig, zSig0, zSig1;
+    uint64_t rem0, rem1, rem2, term0, term1, term2;
     floatx80 z;
 
     aSig = extractFloatx80Frac( a );
@@ -4043,15 +4772,15 @@ floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
     bSign = extractFloatx80Sign( b );
     zSign = aSign ^ bSign;
     if ( aExp == 0x7FFF ) {
-        if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+        if ( (uint64_t) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
         if ( bExp == 0x7FFF ) {
-            if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+            if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
             goto invalid;
         }
         return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( bExp == 0x7FFF ) {
-        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+        if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
         return packFloatx80( zSign, 0, 0 );
     }
     if ( bExp == 0 ) {
@@ -4081,15 +4810,15 @@ floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
     zSig0 = estimateDiv128To64( aSig, rem1, bSig );
     mul64To128( bSig, zSig0, &term0, &term1 );
     sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
-    while ( (sbits64) rem0 < 0 ) {
+    while ( (int64_t) rem0 < 0 ) {
         --zSig0;
         add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
     }
     zSig1 = estimateDiv128To64( rem1, 0, bSig );
-    if ( (bits64) ( zSig1<<1 ) <= 8 ) {
+    if ( (uint64_t) ( zSig1<<1 ) <= 8 ) {
         mul64To128( bSig, zSig1, &term1, &term2 );
         sub128( rem1, 0, term1, term2, &rem1, &rem2 );
-        while ( (sbits64) rem1 < 0 ) {
+        while ( (int64_t) rem1 < 0 ) {
             --zSig1;
             add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
         }
@@ -4111,8 +4840,8 @@ floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, zSign;
     int32 aExp, bExp, expDiff;
-    bits64 aSig0, aSig1, bSig;
-    bits64 q, term0, term1, alternateASig0, alternateASig1;
+    uint64_t aSig0, aSig1, bSig;
+    uint64_t q, term0, term1, alternateASig0, alternateASig1;
     floatx80 z;
 
     aSig0 = extractFloatx80Frac( a );
@@ -4121,14 +4850,14 @@ floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
     bSig = extractFloatx80Frac( b );
     bExp = extractFloatx80Exp( b );
     if ( aExp == 0x7FFF ) {
-        if (    (bits64) ( aSig0<<1 )
-             || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
+        if (    (uint64_t) ( aSig0<<1 )
+             || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) {
             return propagateFloatx80NaN( a, b STATUS_VAR );
         }
         goto invalid;
     }
     if ( bExp == 0x7FFF ) {
-        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
+        if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b STATUS_VAR );
         return a;
     }
     if ( bExp == 0 ) {
@@ -4142,7 +4871,7 @@ floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
         normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
     }
     if ( aExp == 0 ) {
-        if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
+        if ( (uint64_t) ( aSig0<<1 ) == 0 ) return a;
         normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
     }
     bSig |= LIT64( 0x8000000000000000 );
@@ -4207,15 +4936,15 @@ floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, zExp;
-    bits64 aSig0, aSig1, zSig0, zSig1, doubleZSig0;
-    bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
+    uint64_t aSig0, aSig1, zSig0, zSig1, doubleZSig0;
+    uint64_t rem0, rem1, rem2, rem3, term0, term1, term2, term3;
     floatx80 z;
 
     aSig0 = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
     if ( aExp == 0x7FFF ) {
-        if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a STATUS_VAR );
+        if ( (uint64_t) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a STATUS_VAR );
         if ( ! aSign ) return a;
         goto invalid;
     }
@@ -4238,7 +4967,7 @@ floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
     doubleZSig0 = zSig0<<1;
     mul64To128( zSig0, zSig0, &term0, &term1 );
     sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
-    while ( (sbits64) rem0 < 0 ) {
+    while ( (int64_t) rem0 < 0 ) {
         --zSig0;
         doubleZSig0 -= 2;
         add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
@@ -4250,7 +4979,7 @@ floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
         sub128( rem1, 0, term1, term2, &rem1, &rem2 );
         mul64To128( zSig1, zSig1, &term2, &term3 );
         sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
-        while ( (sbits64) rem1 < 0 ) {
+        while ( (int64_t) rem1 < 0 ) {
             --zSig1;
             shortShift128Left( 0, zSig1, 1, &term2, &term3 );
             term3 |= 1;
@@ -4268,31 +4997,28 @@ floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
 }
 
 /*----------------------------------------------------------------------------
-| Returns 1 if the extended double-precision floating-point value `a' is
-| equal to the corresponding value `b', and 0 otherwise.  The comparison is
-| performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
+| Returns 1 if the extended double-precision floating-point value `a' is equal
+| to the corresponding value `b', and 0 otherwise.  The invalid exception is
+| raised if either operand is a NaN.  Otherwise, the comparison is performed
+| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
        ) {
-        if (    floatx80_is_signaling_nan( a )
-             || floatx80_is_signaling_nan( b ) ) {
-            float_raise( float_flag_invalid STATUS_VAR);
-        }
+        float_raise( float_flag_invalid STATUS_VAR);
         return 0;
     }
     return
            ( a.low == b.low )
         && (    ( a.high == b.high )
              || (    ( a.low == 0 )
-                  && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
+                  && ( (uint16_t) ( ( a.high | b.high )<<1 ) == 0 ) )
            );
 
 }
@@ -4300,8 +5026,9 @@ int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
 /*----------------------------------------------------------------------------
 | Returns 1 if the extended double-precision floating-point value `a' is
 | less than or equal to the corresponding value `b', and 0 otherwise.  The
-| comparison is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
+| invalid exception is raised if either operand is a NaN.  The comparison is
+| performed according to the IEC/IEEE Standard for Binary Floating-Point
+| Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
@@ -4309,9 +5036,9 @@ int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
     flag aSign, bSign;
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
        ) {
         float_raise( float_flag_invalid STATUS_VAR);
         return 0;
@@ -4321,7 +5048,7 @@ int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            || (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            || (    ( ( (uint16_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  == 0 );
     }
     return
@@ -4332,9 +5059,9 @@ int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the extended double-precision floating-point value `a' is
-| less than the corresponding value `b', and 0 otherwise.  The comparison
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
+| less than the corresponding value `b', and 0 otherwise.  The invalid
+| exception is raised if either operand is a NaN.  The comparison is performed
+| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
@@ -4342,9 +5069,9 @@ int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
     flag aSign, bSign;
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
        ) {
         float_raise( float_flag_invalid STATUS_VAR);
         return 0;
@@ -4354,7 +5081,7 @@ int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            && (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            && (    ( ( (uint16_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  != 0 );
     }
     return
@@ -4364,28 +5091,50 @@ int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
 }
 
 /*----------------------------------------------------------------------------
-| Returns 1 if the extended double-precision floating-point value `a' is equal
-| to the corresponding value `b', and 0 otherwise.  The invalid exception is
-| raised if either operand is a NaN.  Otherwise, the comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| Returns 1 if the extended double-precision floating-point values `a' and `b'
+| cannot be compared, and 0 otherwise.  The invalid exception is raised if
+| either operand is a NaN.   The comparison is performed according to the
+| IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+int floatx80_unordered( floatx80 a, floatx80 b STATUS_PARAM )
+{
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        float_raise( float_flag_invalid STATUS_VAR);
+        return 1;
+    }
+    return 0;
+}
+
+/*----------------------------------------------------------------------------
+| Returns 1 if the extended double-precision floating-point value `a' is
+| equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
+| cause an exception.  The comparison is performed according to the IEC/IEEE
+| Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-int floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_eq_quiet( floatx80 a, floatx80 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
        ) {
-        float_raise( float_flag_invalid STATUS_VAR);
+        if (    floatx80_is_signaling_nan( a )
+             || floatx80_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
         return 0;
     }
     return
            ( a.low == b.low )
         && (    ( a.high == b.high )
              || (    ( a.low == 0 )
-                  && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
+                  && ( (uint16_t) ( ( a.high | b.high )<<1 ) == 0 ) )
            );
 
 }
@@ -4402,9 +5151,9 @@ int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
     flag aSign, bSign;
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
        ) {
         if (    floatx80_is_signaling_nan( a )
              || floatx80_is_signaling_nan( b ) ) {
@@ -4417,7 +5166,7 @@ int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            || (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            || (    ( ( (uint16_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  == 0 );
     }
     return
@@ -4438,9 +5187,9 @@ int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
     flag aSign, bSign;
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
          || (    ( extractFloatx80Exp( b ) == 0x7FFF )
-              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
        ) {
         if (    floatx80_is_signaling_nan( a )
              || floatx80_is_signaling_nan( b ) ) {
@@ -4453,7 +5202,7 @@ int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            && (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            && (    ( ( (uint16_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  != 0 );
     }
     return
@@ -4462,9 +5211,27 @@ int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
 
 }
 
-#endif
-
-#ifdef FLOAT128
+/*----------------------------------------------------------------------------
+| Returns 1 if the extended double-precision floating-point values `a' and `b'
+| cannot be compared, and 0 otherwise.  Quiet NaNs do not cause an exception.
+| The comparison is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+int floatx80_unordered_quiet( floatx80 a, floatx80 b STATUS_PARAM )
+{
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (uint64_t) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (uint64_t) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        if (    floatx80_is_signaling_nan( a )
+             || floatx80_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return 1;
+    }
+    return 0;
+}
 
 /*----------------------------------------------------------------------------
 | Returns the result of converting the quadruple-precision floating-point
@@ -4480,7 +5247,7 @@ int32 float128_to_int32( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig0, aSig1;
+    uint64_t aSig0, aSig1;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
@@ -4509,8 +5276,8 @@ int32 float128_to_int32_round_to_zero( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig0, aSig1, savedASig;
-    int32 z;
+    uint64_t aSig0, aSig1, savedASig;
+    int32_t z;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
@@ -4534,7 +5301,7 @@ int32 float128_to_int32_round_to_zero( float128 a STATUS_PARAM )
     if ( ( z < 0 ) ^ aSign ) {
  invalid:
         float_raise( float_flag_invalid STATUS_VAR);
-        return aSign ? (sbits32) 0x80000000 : 0x7FFFFFFF;
+        return aSign ? (int32_t) 0x80000000 : 0x7FFFFFFF;
     }
     if ( ( aSig0<<shiftCount ) != savedASig ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
@@ -4557,7 +5324,7 @@ int64 float128_to_int64( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig0, aSig1;
+    uint64_t aSig0, aSig1;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
@@ -4575,7 +5342,7 @@ int64 float128_to_int64( float128 a STATUS_PARAM )
                ) {
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
             }
-            return (sbits64) LIT64( 0x8000000000000000 );
+            return (int64_t) LIT64( 0x8000000000000000 );
         }
         shortShift128Left( aSig0, aSig1, - shiftCount, &aSig0, &aSig1 );
     }
@@ -4600,7 +5367,7 @@ int64 float128_to_int64_round_to_zero( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, shiftCount;
-    bits64 aSig0, aSig1;
+    uint64_t aSig0, aSig1;
     int64 z;
 
     aSig1 = extractFloat128Frac1( a );
@@ -4622,10 +5389,10 @@ int64 float128_to_int64_round_to_zero( float128 a STATUS_PARAM )
                     return LIT64( 0x7FFFFFFFFFFFFFFF );
                 }
             }
-            return (sbits64) LIT64( 0x8000000000000000 );
+            return (int64_t) LIT64( 0x8000000000000000 );
         }
         z = ( aSig0<<shiftCount ) | ( aSig1>>( ( - shiftCount ) & 63 ) );
-        if ( (bits64) ( aSig1<<shiftCount ) ) {
+        if ( (uint64_t) ( aSig1<<shiftCount ) ) {
             STATUS(float_exception_flags) |= float_flag_inexact;
         }
     }
@@ -4638,7 +5405,7 @@ int64 float128_to_int64_round_to_zero( float128 a STATUS_PARAM )
         }
         z = aSig0>>( - shiftCount );
         if (    aSig1
-             || ( shiftCount && (bits64) ( aSig0<<( shiftCount & 63 ) ) ) ) {
+             || ( shiftCount && (uint64_t) ( aSig0<<( shiftCount & 63 ) ) ) ) {
             STATUS(float_exception_flags) |= float_flag_inexact;
         }
     }
@@ -4658,8 +5425,8 @@ float32 float128_to_float32( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 aSig0, aSig1;
-    bits32 zSig;
+    uint64_t aSig0, aSig1;
+    uint32_t zSig;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
@@ -4667,7 +5434,7 @@ float32 float128_to_float32( float128 a STATUS_PARAM )
     aSign = extractFloat128Sign( a );
     if ( aExp == 0x7FFF ) {
         if ( aSig0 | aSig1 ) {
-            return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR ) );
+            return commonNaNToFloat32( float128ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloat32( aSign, 0xFF, 0 );
     }
@@ -4693,7 +5460,7 @@ float64 float128_to_float64( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 aSig0, aSig1;
+    uint64_t aSig0, aSig1;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
@@ -4701,7 +5468,7 @@ float64 float128_to_float64( float128 a STATUS_PARAM )
     aSign = extractFloat128Sign( a );
     if ( aExp == 0x7FFF ) {
         if ( aSig0 | aSig1 ) {
-            return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR ) );
+            return commonNaNToFloat64( float128ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloat64( aSign, 0x7FF, 0 );
     }
@@ -4715,8 +5482,6 @@ float64 float128_to_float64( float128 a STATUS_PARAM )
 
 }
 
-#ifdef FLOATX80
-
 /*----------------------------------------------------------------------------
 | Returns the result of converting the quadruple-precision floating-point
 | value `a' to the extended double-precision floating-point format.  The
@@ -4728,7 +5493,7 @@ floatx80 float128_to_floatx80( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 aSig0, aSig1;
+    uint64_t aSig0, aSig1;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
@@ -4736,7 +5501,7 @@ floatx80 float128_to_floatx80( float128 a STATUS_PARAM )
     aSign = extractFloat128Sign( a );
     if ( aExp == 0x7FFF ) {
         if ( aSig0 | aSig1 ) {
-            return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR ) );
+            return commonNaNToFloatx80( float128ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
@@ -4752,8 +5517,6 @@ floatx80 float128_to_floatx80( float128 a STATUS_PARAM )
 
 }
 
-#endif
-
 /*----------------------------------------------------------------------------
 | Rounds the quadruple-precision floating-point value `a' to an integer, and
 | returns the result as a quadruple-precision floating-point value.  The
@@ -4765,7 +5528,7 @@ float128 float128_round_to_int( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp;
-    bits64 lastBitMask, roundBitsMask;
+    uint64_t lastBitMask, roundBitsMask;
     int8 roundingMode;
     float128 z;
 
@@ -4790,9 +5553,9 @@ float128 float128_round_to_int( float128 a STATUS_PARAM )
                 if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
             }
             else {
-                if ( (sbits64) z.low < 0 ) {
+                if ( (int64_t) z.low < 0 ) {
                     ++z.high;
-                    if ( (bits64) ( z.low<<1 ) == 0 ) z.high &= ~1;
+                    if ( (uint64_t) ( z.low<<1 ) == 0 ) z.high &= ~1;
                 }
             }
         }
@@ -4806,7 +5569,7 @@ float128 float128_round_to_int( float128 a STATUS_PARAM )
     }
     else {
         if ( aExp < 0x3FFF ) {
-            if ( ( ( (bits64) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
+            if ( ( ( (uint64_t) ( a.high<<1 ) ) | a.low ) == 0 ) return a;
             STATUS(float_exception_flags) |= float_flag_inexact;
             aSign = extractFloat128Sign( a );
             switch ( STATUS(float_rounding_mode) ) {
@@ -4868,7 +5631,7 @@ float128 float128_round_to_int( float128 a STATUS_PARAM )
 static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
 {
     int32 aExp, bExp, zExp;
-    bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
+    uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
     int32 expDiff;
 
     aSig1 = extractFloat128Frac1( a );
@@ -4917,7 +5680,12 @@ static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM
         }
         add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
         if ( aExp == 0 ) {
-            if ( STATUS(flush_to_zero) ) return packFloat128( zSign, 0, 0, 0 );
+            if (STATUS(flush_to_zero)) {
+                if (zSig0 | zSig1) {
+                    float_raise(float_flag_output_denormal STATUS_VAR);
+                }
+                return packFloat128(zSign, 0, 0, 0);
+            }
             return packFloat128( zSign, 0, zSig0, zSig1 );
         }
         zSig2 = 0;
@@ -4949,7 +5717,7 @@ static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM
 static float128 subFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM)
 {
     int32 aExp, bExp, zExp;
-    bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
+    uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
     int32 expDiff;
     float128 z;
 
@@ -5074,7 +5842,7 @@ float128 float128_mul( float128 a, float128 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
     int32 aExp, bExp, zExp;
-    bits64 aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
+    uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2, zSig3;
     float128 z;
 
     aSig1 = extractFloat128Frac1( a );
@@ -5138,8 +5906,8 @@ float128 float128_div( float128 a, float128 b STATUS_PARAM )
 {
     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;
+    uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
+    uint64_t rem0, rem1, rem2, rem3, term0, term1, term2, term3;
     float128 z;
 
     aSig1 = extractFloat128Frac1( a );
@@ -5193,7 +5961,7 @@ float128 float128_div( float128 a, float128 b STATUS_PARAM )
     zSig0 = estimateDiv128To64( aSig0, aSig1, bSig0 );
     mul128By64To192( bSig0, bSig1, zSig0, &term0, &term1, &term2 );
     sub192( aSig0, aSig1, 0, term0, term1, term2, &rem0, &rem1, &rem2 );
-    while ( (sbits64) rem0 < 0 ) {
+    while ( (int64_t) rem0 < 0 ) {
         --zSig0;
         add192( rem0, rem1, rem2, 0, bSig0, bSig1, &rem0, &rem1, &rem2 );
     }
@@ -5201,7 +5969,7 @@ float128 float128_div( float128 a, float128 b STATUS_PARAM )
     if ( ( zSig1 & 0x3FFF ) <= 4 ) {
         mul128By64To192( bSig0, bSig1, zSig1, &term1, &term2, &term3 );
         sub192( rem1, rem2, 0, term1, term2, term3, &rem1, &rem2, &rem3 );
-        while ( (sbits64) rem1 < 0 ) {
+        while ( (int64_t) rem1 < 0 ) {
             --zSig1;
             add192( rem1, rem2, rem3, 0, bSig0, bSig1, &rem1, &rem2, &rem3 );
         }
@@ -5222,9 +5990,9 @@ float128 float128_rem( float128 a, float128 b STATUS_PARAM )
 {
     flag aSign, zSign;
     int32 aExp, bExp, expDiff;
-    bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
-    bits64 allZero, alternateASig0, alternateASig1, sigMean1;
-    sbits64 sigMean0;
+    uint64_t aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
+    uint64_t allZero, alternateASig0, alternateASig1, sigMean1;
+    int64_t sigMean0;
     float128 z;
 
     aSig1 = extractFloat128Frac1( a );
@@ -5306,15 +6074,15 @@ float128 float128_rem( float128 a, float128 b STATUS_PARAM )
         alternateASig1 = aSig1;
         ++q;
         sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
-    } while ( 0 <= (sbits64) aSig0 );
+    } while ( 0 <= (int64_t) aSig0 );
     add128(
-        aSig0, aSig1, alternateASig0, alternateASig1, (bits64 *)&sigMean0, &sigMean1 );
+        aSig0, aSig1, alternateASig0, alternateASig1, (uint64_t *)&sigMean0, &sigMean1 );
     if (    ( sigMean0 < 0 )
          || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
         aSig0 = alternateASig0;
         aSig1 = alternateASig1;
     }
-    zSign = ( (sbits64) aSig0 < 0 );
+    zSign = ( (int64_t) aSig0 < 0 );
     if ( zSign ) sub128( 0, 0, aSig0, aSig1, &aSig0, &aSig1 );
     return
         normalizeRoundAndPackFloat128( aSign ^ zSign, bExp - 4, aSig0, aSig1 STATUS_VAR );
@@ -5331,8 +6099,8 @@ float128 float128_sqrt( float128 a STATUS_PARAM )
 {
     flag aSign;
     int32 aExp, zExp;
-    bits64 aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
-    bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
+    uint64_t aSig0, aSig1, zSig0, zSig1, zSig2, doubleZSig0;
+    uint64_t rem0, rem1, rem2, rem3, term0, term1, term2, term3;
     float128 z;
 
     aSig1 = extractFloat128Frac1( a );
@@ -5364,7 +6132,7 @@ float128 float128_sqrt( float128 a STATUS_PARAM )
     doubleZSig0 = zSig0<<1;
     mul64To128( zSig0, zSig0, &term0, &term1 );
     sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
-    while ( (sbits64) rem0 < 0 ) {
+    while ( (int64_t) rem0 < 0 ) {
         --zSig0;
         doubleZSig0 -= 2;
         add128( rem0, rem1, zSig0>>63, doubleZSig0 | 1, &rem0, &rem1 );
@@ -5376,7 +6144,7 @@ float128 float128_sqrt( float128 a STATUS_PARAM )
         sub128( rem1, 0, term1, term2, &rem1, &rem2 );
         mul64To128( zSig1, zSig1, &term2, &term3 );
         sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
-        while ( (sbits64) rem1 < 0 ) {
+        while ( (int64_t) rem1 < 0 ) {
             --zSig1;
             shortShift128Left( 0, zSig1, 1, &term2, &term3 );
             term3 |= 1;
@@ -5392,7 +6160,8 @@ float128 float128_sqrt( float128 a STATUS_PARAM )
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
-| the corresponding value `b', and 0 otherwise.  The comparison is performed
+| the corresponding value `b', and 0 otherwise.  The invalid exception is
+| raised if either operand is a NaN.  Otherwise, the comparison is performed
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
@@ -5404,26 +6173,23 @@ int float128_eq( float128 a, float128 b STATUS_PARAM )
          || (    ( extractFloat128Exp( b ) == 0x7FFF )
               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
        ) {
-        if (    float128_is_signaling_nan( a )
-             || float128_is_signaling_nan( b ) ) {
-            float_raise( float_flag_invalid STATUS_VAR);
-        }
+        float_raise( float_flag_invalid STATUS_VAR);
         return 0;
     }
     return
            ( a.low == b.low )
         && (    ( a.high == b.high )
              || (    ( a.low == 0 )
-                  && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
+                  && ( (uint64_t) ( ( a.high | b.high )<<1 ) == 0 ) )
            );
 
 }
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the quadruple-precision floating-point value `a' is less than
-| or equal to the corresponding value `b', and 0 otherwise.  The comparison
-| is performed according to the IEC/IEEE Standard for Binary Floating-Point
-| Arithmetic.
+| or equal to the corresponding value `b', and 0 otherwise.  The invalid
+| exception is raised if either operand is a NaN.  The comparison is performed
+| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float128_le( float128 a, float128 b STATUS_PARAM )
@@ -5443,7 +6209,7 @@ int float128_le( float128 a, float128 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            || (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            || (    ( ( (uint64_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  == 0 );
     }
     return
@@ -5454,8 +6220,9 @@ int float128_le( float128 a, float128 b STATUS_PARAM )
 
 /*----------------------------------------------------------------------------
 | Returns 1 if the quadruple-precision floating-point value `a' is less than
-| the corresponding value `b', and 0 otherwise.  The comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| the corresponding value `b', and 0 otherwise.  The invalid exception is
+| raised if either operand is a NaN.  The comparison is performed according
+| to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
 int float128_lt( float128 a, float128 b STATUS_PARAM )
@@ -5475,7 +6242,7 @@ int float128_lt( float128 a, float128 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            && (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            && (    ( ( (uint64_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  != 0 );
     }
     return
@@ -5484,14 +6251,34 @@ int float128_lt( float128 a, float128 b STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns 1 if the quadruple-precision floating-point values `a' and `b' cannot
+| be compared, and 0 otherwise.  The invalid exception is raised if either
+| operand is a NaN. The comparison is performed according to the IEC/IEEE
+| Standard for Binary Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+int float128_unordered( float128 a, float128 b STATUS_PARAM )
+{
+    if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
+              && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
+         || (    ( extractFloat128Exp( b ) == 0x7FFF )
+              && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
+       ) {
+        float_raise( float_flag_invalid STATUS_VAR);
+        return 1;
+    }
+    return 0;
+}
+
 /*----------------------------------------------------------------------------
 | Returns 1 if the quadruple-precision floating-point value `a' is equal to
-| the corresponding value `b', and 0 otherwise.  The invalid exception is
-| raised if either operand is a NaN.  Otherwise, the comparison is performed
-| according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
+| the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
+| exception.  The comparison is performed according to the IEC/IEEE Standard
+| for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-int float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
+int float128_eq_quiet( float128 a, float128 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
@@ -5499,14 +6286,17 @@ int float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
          || (    ( extractFloat128Exp( b ) == 0x7FFF )
               && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
        ) {
-        float_raise( float_flag_invalid STATUS_VAR);
+        if (    float128_is_signaling_nan( a )
+             || float128_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
         return 0;
     }
     return
            ( a.low == b.low )
         && (    ( a.high == b.high )
              || (    ( a.low == 0 )
-                  && ( (bits64) ( ( a.high | b.high )<<1 ) == 0 ) )
+                  && ( (uint64_t) ( ( a.high | b.high )<<1 ) == 0 ) )
            );
 
 }
@@ -5538,7 +6328,7 @@ int float128_le_quiet( float128 a, float128 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            || (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            || (    ( ( (uint64_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  == 0 );
     }
     return
@@ -5574,7 +6364,7 @@ int float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
     if ( aSign != bSign ) {
         return
                aSign
-            && (    ( ( (bits64) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+            && (    ( ( (uint64_t) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
                  != 0 );
     }
     return
@@ -5583,23 +6373,44 @@ int float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
 
 }
 
-#endif
+/*----------------------------------------------------------------------------
+| Returns 1 if the quadruple-precision floating-point values `a' and `b' cannot
+| be compared, and 0 otherwise.  Quiet NaNs do not cause an exception.  The
+| comparison is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+int float128_unordered_quiet( float128 a, float128 b STATUS_PARAM )
+{
+    if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
+              && ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) )
+         || (    ( extractFloat128Exp( b ) == 0x7FFF )
+              && ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )
+       ) {
+        if (    float128_is_signaling_nan( a )
+             || float128_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return 1;
+    }
+    return 0;
+}
 
 /* misc functions */
-float32 uint32_to_float32( unsigned int a STATUS_PARAM )
+float32 uint32_to_float32( uint32 a STATUS_PARAM )
 {
     return int64_to_float32(a STATUS_VAR);
 }
 
-float64 uint32_to_float64( unsigned int a STATUS_PARAM )
+float64 uint32_to_float64( uint32 a STATUS_PARAM )
 {
     return int64_to_float64(a STATUS_VAR);
 }
 
-unsigned int float32_to_uint32( float32 a STATUS_PARAM )
+uint32 float32_to_uint32( float32 a STATUS_PARAM )
 {
     int64_t v;
-    unsigned int res;
+    uint32 res;
 
     v = float32_to_int64(a STATUS_VAR);
     if (v < 0) {
@@ -5614,10 +6425,10 @@ unsigned int float32_to_uint32( float32 a STATUS_PARAM )
     return res;
 }
 
-unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
+uint32 float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
 {
     int64_t v;
-    unsigned int res;
+    uint32 res;
 
     v = float32_to_int64_round_to_zero(a STATUS_VAR);
     if (v < 0) {
@@ -5632,10 +6443,28 @@ unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
     return res;
 }
 
-unsigned int float64_to_uint32( float64 a STATUS_PARAM )
+uint_fast16_t float32_to_uint16_round_to_zero(float32 a STATUS_PARAM)
+{
+    int64_t v;
+    uint_fast16_t res;
+
+    v = float32_to_int64_round_to_zero(a STATUS_VAR);
+    if (v < 0) {
+        res = 0;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else if (v > 0xffff) {
+        res = 0xffff;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else {
+        res = v;
+    }
+    return res;
+}
+
+uint32 float64_to_uint32( float64 a STATUS_PARAM )
 {
     int64_t v;
-    unsigned int res;
+    uint32 res;
 
     v = float64_to_int64(a STATUS_VAR);
     if (v < 0) {
@@ -5650,10 +6479,10 @@ unsigned int float64_to_uint32( float64 a STATUS_PARAM )
     return res;
 }
 
-unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
+uint32 float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
 {
     int64_t v;
-    unsigned int res;
+    uint32 res;
 
     v = float64_to_int64_round_to_zero(a STATUS_VAR);
     if (v < 0) {
@@ -5668,6 +6497,24 @@ unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
     return res;
 }
 
+uint_fast16_t float64_to_uint16_round_to_zero(float64 a STATUS_PARAM)
+{
+    int64_t v;
+    uint_fast16_t res;
+
+    v = float64_to_int64_round_to_zero(a STATUS_VAR);
+    if (v < 0) {
+        res = 0;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else if (v > 0xffff) {
+        res = 0xffff;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else {
+        res = v;
+    }
+    return res;
+}
+
 /* FIXME: This looks broken.  */
 uint64_t float64_to_uint64 (float64 a STATUS_PARAM)
 {
@@ -5696,7 +6543,9 @@ INLINE int float ## s ## _compare_internal( float ## s a, float ## s b,      \
                                       int is_quiet STATUS_PARAM )            \
 {                                                                            \
     flag aSign, bSign;                                                       \
-    bits ## s av, bv;                                                        \
+    uint ## s ## _t av, bv;                                                  \
+    a = float ## s ## _squash_input_denormal(a STATUS_VAR);                  \
+    b = float ## s ## _squash_input_denormal(b STATUS_VAR);                  \
                                                                              \
     if (( ( extractFloat ## s ## Exp( a ) == nan_exp ) &&                    \
          extractFloat ## s ## Frac( a ) ) ||                                 \
@@ -5714,7 +6563,7 @@ INLINE int float ## s ## _compare_internal( float ## s a, float ## s b,      \
     av = float ## s ## _val(a);                                              \
     bv = float ## s ## _val(b);                                              \
     if ( aSign != bSign ) {                                                  \
-        if ( (bits ## s) ( ( av | bv )<<1 ) == 0 ) {                         \
+        if ( (uint ## s ## _t) ( ( av | bv )<<1 ) == 0 ) {                   \
             /* zero case */                                                  \
             return float_relation_equal;                                     \
         } else {                                                             \
@@ -5742,6 +6591,52 @@ int float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM )  \
 COMPARE(32, 0xff)
 COMPARE(64, 0x7ff)
 
+INLINE int floatx80_compare_internal( floatx80 a, floatx80 b,
+                                      int is_quiet STATUS_PARAM )
+{
+    flag aSign, bSign;
+
+    if (( ( extractFloatx80Exp( a ) == 0x7fff ) &&
+          ( extractFloatx80Frac( a )<<1 ) ) ||
+        ( ( extractFloatx80Exp( b ) == 0x7fff ) &&
+          ( extractFloatx80Frac( b )<<1 ) )) {
+        if (!is_quiet ||
+            floatx80_is_signaling_nan( a ) ||
+            floatx80_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return float_relation_unordered;
+    }
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign != bSign ) {
+
+        if ( ( ( (uint16_t) ( ( a.high | b.high ) << 1 ) ) == 0) &&
+             ( ( a.low | b.low ) == 0 ) ) {
+            /* zero case */
+            return float_relation_equal;
+        } else {
+            return 1 - (2 * aSign);
+        }
+    } else {
+        if (a.low == b.low && a.high == b.high) {
+            return float_relation_equal;
+        } else {
+            return 1 - 2 * (aSign ^ ( lt128( a.high, a.low, b.high, b.low ) ));
+        }
+    }
+}
+
+int floatx80_compare( floatx80 a, floatx80 b STATUS_PARAM )
+{
+    return floatx80_compare_internal(a, b, 0 STATUS_VAR);
+}
+
+int floatx80_compare_quiet( floatx80 a, floatx80 b STATUS_PARAM )
+{
+    return floatx80_compare_internal(a, b, 1 STATUS_VAR);
+}
+
 INLINE int float128_compare_internal( float128 a, float128 b,
                                       int is_quiet STATUS_PARAM )
 {
@@ -5786,18 +6681,71 @@ int float128_compare_quiet( float128 a, float128 b STATUS_PARAM )
     return float128_compare_internal(a, b, 1 STATUS_VAR);
 }
 
+/* min() and max() functions. These can't be implemented as
+ * 'compare and pick one input' because that would mishandle
+ * NaNs and +0 vs -0.
+ */
+#define MINMAX(s, nan_exp)                                              \
+INLINE float ## s float ## s ## _minmax(float ## s a, float ## s b,     \
+                                        int ismin STATUS_PARAM )        \
+{                                                                       \
+    flag aSign, bSign;                                                  \
+    uint ## s ## _t av, bv;                                             \
+    a = float ## s ## _squash_input_denormal(a STATUS_VAR);             \
+    b = float ## s ## _squash_input_denormal(b STATUS_VAR);             \
+    if (float ## s ## _is_any_nan(a) ||                                 \
+        float ## s ## _is_any_nan(b)) {                                 \
+        return propagateFloat ## s ## NaN(a, b STATUS_VAR);             \
+    }                                                                   \
+    aSign = extractFloat ## s ## Sign(a);                               \
+    bSign = extractFloat ## s ## Sign(b);                               \
+    av = float ## s ## _val(a);                                         \
+    bv = float ## s ## _val(b);                                         \
+    if (aSign != bSign) {                                               \
+        if (ismin) {                                                    \
+            return aSign ? a : b;                                       \
+        } else {                                                        \
+            return aSign ? b : a;                                       \
+        }                                                               \
+    } else {                                                            \
+        if (ismin) {                                                    \
+            return (aSign ^ (av < bv)) ? a : b;                         \
+        } else {                                                        \
+            return (aSign ^ (av < bv)) ? b : a;                         \
+        }                                                               \
+    }                                                                   \
+}                                                                       \
+                                                                        \
+float ## s float ## s ## _min(float ## s a, float ## s b STATUS_PARAM)  \
+{                                                                       \
+    return float ## s ## _minmax(a, b, 1 STATUS_VAR);                   \
+}                                                                       \
+                                                                        \
+float ## s float ## s ## _max(float ## s a, float ## s b STATUS_PARAM)  \
+{                                                                       \
+    return float ## s ## _minmax(a, b, 0 STATUS_VAR);                   \
+}
+
+MINMAX(32, 0xff)
+MINMAX(64, 0x7ff)
+
+
 /* Multiply A by 2 raised to the power N.  */
 float32 float32_scalbn( float32 a, int n STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits32 aSig;
+    int16_t aExp;
+    uint32_t aSig;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
 
     if ( aExp == 0xFF ) {
+        if ( aSig ) {
+            return propagateFloat32NaN( a, a STATUS_VAR );
+        }
         return a;
     }
     if ( aExp != 0 )
@@ -5805,6 +6753,12 @@ float32 float32_scalbn( float32 a, int n STATUS_PARAM )
     else if ( aSig == 0 )
         return a;
 
+    if (n > 0x200) {
+        n = 0x200;
+    } else if (n < -0x200) {
+        n = -0x200;
+    }
+
     aExp += n - 1;
     aSig <<= 7;
     return normalizeRoundAndPackFloat32( aSign, aExp, aSig STATUS_VAR );
@@ -5813,14 +6767,18 @@ float32 float32_scalbn( float32 a, int n STATUS_PARAM )
 float64 float64_scalbn( float64 a, int n STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 aSig;
+    int16_t aExp;
+    uint64_t aSig;
 
+    a = float64_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
     aSign = extractFloat64Sign( a );
 
     if ( aExp == 0x7FF ) {
+        if ( aSig ) {
+            return propagateFloat64NaN( a, a STATUS_VAR );
+        }
         return a;
     }
     if ( aExp != 0 )
@@ -5828,46 +6786,62 @@ float64 float64_scalbn( float64 a, int n STATUS_PARAM )
     else if ( aSig == 0 )
         return a;
 
+    if (n > 0x1000) {
+        n = 0x1000;
+    } else if (n < -0x1000) {
+        n = -0x1000;
+    }
+
     aExp += n - 1;
     aSig <<= 10;
     return normalizeRoundAndPackFloat64( aSign, aExp, aSig STATUS_VAR );
 }
 
-#ifdef FLOATX80
 floatx80 floatx80_scalbn( floatx80 a, int n STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 aSig;
+    int32_t aExp;
+    uint64_t aSig;
 
     aSig = extractFloatx80Frac( a );
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
 
-    if ( aExp == 0x7FF ) {
+    if ( aExp == 0x7FFF ) {
+        if ( aSig<<1 ) {
+            return propagateFloatx80NaN( a, a STATUS_VAR );
+        }
         return a;
     }
+
     if (aExp == 0 && aSig == 0)
         return a;
 
+    if (n > 0x10000) {
+        n = 0x10000;
+    } else if (n < -0x10000) {
+        n = -0x10000;
+    }
+
     aExp += n;
     return normalizeRoundAndPackFloatx80( STATUS(floatx80_rounding_precision),
                                           aSign, aExp, aSig, 0 STATUS_VAR );
 }
-#endif
 
-#ifdef FLOAT128
 float128 float128_scalbn( float128 a, int n STATUS_PARAM )
 {
     flag aSign;
-    int32 aExp;
-    bits64 aSig0, aSig1;
+    int32_t aExp;
+    uint64_t aSig0, aSig1;
 
     aSig1 = extractFloat128Frac1( a );
     aSig0 = extractFloat128Frac0( a );
     aExp = extractFloat128Exp( a );
     aSign = extractFloat128Sign( a );
     if ( aExp == 0x7FFF ) {
+        if ( aSig0 | aSig1 ) {
+            return propagateFloat128NaN( a, a STATUS_VAR );
+        }
         return a;
     }
     if ( aExp != 0 )
@@ -5875,9 +6849,14 @@ float128 float128_scalbn( float128 a, int n STATUS_PARAM )
     else if ( aSig0 == 0 && aSig1 == 0 )
         return a;
 
+    if (n > 0x10000) {
+        n = 0x10000;
+    } else if (n < -0x10000) {
+        n = -0x10000;
+    }
+
     aExp += n - 1;
     return normalizeRoundAndPackFloat128( aSign, aExp, aSig0, aSig1
                                           STATUS_VAR );
 
 }
-#endif