]> git.proxmox.com Git - qemu.git/blobdiff - fpu/softfloat.c
qdev-properties-system.c: Allow vlan or netdev for -device, not both
[qemu.git] / fpu / softfloat.c
index 6f5b05d5fed593bbe924e3fa2343102faf80e782..7ba51b6f3c11bd08f1a41478afdf89d23bfcb4a3 100644 (file)
@@ -1,3 +1,8 @@
+/*
+ * QEMU float support
+ *
+ * Derived from SoftFloat.
+ */
 
 /*============================================================================
 
@@ -30,9 +35,12 @@ 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.  */
-#include "softfloat.h"
+/* softfloat (and in particular the code in softfloat-specialize.h) is
+ * target-dependent and needs the TARGET_* macros.
+ */
+#include "config.h"
+
+#include "fpu/softfloat.h"
 
 /*----------------------------------------------------------------------------
 | Primitive arithmetic functions, including multi-word arithmetic, and
@@ -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 int_fast16_t 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;
@@ -185,7 +218,7 @@ INLINE bits32 extractFloat32Frac( float32 a )
 | Returns the exponent bits of the single-precision floating-point value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE int16 extractFloat32Exp( float32 a )
+INLINE int_fast16_t extractFloat32Exp(float32 a)
 {
 
     return ( float32_val(a)>>23 ) & 0xFF;
@@ -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, int_fast16_t *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, int_fast16_t 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, int_fast16_t 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, int_fast16_t 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 );
@@ -349,7 +400,7 @@ INLINE bits64 extractFloat64Frac( float64 a )
 | Returns the exponent bits of the double-precision floating-point value `a'.
 *----------------------------------------------------------------------------*/
 
-INLINE int16 extractFloat64Exp( float64 a )
+INLINE int_fast16_t extractFloat64Exp(float64 a)
 {
 
     return ( float64_val(a)>>52 ) & 0x7FF;
@@ -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, int_fast16_t *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, int_fast16_t zExp, uint64_t zSig)
 {
 
     return make_float64(
-        ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig);
+        ( ( (uint64_t) zSign )<<63 ) + ( ( (uint64_t) zExp )<<52 ) + zSig);
 
 }
 
@@ -426,11 +492,11 @@ 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, int_fast16_t zExp, uint64_t zSig STATUS_PARAM)
 {
     int8 roundingMode;
     flag roundNearestEven;
-    int16 roundIncrement, roundBits;
+    int_fast16_t roundIncrement, roundBits;
     flag isTiny;
 
     roundingMode = STATUS(float_rounding_mode);
@@ -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, int_fast16_t 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
@@ -1179,7 +1238,7 @@ float32 uint64_to_float32( uint64 a STATUS_PARAM )
     if ( a == 0 ) return float32_zero;
     shiftCount = countLeadingZeros64( a ) - 40;
     if ( 0 <= shiftCount ) {
-        return packFloat32( 1 > 0, 0x95 - shiftCount, a<<shiftCount );
+        return packFloat32(0, 0x95 - shiftCount, a<<shiftCount);
     }
     else {
         shiftCount += 7;
@@ -1189,7 +1248,7 @@ float32 uint64_to_float32( uint64 a STATUS_PARAM )
         else {
             a <<= shiftCount;
         }
-        return roundAndPackFloat32( 1 > 0, 0x9C - shiftCount, a STATUS_VAR );
+        return roundAndPackFloat32(0, 0x9C - shiftCount, a STATUS_VAR);
     }
 }
 
@@ -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 );
@@ -1212,15 +1271,20 @@ float64 int64_to_float64( int64 a STATUS_PARAM )
 
 }
 
-float64 uint64_to_float64( uint64 a STATUS_PARAM )
+float64 uint64_to_float64(uint64 a STATUS_PARAM)
 {
-    if ( a == 0 ) return float64_zero;
-    return normalizeRoundAndPackFloat64( 0, 0x43C, a STATUS_VAR );
+    int exp =  0x43C;
 
+    if (a == 0) {
+        return float64_zero;
+    }
+    if ((int64_t)a < 0) {
+        shift64RightJamming(a, 1, &a);
+        exp += 1;
+    }
+    return normalizeRoundAndPackFloat64(0, exp, a STATUS_VAR);
 }
 
-#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 +1306,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 +1318,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,7 +1339,13 @@ float128 int64_to_float128( int64 a STATUS_PARAM )
 
 }
 
-#endif
+float128 uint64_to_float128(uint64 a STATUS_PARAM)
+{
+    if (a == 0) {
+        return float128_zero;
+    }
+    return normalizeRoundAndPackFloat128(0, 0x406E, a, 0 STATUS_VAR);
+}
 
 /*----------------------------------------------------------------------------
 | Returns the result of converting the single-precision floating-point value
@@ -1294,10 +1360,11 @@ float128 int64_to_float128( int64 a STATUS_PARAM )
 int32 float32_to_int32( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits32 aSig;
-    bits64 aSig64;
+    int_fast16_t aExp, shiftCount;
+    uint32_t aSig;
+    uint64_t aSig64;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
@@ -1324,9 +1391,10 @@ int32 float32_to_int32( float32 a STATUS_PARAM )
 int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits32 aSig;
-    int32 z;
+    int_fast16_t aExp, shiftCount;
+    uint32_t aSig;
+    int32_t z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1337,7 +1405,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 +1413,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;
@@ -1363,11 +1431,11 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
 | returned.
 *----------------------------------------------------------------------------*/
 
-int16 float32_to_int16_round_to_zero( float32 a STATUS_PARAM )
+int_fast16_t float32_to_int16_round_to_zero(float32 a STATUS_PARAM)
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits32 aSig;
+    int_fast16_t aExp, shiftCount;
+    uint32_t aSig;
     int32 z;
 
     aSig = extractFloat32Frac( a );
@@ -1381,7 +1449,7 @@ int16 float32_to_int16_round_to_zero( float32 a STATUS_PARAM )
                 return 0x7FFF;
             }
         }
-        return (sbits32) 0xffff8000;
+        return (int32_t) 0xffff8000;
     }
     else if ( aExp <= 0x7E ) {
         if ( aExp | aSig ) {
@@ -1392,7 +1460,7 @@ int16 float32_to_int16_round_to_zero( float32 a STATUS_PARAM )
     shiftCount -= 0x10;
     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 ) {
@@ -1415,9 +1483,10 @@ int16 float32_to_int16_round_to_zero( float32 a STATUS_PARAM )
 int64 float32_to_int64( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits32 aSig;
-    bits64 aSig64, aSigExtra;
+    int_fast16_t aExp, shiftCount;
+    uint32_t aSig;
+    uint64_t aSig64, aSigExtra;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1428,7 +1497,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;
@@ -1451,10 +1520,11 @@ int64 float32_to_int64( float32 a STATUS_PARAM )
 int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits32 aSig;
-    bits64 aSig64;
+    int_fast16_t aExp, shiftCount;
+    uint32_t aSig;
+    uint64_t aSig64;
     int64 z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1467,7 +1537,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;
@@ -1476,7 +1546,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;
@@ -1494,14 +1564,15 @@ int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
 float64 float32_to_float64( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits32 aSig;
+    int_fast16_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 commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ));
+        if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloat64( aSign, 0x7FF, 0 );
     }
     if ( aExp == 0 ) {
@@ -1509,12 +1580,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
@@ -1525,14 +1594,15 @@ float64 float32_to_float64( float32 a STATUS_PARAM )
 floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits32 aSig;
+    int_fast16_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 commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( aExp == 0 ) {
@@ -1540,14 +1610,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
@@ -1558,14 +1624,15 @@ floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
 float128 float32_to_float128( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits32 aSig;
+    int_fast16_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 commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloat128( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloat128( aSign, 0x7FFF, 0, 0 );
     }
     if ( aExp == 0 ) {
@@ -1573,12 +1640,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
@@ -1589,10 +1654,11 @@ float128 float32_to_float128( float32 a STATUS_PARAM )
 float32 float32_round_to_int( float32 a STATUS_PARAM)
 {
     flag aSign;
-    int16 aExp;
-    bits32 lastBitMask, roundBitsMask;
+    int_fast16_t aExp;
+    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 ) {
@@ -1602,7 +1668,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) ) {
@@ -1648,9 +1714,9 @@ 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;
-    int16 expDiff;
+    int_fast16_t aExp, bExp, zExp;
+    uint32_t aSig, bSig, zSig;
+    int_fast16_t expDiff;
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1693,7 +1759,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;
@@ -1703,7 +1774,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;
     }
@@ -1722,9 +1793,9 @@ 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;
-    int16 expDiff;
+    int_fast16_t aExp, bExp, zExp;
+    uint32_t aSig, bSig, zSig;
+    int_fast16_t expDiff;
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1796,6 +1867,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 );
@@ -1817,6 +1890,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 );
@@ -1838,10 +1913,13 @@ float32 float32_sub( float32 a, float32 b STATUS_PARAM )
 float32 float32_mul( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
-    int16 aExp, bExp, zExp;
-    bits32 aSig, bSig;
-    bits64 zSig64;
-    bits32 zSig;
+    int_fast16_t aExp, bExp, zExp;
+    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 );
@@ -1879,9 +1957,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;
     }
@@ -1898,8 +1976,10 @@ float32 float32_mul( float32 a, float32 b STATUS_PARAM )
 float32 float32_div( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
-    int16 aExp, bExp, zExp;
-    bits32 aSig, bSig, zSig;
+    int_fast16_t aExp, bExp, zExp;
+    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 );
@@ -1943,9 +2023,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 );
 
@@ -1960,12 +2040,14 @@ float32 float32_div( float32 a, float32 b STATUS_PARAM )
 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;
+    int_fast16_t aExp, bExp, expDiff;
+    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 );
@@ -2007,7 +2089,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;
@@ -2019,8 +2101,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 );
@@ -2039,17 +2121,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;
+    int_fast16_t 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 packFloat32(cSign ^ signflip, cExp, cSig);
+    }
+
+    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
@@ -2059,9 +2348,10 @@ float32 float32_rem( float32 a, float32 b STATUS_PARAM )
 float32 float32_sqrt( float32 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, zExp;
-    bits32 aSig, zSig;
-    bits64 rem, term;
+    int_fast16_t aExp, zExp;
+    uint32_t aSig, zSig;
+    uint64_t rem, term;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -2090,11 +2380,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 );
     }
@@ -2124,30 +2414,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;
+    int_fast16_t aExp;
+    uint32_t aSig;
     float64 r, x, xn;
     int i;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -2191,9 +2482,10 @@ float32 float32_exp2( float32 a STATUS_PARAM )
 float32 float32_log2( float32 a STATUS_PARAM )
 {
     flag aSign, zSign;
-    int16 aExp;
-    bits32 aSig, zSig, i;
+    int_fast16_t aExp;
+    uint32_t aSig, zSig, i;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
@@ -2217,7 +2509,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;
@@ -2232,37 +2524,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 ) )
@@ -2274,21 +2570,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 ) )
@@ -2300,32 +2599,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 );
 }
 
 /*----------------------------------------------------------------------------
@@ -2338,7 +2659,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 ) )
@@ -2352,7 +2675,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 ) );
 
 }
@@ -2367,7 +2690,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 ) )
@@ -2381,11 +2706,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
@@ -2399,8 +2747,9 @@ int float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
 int32 float64_to_int32( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits64 aSig;
+    int_fast16_t aExp, shiftCount;
+    uint64_t aSig;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2426,9 +2775,10 @@ int32 float64_to_int32( float64 a STATUS_PARAM )
 int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits64 aSig, savedASig;
-    int32 z;
+    int_fast16_t aExp, shiftCount;
+    uint64_t aSig, savedASig;
+    int32_t z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2450,7 +2800,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;
@@ -2469,11 +2819,11 @@ int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
 | returned.
 *----------------------------------------------------------------------------*/
 
-int16 float64_to_int16_round_to_zero( float64 a STATUS_PARAM )
+int_fast16_t float64_to_int16_round_to_zero(float64 a STATUS_PARAM)
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits64 aSig, savedASig;
+    int_fast16_t aExp, shiftCount;
+    uint64_t aSig, savedASig;
     int32 z;
 
     aSig = extractFloat64Frac( a );
@@ -2502,7 +2852,7 @@ int16 float64_to_int16_round_to_zero( float64 a STATUS_PARAM )
     if ( ( (int16_t)z < 0 ) ^ aSign ) {
  invalid:
         float_raise( float_flag_invalid STATUS_VAR);
-        return aSign ? (sbits32) 0xffff8000 : 0x7FFF;
+        return aSign ? (int32_t) 0xffff8000 : 0x7FFF;
     }
     if ( ( aSig<<shiftCount ) != savedASig ) {
         STATUS(float_exception_flags) |= float_flag_inexact;
@@ -2523,8 +2873,9 @@ int16 float64_to_int16_round_to_zero( float64 a STATUS_PARAM )
 int64 float64_to_int64( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits64 aSig, aSigExtra;
+    int_fast16_t aExp, shiftCount;
+    uint64_t aSig, aSigExtra;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2540,7 +2891,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;
@@ -2565,9 +2916,10 @@ int64 float64_to_int64( float64 a STATUS_PARAM )
 int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, shiftCount;
-    bits64 aSig;
+    int_fast16_t aExp, shiftCount;
+    uint64_t aSig;
     int64 z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2585,7 +2937,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;
     }
@@ -2595,7 +2947,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;
         }
     }
@@ -2614,15 +2966,16 @@ int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
 float32 float64_to_float32( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 aSig;
-    bits32 zSig;
+    int_fast16_t aExp;
+    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 );
@@ -2646,31 +2999,30 @@ 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, int_fast16_t 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;
+    int_fast16_t aExp;
+    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);
+        return packFloat32(aSign, 0xff, 0);
     }
     if (aExp == 0) {
         int8 shiftCount;
@@ -2686,38 +3038,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;
+    int_fast16_t aExp;
+    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;
         }
-        return packFloat16(aSign, 0x1f, aSig >> 13);
+        /* Infinity */
+        if (!ieee) {
+            float_raise(float_flag_invalid STATUS_VAR);
+            return packFloat16(aSign, 0x1f, 0x3ff);
+        }
+        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;
@@ -2759,7 +3118,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);
         }
     }
@@ -2773,8 +3132,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
@@ -2785,14 +3142,15 @@ bits16 float32_to_float16( float32 a, flag ieee STATUS_PARAM)
 floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 aSig;
+    int_fast16_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 commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) );
+        if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
     }
     if ( aExp == 0 ) {
@@ -2805,10 +3163,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
@@ -2819,14 +3173,15 @@ floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
 float128 float64_to_float128( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 aSig, zSig0, zSig1;
+    int_fast16_t aExp;
+    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 ) {
@@ -2839,8 +3194,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
@@ -2851,10 +3204,11 @@ float128 float64_to_float128( float64 a STATUS_PARAM )
 float64 float64_round_to_int( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 lastBitMask, roundBitsMask;
+    int_fast16_t aExp;
+    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 ) {
@@ -2864,7 +3218,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) ) {
@@ -2923,9 +3277,9 @@ 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;
-    int16 expDiff;
+    int_fast16_t aExp, bExp, zExp;
+    uint64_t aSig, bSig, zSig;
+    int_fast16_t expDiff;
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2968,7 +3322,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;
@@ -2978,7 +3337,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;
     }
@@ -2997,9 +3356,9 @@ 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;
-    int16 expDiff;
+    int_fast16_t aExp, bExp, zExp;
+    uint64_t aSig, bSig, zSig;
+    int_fast16_t expDiff;
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -3071,6 +3430,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 );
@@ -3092,6 +3453,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 );
@@ -3113,8 +3476,11 @@ float64 float64_sub( float64 a, float64 b STATUS_PARAM )
 float64 float64_mul( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign, zSign;
-    int16 aExp, bExp, zExp;
-    bits64 aSig, bSig, zSig0, zSig1;
+    int_fast16_t aExp, bExp, zExp;
+    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 );
@@ -3154,7 +3520,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;
     }
@@ -3171,10 +3537,12 @@ float64 float64_mul( float64 a, float64 b STATUS_PARAM )
 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;
+    int_fast16_t aExp, bExp, zExp;
+    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 );
@@ -3222,7 +3590,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 );
         }
@@ -3241,11 +3609,13 @@ float64 float64_div( float64 a, float64 b STATUS_PARAM )
 float64 float64_rem( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, zSign;
-    int16 aExp, bExp, expDiff;
-    bits64 aSig, bSig;
-    bits64 q, alternateASig;
-    sbits64 sigMean;
+    int_fast16_t aExp, bExp, expDiff;
+    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 );
@@ -3305,17 +3675,243 @@ 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 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;
+    int_fast16_t 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;
+    }
+
+    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 packFloat64(cSign ^ signflip, cExp, cSig);
+    }
+
+    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);
+            if (shiftcount == 0) {
+                zSig0 = (zSig1 >> 1) | (zSig1 & 1);
+                zExp -= 63;
+            } else {
+                shiftcount--;
+                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
@@ -3325,9 +3921,10 @@ float64 float64_rem( float64 a, float64 b STATUS_PARAM )
 float64 float64_sqrt( float64 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp, zExp;
-    bits64 aSig, zSig, doubleZSig;
-    bits64 rem0, rem1, term0, term1;
+    int_fast16_t aExp, zExp;
+    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 );
@@ -3356,7 +3953,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 );
@@ -3375,8 +3972,9 @@ float64 float64_sqrt( float64 a STATUS_PARAM )
 float64 float64_log2( float64 a STATUS_PARAM )
 {
     flag aSign, zSign;
-    int16 aExp;
-    bits64 aSig, aSig0, aSig1, zSig, i;
+    int_fast16_t aExp;
+    uint64_t aSig, aSig0, aSig1, zSig, i;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -3398,7 +3996,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 );
@@ -3415,39 +4013,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 ) )
@@ -3459,22 +4060,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 ) )
        ) {
@@ -3485,31 +4089,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 );
 
 }
 
@@ -3523,7 +4152,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 ) )
@@ -3537,7 +4168,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 ) );
 
 }
@@ -3552,7 +4183,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 ) )
@@ -3566,12 +4199,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-
@@ -3587,12 +4241,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 );
@@ -3614,14 +4268,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 ) {
@@ -3636,7 +4290,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;
@@ -3659,7 +4313,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 );
@@ -3674,7 +4328,7 @@ int64 floatx80_to_int64( floatx80 a STATUS_PARAM )
                ) {
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
             }
-            return (sbits64) LIT64( 0x8000000000000000 );
+            return (int64_t) LIT64( 0x8000000000000000 );
         }
         aSigExtra = 0;
     }
@@ -3699,7 +4353,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 );
@@ -3714,14 +4368,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;
@@ -3740,14 +4394,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 );
     }
@@ -3768,14 +4422,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 );
     }
@@ -3785,8 +4439,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
@@ -3797,22 +4449,20 @@ float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
 float128 floatx80_to_float128( floatx80 a STATUS_PARAM )
 {
     flag aSign;
-    int16 aExp;
-    bits64 aSig, zSig0, zSig1;
+    int_fast16_t aExp;
+    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
@@ -3824,27 +4474,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 ) );
@@ -3897,7 +4547,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 );
@@ -3907,7 +4557,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;
@@ -3916,7 +4566,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;
@@ -3925,7 +4575,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;
@@ -3940,7 +4590,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 );
@@ -3963,7 +4613,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;
 
@@ -3975,7 +4625,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);
@@ -3993,7 +4643,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;
@@ -4005,7 +4655,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;
@@ -4072,7 +4722,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 );
@@ -4083,15 +4733,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);
@@ -4111,7 +4761,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;
     }
@@ -4131,8 +4781,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 );
@@ -4143,15 +4793,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 ) {
@@ -4181,15 +4831,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 );
         }
@@ -4211,8 +4861,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 );
@@ -4221,14 +4871,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 ) {
@@ -4242,7 +4892,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 );
@@ -4307,15 +4957,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;
     }
@@ -4338,7 +4988,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 );
@@ -4350,7 +5000,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;
@@ -4368,31 +5018,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 ) )
            );
 
 }
@@ -4400,8 +5047,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 )
@@ -4409,9 +5057,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;
@@ -4421,7 +5069,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
@@ -4432,9 +5080,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 )
@@ -4442,9 +5090,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;
@@ -4454,7 +5102,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
@@ -4464,28 +5112,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 ) )
            );
 
 }
@@ -4502,9 +5172,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 ) ) {
@@ -4517,7 +5187,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
@@ -4538,9 +5208,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 ) ) {
@@ -4553,7 +5223,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
@@ -4562,9 +5232,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
@@ -4580,7 +5268,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 );
@@ -4609,8 +5297,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 );
@@ -4634,7 +5322,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;
@@ -4657,7 +5345,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 );
@@ -4675,7 +5363,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 );
     }
@@ -4700,7 +5388,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 );
@@ -4722,10 +5410,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;
         }
     }
@@ -4738,7 +5426,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;
         }
     }
@@ -4758,8 +5446,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 );
@@ -4767,7 +5455,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 );
     }
@@ -4793,7 +5481,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 );
@@ -4801,7 +5489,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 );
     }
@@ -4815,8 +5503,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
@@ -4828,7 +5514,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 );
@@ -4836,7 +5522,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 ) );
     }
@@ -4852,8 +5538,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
@@ -4865,7 +5549,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;
 
@@ -4890,9 +5574,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;
                 }
             }
         }
@@ -4906,7 +5590,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) ) {
@@ -4968,7 +5652,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 );
@@ -5017,7 +5701,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;
@@ -5049,7 +5738,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;
 
@@ -5174,7 +5863,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 );
@@ -5238,8 +5927,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 );
@@ -5293,7 +5982,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 );
     }
@@ -5301,7 +5990,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 );
         }
@@ -5322,9 +6011,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 );
@@ -5406,15 +6095,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 );
@@ -5431,8 +6120,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 );
@@ -5464,7 +6153,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 );
@@ -5476,7 +6165,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;
@@ -5492,7 +6181,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.
 *----------------------------------------------------------------------------*/
 
@@ -5504,26 +6194,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 )
@@ -5543,7 +6230,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
@@ -5554,8 +6241,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 )
@@ -5575,7 +6263,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
@@ -5584,14 +6272,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 )
@@ -5599,14 +6307,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 ) )
            );
 
 }
@@ -5638,7 +6349,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
@@ -5674,7 +6385,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
@@ -5683,23 +6394,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) {
@@ -5714,10 +6446,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) {
@@ -5732,10 +6464,10 @@ unsigned int float32_to_uint32_round_to_zero( float32 a STATUS_PARAM )
     return res;
 }
 
-unsigned int float32_to_uint16_round_to_zero( float32 a STATUS_PARAM )
+uint_fast16_t float32_to_uint16_round_to_zero(float32 a STATUS_PARAM)
 {
     int64_t v;
-    unsigned int res;
+    uint_fast16_t res;
 
     v = float32_to_int64_round_to_zero(a STATUS_VAR);
     if (v < 0) {
@@ -5750,10 +6482,10 @@ unsigned int float32_to_uint16_round_to_zero( float32 a STATUS_PARAM )
     return res;
 }
 
-unsigned int float64_to_uint32( float64 a STATUS_PARAM )
+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) {
@@ -5768,10 +6500,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) {
@@ -5786,10 +6518,10 @@ unsigned int float64_to_uint32_round_to_zero( float64 a STATUS_PARAM )
     return res;
 }
 
-unsigned int float64_to_uint16_round_to_zero( float64 a STATUS_PARAM )
+uint_fast16_t float64_to_uint16_round_to_zero(float64 a STATUS_PARAM)
 {
     int64_t v;
-    unsigned int res;
+    uint_fast16_t res;
 
     v = float64_to_int64_round_to_zero(a STATUS_VAR);
     if (v < 0) {
@@ -5832,7 +6564,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 ) ) ||                                 \
@@ -5850,7 +6584,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 {                                                             \
@@ -5878,6 +6612,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 )
 {
@@ -5922,18 +6702,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 )
@@ -5941,6 +6774,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 );
@@ -5949,14 +6788,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 )
@@ -5964,46 +6807,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 )
@@ -6011,9 +6870,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