]> git.proxmox.com Git - qemu.git/blobdiff - fpu/softfloat.c
softfloat: Correctly handle NaNs in float16_to_float32()
[qemu.git] / fpu / softfloat.c
index 5e846200ae7470a6df6b9c5d3872120eba73c924..3abd1708e47b35db3863e9ac50dae66400908c30 100644 (file)
@@ -66,6 +66,33 @@ void set_floatx80_rounding_precision(int val STATUS_PARAM)
 }
 #endif
 
+/*----------------------------------------------------------------------------
+| Returns the fraction bits of the half-precision floating-point value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE uint32_t extractFloat16Frac(float16 a)
+{
+    return float16_val(a) & 0x3ff;
+}
+
+/*----------------------------------------------------------------------------
+| Returns the exponent bits of the half-precision floating-point value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE int16 extractFloat16Exp(float16 a)
+{
+    return (float16_val(a) >> 10) & 0x1f;
+}
+
+/*----------------------------------------------------------------------------
+| Returns the sign bit of the single-precision floating-point value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE flag extractFloat16Sign(float16 a)
+{
+    return float16_val(a)>>15;
+}
+
 /*----------------------------------------------------------------------------
 | Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
 | and 7, and returns the properly rounded 32-bit integer corresponding to the
@@ -175,7 +202,7 @@ static int64 roundAndPackInt64( flag zSign, bits64 absZ0, bits64 absZ1 STATUS_PA
 INLINE bits32 extractFloat32Frac( float32 a )
 {
 
-    return a & 0x007FFFFF;
+    return float32_val(a) & 0x007FFFFF;
 
 }
 
@@ -186,7 +213,7 @@ INLINE bits32 extractFloat32Frac( float32 a )
 INLINE int16 extractFloat32Exp( float32 a )
 {
 
-    return ( a>>23 ) & 0xFF;
+    return ( float32_val(a)>>23 ) & 0xFF;
 
 }
 
@@ -197,8 +224,23 @@ INLINE int16 extractFloat32Exp( float32 a )
 INLINE flag extractFloat32Sign( float32 a )
 {
 
-    return a>>31;
+    return float32_val(a)>>31;
+
+}
 
+/*----------------------------------------------------------------------------
+| 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;
 }
 
 /*----------------------------------------------------------------------------
@@ -233,7 +275,8 @@ static void
 INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
 {
 
-    return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
+    return make_float32(
+          ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig);
 
 }
 
@@ -290,9 +333,10 @@ static float32 roundAndPackFloat32( flag zSign, int16 zExp, bits32 zSig STATUS_P
                   && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
            ) {
             float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
-            return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
+            return packFloat32( zSign, 0xFF, - ( roundIncrement == 0 ));
         }
         if ( zExp < 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloat32( zSign, 0, 0 );
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < -1 )
@@ -337,7 +381,7 @@ static float32
 INLINE bits64 extractFloat64Frac( float64 a )
 {
 
-    return a & LIT64( 0x000FFFFFFFFFFFFF );
+    return float64_val(a) & LIT64( 0x000FFFFFFFFFFFFF );
 
 }
 
@@ -348,7 +392,7 @@ INLINE bits64 extractFloat64Frac( float64 a )
 INLINE int16 extractFloat64Exp( float64 a )
 {
 
-    return ( a>>52 ) & 0x7FF;
+    return ( float64_val(a)>>52 ) & 0x7FF;
 
 }
 
@@ -359,8 +403,23 @@ INLINE int16 extractFloat64Exp( float64 a )
 INLINE flag extractFloat64Sign( float64 a )
 {
 
-    return a>>63;
+    return float64_val(a)>>63;
+
+}
 
+/*----------------------------------------------------------------------------
+| 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;
 }
 
 /*----------------------------------------------------------------------------
@@ -395,7 +454,8 @@ static void
 INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
 {
 
-    return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig;
+    return make_float64(
+        ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig);
 
 }
 
@@ -452,9 +512,10 @@ static float64 roundAndPackFloat64( flag zSign, int16 zExp, bits64 zSig STATUS_P
                   && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
            ) {
             float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
-            return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );
+            return packFloat64( zSign, 0x7FF, - ( roundIncrement == 0 ));
         }
         if ( zExp < 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloat64( zSign, 0, 0 );
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < -1 )
@@ -633,6 +694,7 @@ static floatx80
             goto overflow;
         }
         if ( zExp <= 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloatx80( zSign, 0, 0 );
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < 0 )
@@ -963,6 +1025,7 @@ static float128
             return packFloat128( zSign, 0x7FFF, 0, 0 );
         }
         if ( zExp < 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloat128( zSign, 0, 0, 0 );
             isTiny =
                    ( STATUS(float_detect_tininess) == float_tininess_before_rounding )
                 || ( zExp < -1 )
@@ -1050,7 +1113,7 @@ float32 int32_to_float32( int32 a STATUS_PARAM )
 {
     flag zSign;
 
-    if ( a == 0 ) return 0;
+    if ( a == 0 ) return float32_zero;
     if ( a == (sbits32) 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
     zSign = ( a < 0 );
     return normalizeRoundAndPackFloat32( zSign, 0x9C, zSign ? - a : a STATUS_VAR );
@@ -1070,7 +1133,7 @@ float64 int32_to_float64( int32 a STATUS_PARAM )
     int8 shiftCount;
     bits64 zSig;
 
-    if ( a == 0 ) return 0;
+    if ( a == 0 ) return float64_zero;
     zSign = ( a < 0 );
     absA = zSign ? - a : a;
     shiftCount = countLeadingZeros32( absA ) + 21;
@@ -1144,7 +1207,7 @@ float32 int64_to_float32( int64 a STATUS_PARAM )
     uint64 absA;
     int8 shiftCount;
 
-    if ( a == 0 ) return 0;
+    if ( a == 0 ) return float32_zero;
     zSign = ( a < 0 );
     absA = zSign ? - a : a;
     shiftCount = countLeadingZeros64( absA ) - 40;
@@ -1164,6 +1227,27 @@ float32 int64_to_float32( int64 a STATUS_PARAM )
 
 }
 
+float32 uint64_to_float32( uint64 a STATUS_PARAM )
+{
+    int8 shiftCount;
+
+    if ( a == 0 ) return float32_zero;
+    shiftCount = countLeadingZeros64( a ) - 40;
+    if ( 0 <= shiftCount ) {
+        return packFloat32( 1 > 0, 0x95 - shiftCount, a<<shiftCount );
+    }
+    else {
+        shiftCount += 7;
+        if ( shiftCount < 0 ) {
+            shift64RightJamming( a, - shiftCount, &a );
+        }
+        else {
+            a <<= shiftCount;
+        }
+        return roundAndPackFloat32( 1 > 0, 0x9C - shiftCount, a STATUS_VAR );
+    }
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of converting the 64-bit two's complement integer `a'
 | to the double-precision floating-point format.  The conversion is performed
@@ -1174,7 +1258,7 @@ float64 int64_to_float64( int64 a STATUS_PARAM )
 {
     flag zSign;
 
-    if ( a == 0 ) return 0;
+    if ( a == 0 ) return float64_zero;
     if ( a == (sbits64) LIT64( 0x8000000000000000 ) ) {
         return packFloat64( 1, 0x43E, 0 );
     }
@@ -1183,6 +1267,13 @@ float64 int64_to_float64( int64 a STATUS_PARAM )
 
 }
 
+float64 uint64_to_float64( uint64 a STATUS_PARAM )
+{
+    if ( a == 0 ) return float64_zero;
+    return normalizeRoundAndPackFloat64( 0, 0x43C, a STATUS_VAR );
+
+}
+
 #ifdef FLOATX80
 
 /*----------------------------------------------------------------------------
@@ -1262,6 +1353,7 @@ int32 float32_to_int32( float32 a STATUS_PARAM )
     bits32 aSig;
     bits64 aSig64;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
@@ -1291,13 +1383,14 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
     int16 aExp, shiftCount;
     bits32 aSig;
     int32 z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     shiftCount = aExp - 0x9E;
     if ( 0 <= shiftCount ) {
-        if ( a != 0xCF000000 ) {
+        if ( float32_val(a) != 0xCF000000 ) {
             float_raise( float_flag_invalid STATUS_VAR);
             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
         }
@@ -1317,6 +1410,55 @@ int32 float32_to_int32_round_to_zero( float32 a STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the result of converting the single-precision floating-point value
+| `a' to the 16-bit two's complement integer format.  The conversion is
+| performed according to the IEC/IEEE Standard for Binary Floating-Point
+| Arithmetic, except that the conversion is always rounded toward zero.
+| If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
+| the conversion overflows, the largest integer with the same sign as `a' is
+| returned.
+*----------------------------------------------------------------------------*/
+
+int16 float32_to_int16_round_to_zero( float32 a STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits32 aSig;
+    int32 z;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    shiftCount = aExp - 0x8E;
+    if ( 0 <= shiftCount ) {
+        if ( float32_val(a) != 0xC7000000 ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+            if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
+                return 0x7FFF;
+            }
+        }
+        return (sbits32) 0xffff8000;
+    }
+    else if ( aExp <= 0x7E ) {
+        if ( aExp | aSig ) {
+            STATUS(float_exception_flags) |= float_flag_inexact;
+        }
+        return 0;
+    }
+    shiftCount -= 0x10;
+    aSig = ( aSig | 0x00800000 )<<8;
+    z = aSig>>( - shiftCount );
+    if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
+        STATUS(float_exception_flags) |= float_flag_inexact;
+    }
+    if ( aSign ) {
+        z = - z;
+    }
+    return z;
+
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of converting the single-precision floating-point value
 | `a' to the 64-bit two's complement integer format.  The conversion is
@@ -1333,6 +1475,7 @@ int64 float32_to_int64( float32 a STATUS_PARAM )
     int16 aExp, shiftCount;
     bits32 aSig;
     bits64 aSig64, aSigExtra;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
@@ -1370,13 +1513,14 @@ int64 float32_to_int64_round_to_zero( float32 a STATUS_PARAM )
     bits32 aSig;
     bits64 aSig64;
     int64 z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     shiftCount = aExp - 0xBE;
     if ( 0 <= shiftCount ) {
-        if ( a != 0xDF000000 ) {
+        if ( float32_val(a) != 0xDF000000 ) {
             float_raise( float_flag_invalid STATUS_VAR);
             if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) {
                 return LIT64( 0x7FFFFFFFFFFFFFFF );
@@ -1411,12 +1555,13 @@ float64 float32_to_float64( float32 a STATUS_PARAM )
     flag aSign;
     int16 aExp;
     bits32 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 ) {
@@ -1443,11 +1588,12 @@ floatx80 float32_to_floatx80( float32 a STATUS_PARAM )
     int16 aExp;
     bits32 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 ) {
@@ -1476,11 +1622,12 @@ float128 float32_to_float128( float32 a STATUS_PARAM )
     int16 aExp;
     bits32 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 ) {
@@ -1507,7 +1654,8 @@ float32 float32_round_to_int( float32 a STATUS_PARAM)
     int16 aExp;
     bits32 lastBitMask, roundBitsMask;
     int8 roundingMode;
-    float32 z;
+    bits32 z;
+    a = float32_squash_input_denormal(a STATUS_VAR);
 
     aExp = extractFloat32Exp( a );
     if ( 0x96 <= aExp ) {
@@ -1517,7 +1665,7 @@ float32 float32_round_to_int( float32 a STATUS_PARAM)
         return a;
     }
     if ( aExp <= 0x7E ) {
-        if ( (bits32) ( a<<1 ) == 0 ) return a;
+        if ( (bits32) ( float32_val(a)<<1 ) == 0 ) return a;
         STATUS(float_exception_flags) |= float_flag_inexact;
         aSign = extractFloat32Sign( a );
         switch ( STATUS(float_rounding_mode) ) {
@@ -1527,29 +1675,29 @@ float32 float32_round_to_int( float32 a STATUS_PARAM)
             }
             break;
          case float_round_down:
-            return aSign ? 0xBF800000 : 0;
+            return make_float32(aSign ? 0xBF800000 : 0);
          case float_round_up:
-            return aSign ? 0x80000000 : 0x3F800000;
+            return make_float32(aSign ? 0x80000000 : 0x3F800000);
         }
         return packFloat32( aSign, 0, 0 );
     }
     lastBitMask = 1;
     lastBitMask <<= 0x96 - aExp;
     roundBitsMask = lastBitMask - 1;
-    z = a;
+    z = float32_val(a);
     roundingMode = STATUS(float_rounding_mode);
     if ( roundingMode == float_round_nearest_even ) {
         z += lastBitMask>>1;
         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
     }
     else if ( roundingMode != float_round_to_zero ) {
-        if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
+        if ( extractFloat32Sign( make_float32(z) ) ^ ( roundingMode == float_round_up ) ) {
             z += roundBitsMask;
         }
     }
     z &= ~ roundBitsMask;
-    if ( z != a ) STATUS(float_exception_flags) |= float_flag_inexact;
-    return z;
+    if ( z != float32_val(a) ) STATUS(float_exception_flags) |= float_flag_inexact;
+    return make_float32(z);
 
 }
 
@@ -1607,7 +1755,10 @@ static float32 addFloat32Sigs( float32 a, float32 b, flag zSign STATUS_PARAM)
             if ( aSig | bSig ) return propagateFloat32NaN( a, b STATUS_VAR );
             return a;
         }
-        if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
+        if ( aExp == 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloat32( zSign, 0, 0 );
+            return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
+        }
         zSig = 0x40000000 + aSig + bSig;
         zExp = aExp;
         goto roundAndPack;
@@ -1708,6 +1859,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 );
@@ -1729,6 +1882,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 );
@@ -1755,6 +1910,9 @@ float32 float32_mul( float32 a, float32 b STATUS_PARAM )
     bits64 zSig64;
     bits32 zSig;
 
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
+
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
@@ -1812,6 +1970,8 @@ float32 float32_div( float32 a, float32 b STATUS_PARAM )
     flag aSign, bSign, zSign;
     int16 aExp, bExp, zExp;
     bits32 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 );
@@ -1871,20 +2031,21 @@ float32 float32_div( float32 a, float32 b STATUS_PARAM )
 
 float32 float32_rem( float32 a, float32 b STATUS_PARAM )
 {
-    flag aSign, bSign, zSign;
+    flag aSign, zSign;
     int16 aExp, bExp, expDiff;
     bits32 aSig, bSig;
     bits32 q;
     bits64 aSig64, bSig64, q64;
     bits32 alternateASig;
     sbits32 sigMean;
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    b = float32_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat32Frac( a );
     aExp = extractFloat32Exp( a );
     aSign = extractFloat32Sign( a );
     bSig = extractFloat32Frac( b );
     bExp = extractFloat32Exp( b );
-    bSign = extractFloat32Sign( b );
     if ( aExp == 0xFF ) {
         if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
             return propagateFloat32NaN( a, b STATUS_VAR );
@@ -1975,12 +2136,13 @@ float32 float32_sqrt( float32 a STATUS_PARAM )
     int16 aExp, zExp;
     bits32 aSig, zSig;
     bits64 rem, term;
+    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, 0 STATUS_VAR );
+        if ( aSig ) return propagateFloat32NaN( a, float32_zero STATUS_VAR );
         if ( ! aSign ) return a;
         float_raise( float_flag_invalid STATUS_VAR);
         return float32_default_nan;
@@ -1991,7 +2153,7 @@ float32 float32_sqrt( float32 a STATUS_PARAM )
         return float32_default_nan;
     }
     if ( aExp == 0 ) {
-        if ( aSig == 0 ) return 0;
+        if ( aSig == 0 ) return float32_zero;
         normalizeFloat32Subnormal( aSig, &aExp, &aSig );
     }
     zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
@@ -2017,14 +2179,144 @@ float32 float32_sqrt( float32 a STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the binary exponential of the single-precision floating-point value
+| `a'. The operation is performed according to the IEC/IEEE Standard for
+| Binary Floating-Point Arithmetic.
+|
+| Uses the following identities:
+|
+| 1. -------------------------------------------------------------------------
+|      x    x*ln(2)
+|     2  = e
+|
+| 2. -------------------------------------------------------------------------
+|                      2     3     4     5           n
+|      x        x     x     x     x     x           x
+|     e  = 1 + --- + --- + --- + --- + --- + ... + --- + ...
+|               1!    2!    3!    4!    5!          n!
+*----------------------------------------------------------------------------*/
+
+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 */
+};
+
+float32 float32_exp2( float32 a STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp;
+    bits32 aSig;
+    float64 r, x, xn;
+    int i;
+    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, float32_zero STATUS_VAR );
+        return (aSign) ? float32_zero : a;
+    }
+    if (aExp == 0) {
+        if (aSig == 0) return float32_one;
+    }
+
+    float_raise( float_flag_inexact STATUS_VAR);
+
+    /* ******************************* */
+    /* using float64 for approximation */
+    /* ******************************* */
+    x = float32_to_float64(a STATUS_VAR);
+    x = float64_mul(x, float64_ln2 STATUS_VAR);
+
+    xn = x;
+    r = float64_one;
+    for (i = 0 ; i < 15 ; i++) {
+        float64 f;
+
+        f = float64_mul(xn, float32_exp2_coefficients[i] STATUS_VAR);
+        r = float64_add(r, f STATUS_VAR);
+
+        xn = float64_mul(xn, x STATUS_VAR);
+    }
+
+    return float64_to_float32(r, status);
+}
+
+/*----------------------------------------------------------------------------
+| Returns the binary log of the single-precision floating-point value `a'.
+| The operation is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+float32 float32_log2( float32 a STATUS_PARAM )
+{
+    flag aSign, zSign;
+    int16 aExp;
+    bits32 aSig, zSig, i;
+
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat32( 1, 0xFF, 0 );
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+    }
+    if ( aSign ) {
+        float_raise( float_flag_invalid STATUS_VAR);
+        return float32_default_nan;
+    }
+    if ( aExp == 0xFF ) {
+        if ( aSig ) return propagateFloat32NaN( a, float32_zero STATUS_VAR );
+        return a;
+    }
+
+    aExp -= 0x7F;
+    aSig |= 0x00800000;
+    zSign = aExp < 0;
+    zSig = aExp << 23;
+
+    for (i = 1 << 22; i > 0; i >>= 1) {
+        aSig = ( (bits64)aSig * aSig ) >> 23;
+        if ( aSig & 0x01000000 ) {
+            aSig >>= 1;
+            zSig |= i;
+        }
+    }
+
+    if ( zSign )
+        zSig = -zSig;
+
+    return normalizeRoundAndPackFloat32( zSign, 0x85, zSig STATUS_VAR );
+}
+
 /*----------------------------------------------------------------------------
 | Returns 1 if the single-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.
 *----------------------------------------------------------------------------*/
 
-flag float32_eq( float32 a, float32 b STATUS_PARAM )
+int float32_eq( 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 ) )
@@ -2034,7 +2326,8 @@ flag float32_eq( float32 a, float32 b STATUS_PARAM )
         }
         return 0;
     }
-    return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
+    return ( float32_val(a) == float32_val(b) ) ||
+            ( (bits32) ( ( float32_val(a) | float32_val(b) )<<1 ) == 0 );
 
 }
 
@@ -2045,9 +2338,12 @@ flag float32_eq( float32 a, float32 b STATUS_PARAM )
 | Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float32_le( float32 a, float32 b STATUS_PARAM )
+int float32_le( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -2057,8 +2353,10 @@ flag float32_le( float32 a, float32 b STATUS_PARAM )
     }
     aSign = extractFloat32Sign( a );
     bSign = extractFloat32Sign( b );
-    if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
-    return ( a == b ) || ( aSign ^ ( a < b ) );
+    av = float32_val(a);
+    bv = float32_val(b);
+    if ( aSign != bSign ) return aSign || ( (bits32) ( ( av | bv )<<1 ) == 0 );
+    return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
 
@@ -2068,9 +2366,12 @@ flag float32_le( float32 a, float32 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float32_lt( float32 a, float32 b STATUS_PARAM )
+int float32_lt( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -2080,8 +2381,10 @@ flag float32_lt( float32 a, float32 b STATUS_PARAM )
     }
     aSign = extractFloat32Sign( a );
     bSign = extractFloat32Sign( b );
-    if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
-    return ( a != b ) && ( aSign ^ ( a < b ) );
+    av = float32_val(a);
+    bv = float32_val(b);
+    if ( aSign != bSign ) return aSign && ( (bits32) ( ( av | bv )<<1 ) != 0 );
+    return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
@@ -2092,8 +2395,11 @@ flag float32_lt( float32 a, float32 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
+int float32_eq_signaling( 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 ) )
@@ -2101,7 +2407,9 @@ flag float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
         float_raise( float_flag_invalid STATUS_VAR);
         return 0;
     }
-    return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
+    av = float32_val(a);
+    bv = float32_val(b);
+    return ( av == bv ) || ( (bits32) ( ( av | bv )<<1 ) == 0 );
 
 }
 
@@ -2112,9 +2420,12 @@ flag float32_eq_signaling( float32 a, float32 b STATUS_PARAM )
 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float32_le_quiet( float32 a, float32 b STATUS_PARAM )
+int float32_le_quiet( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -2126,8 +2437,10 @@ flag float32_le_quiet( float32 a, float32 b STATUS_PARAM )
     }
     aSign = extractFloat32Sign( a );
     bSign = extractFloat32Sign( b );
-    if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
-    return ( a == b ) || ( aSign ^ ( a < b ) );
+    av = float32_val(a);
+    bv = float32_val(b);
+    if ( aSign != bSign ) return aSign || ( (bits32) ( ( av | bv )<<1 ) == 0 );
+    return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
 
@@ -2138,9 +2451,12 @@ flag float32_le_quiet( float32 a, float32 b STATUS_PARAM )
 | Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
+int float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -2152,8 +2468,10 @@ flag float32_lt_quiet( float32 a, float32 b STATUS_PARAM )
     }
     aSign = extractFloat32Sign( a );
     bSign = extractFloat32Sign( b );
-    if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
-    return ( a != b ) && ( aSign ^ ( a < b ) );
+    av = float32_val(a);
+    bv = float32_val(b);
+    if ( aSign != bSign ) return aSign && ( (bits32) ( ( av | bv )<<1 ) != 0 );
+    return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
@@ -2172,6 +2490,7 @@ int32 float64_to_int32( float64 a STATUS_PARAM )
     flag aSign;
     int16 aExp, shiftCount;
     bits64 aSig;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2200,6 +2519,7 @@ int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
     int16 aExp, shiftCount;
     bits64 aSig, savedASig;
     int32 z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2230,6 +2550,57 @@ int32 float64_to_int32_round_to_zero( float64 a STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the result of converting the double-precision floating-point value
+| `a' to the 16-bit two's complement integer format.  The conversion is
+| performed according to the IEC/IEEE Standard for Binary Floating-Point
+| Arithmetic, except that the conversion is always rounded toward zero.
+| If `a' is a NaN, the largest positive integer is returned.  Otherwise, if
+| the conversion overflows, the largest integer with the same sign as `a' is
+| returned.
+*----------------------------------------------------------------------------*/
+
+int16 float64_to_int16_round_to_zero( float64 a STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits64 aSig, savedASig;
+    int32 z;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    if ( 0x40E < aExp ) {
+        if ( ( aExp == 0x7FF ) && aSig ) {
+            aSign = 0;
+        }
+        goto invalid;
+    }
+    else if ( aExp < 0x3FF ) {
+        if ( aExp || aSig ) {
+            STATUS(float_exception_flags) |= float_flag_inexact;
+        }
+        return 0;
+    }
+    aSig |= LIT64( 0x0010000000000000 );
+    shiftCount = 0x433 - aExp;
+    savedASig = aSig;
+    aSig >>= shiftCount;
+    z = aSig;
+    if ( aSign ) {
+        z = - z;
+    }
+    if ( ( (int16_t)z < 0 ) ^ aSign ) {
+ invalid:
+        float_raise( float_flag_invalid STATUS_VAR);
+        return aSign ? (sbits32) 0xffff8000 : 0x7FFF;
+    }
+    if ( ( aSig<<shiftCount ) != savedASig ) {
+        STATUS(float_exception_flags) |= float_flag_inexact;
+    }
+    return z;
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of converting the double-precision floating-point value
 | `a' to the 64-bit two's complement integer format.  The conversion is
@@ -2245,6 +2616,7 @@ int64 float64_to_int64( float64 a STATUS_PARAM )
     flag aSign;
     int16 aExp, shiftCount;
     bits64 aSig, aSigExtra;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2288,6 +2660,7 @@ int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
     int16 aExp, shiftCount;
     bits64 aSig;
     int64 z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2296,7 +2669,7 @@ int64 float64_to_int64_round_to_zero( float64 a STATUS_PARAM )
     shiftCount = aExp - 0x433;
     if ( 0 <= shiftCount ) {
         if ( 0x43E <= aExp ) {
-            if ( a != LIT64( 0xC3E0000000000000 ) ) {
+            if ( float64_val(a) != LIT64( 0xC3E0000000000000 ) ) {
                 float_raise( float_flag_invalid STATUS_VAR);
                 if (    ! aSign
                      || (    ( aExp == 0x7FF )
@@ -2337,12 +2710,13 @@ float32 float64_to_float32( float64 a STATUS_PARAM )
     int16 aExp;
     bits64 aSig;
     bits32 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 );
@@ -2355,6 +2729,150 @@ float32 float64_to_float32( float64 a STATUS_PARAM )
 
 }
 
+
+/*----------------------------------------------------------------------------
+| Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
+| half-precision floating-point value, returning the result.  After being
+| shifted into the proper positions, the three fields are simply added
+| together to form the result.  This means that any integer portion of `zSig'
+| will be added into the exponent.  Since a properly normalized significand
+| will have an integer portion equal to 1, the `zExp' input should be 1 less
+| than the desired result exponent whenever `zSig' is a complete, normalized
+| significand.
+*----------------------------------------------------------------------------*/
+static float16 packFloat16(flag zSign, int16 zExp, bits16 zSig)
+{
+    return make_float16(
+        (((bits32)zSign) << 15) + (((bits32)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(float16 a, flag ieee STATUS_PARAM)
+{
+    flag aSign;
+    int16 aExp;
+    bits32 aSig;
+
+    aSign = extractFloat16Sign(a);
+    aExp = extractFloat16Exp(a);
+    aSig = extractFloat16Frac(a);
+
+    if (aExp == 0x1f && ieee) {
+        if (aSig) {
+            return commonNaNToFloat32(float16ToCommonNaN(a STATUS_VAR) STATUS_VAR);
+        }
+        return packFloat32(aSign, 0xff, aSig << 13);
+    }
+    if (aExp == 0) {
+        int8 shiftCount;
+
+        if (aSig == 0) {
+            return packFloat32(aSign, 0, 0);
+        }
+
+        shiftCount = countLeadingZeros32( aSig ) - 21;
+        aSig = aSig << shiftCount;
+        aExp = -shiftCount;
+    }
+    return packFloat32( aSign, aExp + 0x70, aSig << 13);
+}
+
+float16 float32_to_float16(float32 a, flag ieee STATUS_PARAM)
+{
+    flag aSign;
+    int16 aExp;
+    bits32 aSig;
+    bits32 mask;
+    bits32 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) {
+            /* Input is a NaN */
+            float16 r = commonNaNToFloat16( float32ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
+            if (!ieee) {
+                return packFloat16(aSign, 0, 0);
+            }
+            return r;
+        }
+        /* Infinity */
+        if (!ieee) {
+            float_raise(float_flag_invalid STATUS_VAR);
+            return packFloat16(aSign, 0x1f, 0x3ff);
+        }
+        return packFloat16(aSign, 0x1f, 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 = 0x00ffffff;
+        if (aExp >= -24) {
+            mask >>= 25 + aExp;
+        }
+    } else {
+        mask = 0x00001fff;
+    }
+    if (aSig & mask) {
+        float_raise( float_flag_underflow STATUS_VAR );
+        roundingMode = STATUS(float_rounding_mode);
+        switch (roundingMode) {
+        case float_round_nearest_even:
+            increment = (mask + 1) >> 1;
+            if ((aSig & mask) == increment) {
+                increment = aSig & (increment << 1);
+            }
+            break;
+        case float_round_up:
+            increment = aSign ? 0 : mask;
+            break;
+        case float_round_down:
+            increment = aSign ? mask : 0;
+            break;
+        default: /* round_to_zero */
+            increment = 0;
+            break;
+        }
+        aSig += increment;
+        if (aSig >= 0x01000000) {
+            aSig >>= 1;
+            aExp++;
+        }
+    } else if (aExp < -14
+          && STATUS(float_detect_tininess) == float_tininess_before_rounding) {
+        float_raise( float_flag_underflow STATUS_VAR);
+    }
+
+    if (ieee) {
+        if (aExp > 15) {
+            float_raise( float_flag_overflow | float_flag_inexact STATUS_VAR);
+            return packFloat16(aSign, 0x1f, 0);
+        }
+    } else {
+        if (aExp > 16) {
+            float_raise(float_flag_invalid | float_flag_inexact STATUS_VAR);
+            return packFloat16(aSign, 0x1f, 0x3ff);
+        }
+    }
+    if (aExp < -24) {
+        return packFloat16(aSign, 0, 0);
+    }
+    if (aExp < -14) {
+        aSig >>= -14 - aExp;
+        aExp = -14;
+    }
+    return packFloat16(aSign, aExp + 14, aSig >> 13);
+}
+
 #ifdef FLOATX80
 
 /*----------------------------------------------------------------------------
@@ -2370,11 +2888,12 @@ floatx80 float64_to_floatx80( float64 a STATUS_PARAM )
     int16 aExp;
     bits64 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 ) {
@@ -2404,11 +2923,12 @@ float128 float64_to_float128( float64 a STATUS_PARAM )
     int16 aExp;
     bits64 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 ) {
@@ -2436,7 +2956,8 @@ float64 float64_round_to_int( float64 a STATUS_PARAM )
     int16 aExp;
     bits64 lastBitMask, roundBitsMask;
     int8 roundingMode;
-    float64 z;
+    bits64 z;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aExp = extractFloat64Exp( a );
     if ( 0x433 <= aExp ) {
@@ -2446,7 +2967,7 @@ float64 float64_round_to_int( float64 a STATUS_PARAM )
         return a;
     }
     if ( aExp < 0x3FF ) {
-        if ( (bits64) ( a<<1 ) == 0 ) return a;
+        if ( (bits64) ( float64_val(a)<<1 ) == 0 ) return a;
         STATUS(float_exception_flags) |= float_flag_inexact;
         aSign = extractFloat64Sign( a );
         switch ( STATUS(float_rounding_mode) ) {
@@ -2456,33 +2977,45 @@ float64 float64_round_to_int( float64 a STATUS_PARAM )
             }
             break;
          case float_round_down:
-            return aSign ? LIT64( 0xBFF0000000000000 ) : 0;
+            return make_float64(aSign ? LIT64( 0xBFF0000000000000 ) : 0);
          case float_round_up:
-            return
-            aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );
+            return make_float64(
+            aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 ));
         }
         return packFloat64( aSign, 0, 0 );
     }
     lastBitMask = 1;
     lastBitMask <<= 0x433 - aExp;
     roundBitsMask = lastBitMask - 1;
-    z = a;
+    z = float64_val(a);
     roundingMode = STATUS(float_rounding_mode);
     if ( roundingMode == float_round_nearest_even ) {
         z += lastBitMask>>1;
         if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
     }
     else if ( roundingMode != float_round_to_zero ) {
-        if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {
+        if ( extractFloat64Sign( make_float64(z) ) ^ ( roundingMode == float_round_up ) ) {
             z += roundBitsMask;
         }
     }
     z &= ~ roundBitsMask;
-    if ( z != a ) STATUS(float_exception_flags) |= float_flag_inexact;
-    return z;
+    if ( z != float64_val(a) )
+        STATUS(float_exception_flags) |= float_flag_inexact;
+    return make_float64(z);
 
 }
 
+float64 float64_trunc_to_int( float64 a STATUS_PARAM)
+{
+    int oldmode;
+    float64 res;
+    oldmode = STATUS(float_rounding_mode);
+    STATUS(float_rounding_mode) = float_round_to_zero;
+    res = float64_round_to_int(a STATUS_VAR);
+    STATUS(float_rounding_mode) = oldmode;
+    return res;
+}
+
 /*----------------------------------------------------------------------------
 | Returns the result of adding the absolute values of the double-precision
 | floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
@@ -2537,7 +3070,10 @@ static float64 addFloat64Sigs( float64 a, float64 b, flag zSign STATUS_PARAM )
             if ( aSig | bSig ) return propagateFloat64NaN( a, b STATUS_VAR );
             return a;
         }
-        if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
+        if ( aExp == 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloat64( zSign, 0, 0 );
+            return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
+        }
         zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
         zExp = aExp;
         goto roundAndPack;
@@ -2638,6 +3174,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 );
@@ -2659,6 +3197,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 );
@@ -2683,6 +3223,9 @@ float64 float64_mul( float64 a, float64 b STATUS_PARAM )
     int16 aExp, bExp, zExp;
     bits64 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 );
     aSign = extractFloat64Sign( a );
@@ -2742,6 +3285,8 @@ float64 float64_div( float64 a, float64 b STATUS_PARAM )
     bits64 aSig, bSig, zSig;
     bits64 rem0, rem1;
     bits64 term0, term1;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    b = float64_squash_input_denormal(b STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2807,18 +3352,19 @@ float64 float64_div( float64 a, float64 b STATUS_PARAM )
 
 float64 float64_rem( float64 a, float64 b STATUS_PARAM )
 {
-    flag aSign, bSign, zSign;
+    flag aSign, zSign;
     int16 aExp, bExp, expDiff;
     bits64 aSig, bSig;
     bits64 q, alternateASig;
     sbits64 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 );
     bSig = extractFloat64Frac( b );
     bExp = extractFloat64Exp( b );
-    bSign = extractFloat64Sign( b );
     if ( aExp == 0x7FF ) {
         if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
             return propagateFloat64NaN( a, b STATUS_VAR );
@@ -2896,6 +3442,7 @@ float64 float64_sqrt( float64 a STATUS_PARAM )
     int16 aExp, zExp;
     bits64 aSig, zSig, doubleZSig;
     bits64 rem0, rem1, term0, term1;
+    a = float64_squash_input_denormal(a STATUS_VAR);
 
     aSig = extractFloat64Frac( a );
     aExp = extractFloat64Exp( a );
@@ -2912,7 +3459,7 @@ float64 float64_sqrt( float64 a STATUS_PARAM )
         return float64_default_nan;
     }
     if ( aExp == 0 ) {
-        if ( aSig == 0 ) return 0;
+        if ( aSig == 0 ) return float64_zero;
         normalizeFloat64Subnormal( aSig, &aExp, &aSig );
     }
     zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
@@ -2935,14 +3482,64 @@ float64 float64_sqrt( float64 a STATUS_PARAM )
 
 }
 
+/*----------------------------------------------------------------------------
+| Returns the binary log of the double-precision floating-point value `a'.
+| The operation is performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+float64 float64_log2( float64 a STATUS_PARAM )
+{
+    flag aSign, zSign;
+    int16 aExp;
+    bits64 aSig, aSig0, aSig1, zSig, i;
+    a = float64_squash_input_denormal(a STATUS_VAR);
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat64( 1, 0x7FF, 0 );
+        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
+    }
+    if ( aSign ) {
+        float_raise( float_flag_invalid STATUS_VAR);
+        return float64_default_nan;
+    }
+    if ( aExp == 0x7FF ) {
+        if ( aSig ) return propagateFloat64NaN( a, float64_zero STATUS_VAR );
+        return a;
+    }
+
+    aExp -= 0x3FF;
+    aSig |= LIT64( 0x0010000000000000 );
+    zSign = aExp < 0;
+    zSig = (bits64)aExp << 52;
+    for (i = 1LL << 51; i > 0; i >>= 1) {
+        mul64To128( aSig, aSig, &aSig0, &aSig1 );
+        aSig = ( aSig0 << 12 ) | ( aSig1 >> 52 );
+        if ( aSig & LIT64( 0x0020000000000000 ) ) {
+            aSig >>= 1;
+            zSig |= i;
+        }
+    }
+
+    if ( zSign )
+        zSig = -zSig;
+    return normalizeRoundAndPackFloat64( zSign, 0x408, zSig STATUS_VAR );
+}
+
 /*----------------------------------------------------------------------------
 | Returns 1 if the 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.
 *----------------------------------------------------------------------------*/
 
-flag float64_eq( float64 a, float64 b STATUS_PARAM )
+int float64_eq( 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 ) )
@@ -2952,7 +3549,9 @@ flag float64_eq( float64 a, float64 b STATUS_PARAM )
         }
         return 0;
     }
-    return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
+    av = float64_val(a);
+    bv = float64_val(b);
+    return ( av == bv ) || ( (bits64) ( ( av | bv )<<1 ) == 0 );
 
 }
 
@@ -2963,9 +3562,12 @@ flag float64_eq( float64 a, float64 b STATUS_PARAM )
 | Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float64_le( float64 a, float64 b STATUS_PARAM )
+int float64_le( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -2975,8 +3577,10 @@ flag float64_le( float64 a, float64 b STATUS_PARAM )
     }
     aSign = extractFloat64Sign( a );
     bSign = extractFloat64Sign( b );
-    if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
-    return ( a == b ) || ( aSign ^ ( a < b ) );
+    av = float64_val(a);
+    bv = float64_val(b);
+    if ( aSign != bSign ) return aSign || ( (bits64) ( ( av | bv )<<1 ) == 0 );
+    return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
 
@@ -2986,10 +3590,13 @@ flag float64_le( float64 a, float64 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float64_lt( float64 a, float64 b STATUS_PARAM )
+int float64_lt( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
        ) {
@@ -2998,8 +3605,10 @@ flag float64_lt( float64 a, float64 b STATUS_PARAM )
     }
     aSign = extractFloat64Sign( a );
     bSign = extractFloat64Sign( b );
-    if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
-    return ( a != b ) && ( aSign ^ ( a < b ) );
+    av = float64_val(a);
+    bv = float64_val(b);
+    if ( aSign != bSign ) return aSign && ( (bits64) ( ( av | bv )<<1 ) != 0 );
+    return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
@@ -3010,8 +3619,11 @@ flag float64_lt( float64 a, float64 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
+int float64_eq_signaling( 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 ) )
@@ -3019,7 +3631,9 @@ flag float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
         float_raise( float_flag_invalid STATUS_VAR);
         return 0;
     }
-    return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
+    av = float64_val(a);
+    bv = float64_val(b);
+    return ( av == bv ) || ( (bits64) ( ( av | bv )<<1 ) == 0 );
 
 }
 
@@ -3030,9 +3644,12 @@ flag float64_eq_signaling( float64 a, float64 b STATUS_PARAM )
 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float64_le_quiet( float64 a, float64 b STATUS_PARAM )
+int float64_le_quiet( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -3044,8 +3661,10 @@ flag float64_le_quiet( float64 a, float64 b STATUS_PARAM )
     }
     aSign = extractFloat64Sign( a );
     bSign = extractFloat64Sign( b );
-    if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
-    return ( a == b ) || ( aSign ^ ( a < b ) );
+    av = float64_val(a);
+    bv = float64_val(b);
+    if ( aSign != bSign ) return aSign || ( (bits64) ( ( av | bv )<<1 ) == 0 );
+    return ( av == bv ) || ( aSign ^ ( av < bv ) );
 
 }
 
@@ -3056,9 +3675,12 @@ flag float64_le_quiet( float64 a, float64 b STATUS_PARAM )
 | Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
+int float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
 {
     flag aSign, bSign;
+    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 ) )
@@ -3070,8 +3692,10 @@ flag float64_lt_quiet( float64 a, float64 b STATUS_PARAM )
     }
     aSign = extractFloat64Sign( a );
     bSign = extractFloat64Sign( b );
-    if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
-    return ( a != b ) && ( aSign ^ ( a < b ) );
+    av = float64_val(a);
+    bv = float64_val(b);
+    if ( aSign != bSign ) return aSign && ( (bits64) ( ( av | bv )<<1 ) != 0 );
+    return ( av != bv ) && ( aSign ^ ( av < bv ) );
 
 }
 
@@ -3251,7 +3875,7 @@ float32 floatx80_to_float32( floatx80 a STATUS_PARAM )
     aSign = extractFloatx80Sign( a );
     if ( aExp == 0x7FFF ) {
         if ( (bits64) ( aSig<<1 ) ) {
-            return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) );
+            return commonNaNToFloat32( floatx80ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloat32( aSign, 0xFF, 0 );
     }
@@ -3279,7 +3903,7 @@ float64 floatx80_to_float64( floatx80 a STATUS_PARAM )
     aSign = extractFloatx80Sign( a );
     if ( aExp == 0x7FFF ) {
         if ( (bits64) ( aSig<<1 ) ) {
-            return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) );
+            return commonNaNToFloat64( floatx80ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
         }
         return packFloat64( aSign, 0x7FF, 0 );
     }
@@ -3308,7 +3932,7 @@ float128 floatx80_to_float128( floatx80 a STATUS_PARAM )
     aExp = extractFloatx80Exp( a );
     aSign = extractFloatx80Sign( a );
     if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) {
-        return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) );
+        return commonNaNToFloat128( floatx80ToCommonNaN( a STATUS_VAR ) STATUS_VAR );
     }
     shift128Right( aSig<<1, 0, 16, &zSig0, &zSig1 );
     return packFloat128( aSign, aExp, zSig0, zSig1 );
@@ -3713,7 +4337,7 @@ floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM )
 
 floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
 {
-    flag aSign, bSign, zSign;
+    flag aSign, zSign;
     int32 aExp, bExp, expDiff;
     bits64 aSig0, aSig1, bSig;
     bits64 q, term0, term1, alternateASig0, alternateASig1;
@@ -3724,7 +4348,6 @@ floatx80 floatx80_rem( floatx80 a, floatx80 b STATUS_PARAM )
     aSign = extractFloatx80Sign( a );
     bSig = extractFloatx80Frac( b );
     bExp = extractFloatx80Exp( b );
-    bSign = extractFloatx80Sign( b );
     if ( aExp == 0x7FFF ) {
         if (    (bits64) ( aSig0<<1 )
              || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
@@ -3879,7 +4502,7 @@ floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM )
 | Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
@@ -3909,7 +4532,7 @@ flag floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM )
 | Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -3942,7 +4565,7 @@ flag floatx80_le( floatx80 a, floatx80 b STATUS_PARAM )
 | Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -3975,7 +4598,7 @@ flag floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
@@ -4002,7 +4625,7 @@ flag floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM )
 | to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -4038,7 +4661,7 @@ flag floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM )
 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
+int floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -4272,7 +4895,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 );
     }
@@ -4306,7 +4929,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 );
     }
@@ -4341,7 +4964,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 ) );
     }
@@ -4521,7 +5144,10 @@ static float128 addFloat128Sigs( float128 a, float128 b, flag zSign STATUS_PARAM
             return a;
         }
         add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
-        if ( aExp == 0 ) return packFloat128( zSign, 0, zSig0, zSig1 );
+        if ( aExp == 0 ) {
+            if ( STATUS(flush_to_zero) ) return packFloat128( zSign, 0, 0, 0 );
+            return packFloat128( zSign, 0, zSig0, zSig1 );
+        }
         zSig2 = 0;
         zSig0 |= LIT64( 0x0002000000000000 );
         zExp = aExp;
@@ -4822,7 +5448,7 @@ float128 float128_div( float128 a, float128 b STATUS_PARAM )
 
 float128 float128_rem( float128 a, float128 b STATUS_PARAM )
 {
-    flag aSign, bSign, zSign;
+    flag aSign, zSign;
     int32 aExp, bExp, expDiff;
     bits64 aSig0, aSig1, bSig0, bSig1, q, term0, term1, term2;
     bits64 allZero, alternateASig0, alternateASig1, sigMean1;
@@ -4836,7 +5462,6 @@ float128 float128_rem( float128 a, float128 b STATUS_PARAM )
     bSig1 = extractFloat128Frac1( b );
     bSig0 = extractFloat128Frac0( b );
     bExp = extractFloat128Exp( b );
-    bSign = extractFloat128Sign( b );
     if ( aExp == 0x7FFF ) {
         if (    ( aSig0 | aSig1 )
              || ( ( bExp == 0x7FFF ) && ( bSig0 | bSig1 ) ) ) {
@@ -4911,7 +5536,7 @@ float128 float128_rem( float128 a, float128 b STATUS_PARAM )
         sub128( aSig0, aSig1, bSig0, bSig1, &aSig0, &aSig1 );
     } while ( 0 <= (sbits64) aSig0 );
     add128(
-        aSig0, aSig1, alternateASig0, alternateASig1, &sigMean0, &sigMean1 );
+        aSig0, aSig1, alternateASig0, alternateASig1, (bits64 *)&sigMean0, &sigMean1 );
     if (    ( sigMean0 < 0 )
          || ( ( ( sigMean0 | sigMean1 ) == 0 ) && ( q & 1 ) ) ) {
         aSig0 = alternateASig0;
@@ -4999,7 +5624,7 @@ float128 float128_sqrt( float128 a STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float128_eq( float128 a, float128 b STATUS_PARAM )
+int float128_eq( float128 a, float128 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
@@ -5029,7 +5654,7 @@ flag float128_eq( float128 a, float128 b STATUS_PARAM )
 | Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float128_le( float128 a, float128 b STATUS_PARAM )
+int float128_le( float128 a, float128 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -5061,7 +5686,7 @@ flag float128_le( float128 a, float128 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float128_lt( float128 a, float128 b STATUS_PARAM )
+int float128_lt( float128 a, float128 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -5094,7 +5719,7 @@ flag float128_lt( float128 a, float128 b STATUS_PARAM )
 | according to the IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
+int float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
 {
 
     if (    (    ( extractFloat128Exp( a ) == 0x7FFF )
@@ -5121,7 +5746,7 @@ flag float128_eq_signaling( float128 a, float128 b STATUS_PARAM )
 | IEC/IEEE Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float128_le_quiet( float128 a, float128 b STATUS_PARAM )
+int float128_le_quiet( float128 a, float128 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -5157,7 +5782,7 @@ flag float128_le_quiet( float128 a, float128 b STATUS_PARAM )
 | Standard for Binary Floating-Point Arithmetic.
 *----------------------------------------------------------------------------*/
 
-flag float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
+int float128_lt_quiet( float128 a, float128 b STATUS_PARAM )
 {
     flag aSign, bSign;
 
@@ -5235,6 +5860,24 @@ 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 )
+{
+    int64_t v;
+    unsigned int res;
+
+    v = float32_to_int64_round_to_zero(a STATUS_VAR);
+    if (v < 0) {
+        res = 0;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else if (v > 0xffff) {
+        res = 0xffff;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else {
+        res = v;
+    }
+    return res;
+}
+
 unsigned int float64_to_uint32( float64 a STATUS_PARAM )
 {
     int64_t v;
@@ -5271,11 +5914,55 @@ 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 )
+{
+    int64_t v;
+    unsigned int res;
+
+    v = float64_to_int64_round_to_zero(a STATUS_VAR);
+    if (v < 0) {
+        res = 0;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else if (v > 0xffff) {
+        res = 0xffff;
+        float_raise( float_flag_invalid STATUS_VAR);
+    } else {
+        res = v;
+    }
+    return res;
+}
+
+/* FIXME: This looks broken.  */
+uint64_t float64_to_uint64 (float64 a STATUS_PARAM)
+{
+    int64_t v;
+
+    v = float64_val(int64_to_float64(INT64_MIN STATUS_VAR));
+    v += float64_val(a);
+    v = float64_to_int64(make_float64(v) STATUS_VAR);
+
+    return v - INT64_MIN;
+}
+
+uint64_t float64_to_uint64_round_to_zero (float64 a STATUS_PARAM)
+{
+    int64_t v;
+
+    v = float64_val(int64_to_float64(INT64_MIN STATUS_VAR));
+    v += float64_val(a);
+    v = float64_to_int64_round_to_zero(make_float64(v) STATUS_VAR);
+
+    return v - INT64_MIN;
+}
+
 #define COMPARE(s, nan_exp)                                                  \
-INLINE char float ## s ## _compare_internal( float ## s a, float ## s b,     \
+INLINE int float ## s ## _compare_internal( float ## s a, float ## s b,      \
                                       int is_quiet STATUS_PARAM )            \
 {                                                                            \
     flag aSign, bSign;                                                       \
+    bits ## s 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 ) ) ||                                 \
@@ -5290,31 +5977,175 @@ INLINE char float ## s ## _compare_internal( float ## s a, float ## s b,     \
     }                                                                        \
     aSign = extractFloat ## s ## Sign( a );                                  \
     bSign = extractFloat ## s ## Sign( b );                                  \
+    av = float ## s ## _val(a);                                              \
+    bv = float ## s ## _val(b);                                              \
     if ( aSign != bSign ) {                                                  \
-        if ( (bits ## s) ( ( a | b )<<1 ) == 0 ) {                           \
+        if ( (bits ## s) ( ( av | bv )<<1 ) == 0 ) {                         \
             /* zero case */                                                  \
             return float_relation_equal;                                     \
         } else {                                                             \
             return 1 - (2 * aSign);                                          \
         }                                                                    \
     } else {                                                                 \
-        if (a == b) {                                                        \
+        if (av == bv) {                                                      \
             return float_relation_equal;                                     \
         } else {                                                             \
-            return 1 - 2 * (aSign ^ ( a < b ));                              \
+            return 1 - 2 * (aSign ^ ( av < bv ));                            \
         }                                                                    \
     }                                                                        \
 }                                                                            \
                                                                              \
-char float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM )       \
+int float ## s ## _compare( float ## s a, float ## s b STATUS_PARAM )        \
 {                                                                            \
     return float ## s ## _compare_internal(a, b, 0 STATUS_VAR);              \
 }                                                                            \
                                                                              \
-char float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM ) \
+int float ## s ## _compare_quiet( float ## s a, float ## s b STATUS_PARAM )  \
 {                                                                            \
     return float ## s ## _compare_internal(a, b, 1 STATUS_VAR);              \
 }
 
 COMPARE(32, 0xff)
 COMPARE(64, 0x7ff)
+
+INLINE int float128_compare_internal( float128 a, float128 b,
+                                      int is_quiet STATUS_PARAM )
+{
+    flag aSign, bSign;
+
+    if (( ( extractFloat128Exp( a ) == 0x7fff ) &&
+          ( extractFloat128Frac0( a ) | extractFloat128Frac1( a ) ) ) ||
+        ( ( extractFloat128Exp( b ) == 0x7fff ) &&
+          ( extractFloat128Frac0( b ) | extractFloat128Frac1( b ) ) )) {
+        if (!is_quiet ||
+            float128_is_signaling_nan( a ) ||
+            float128_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid STATUS_VAR);
+        }
+        return float_relation_unordered;
+    }
+    aSign = extractFloat128Sign( a );
+    bSign = extractFloat128Sign( b );
+    if ( aSign != bSign ) {
+        if ( ( ( ( a.high | b.high )<<1 ) | 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 float128_compare( float128 a, float128 b STATUS_PARAM )
+{
+    return float128_compare_internal(a, b, 0 STATUS_VAR);
+}
+
+int float128_compare_quiet( float128 a, float128 b STATUS_PARAM )
+{
+    return float128_compare_internal(a, b, 1 STATUS_VAR);
+}
+
+/* Multiply A by 2 raised to the power N.  */
+float32 float32_scalbn( float32 a, int n STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp;
+    bits32 aSig;
+
+    a = float32_squash_input_denormal(a STATUS_VAR);
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+
+    if ( aExp == 0xFF ) {
+        return a;
+    }
+    if ( aExp != 0 )
+        aSig |= 0x00800000;
+    else if ( aSig == 0 )
+        return a;
+
+    aExp += n - 1;
+    aSig <<= 7;
+    return normalizeRoundAndPackFloat32( aSign, aExp, aSig STATUS_VAR );
+}
+
+float64 float64_scalbn( float64 a, int n STATUS_PARAM )
+{
+    flag aSign;
+    int16 aExp;
+    bits64 aSig;
+
+    a = float64_squash_input_denormal(a STATUS_VAR);
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+
+    if ( aExp == 0x7FF ) {
+        return a;
+    }
+    if ( aExp != 0 )
+        aSig |= LIT64( 0x0010000000000000 );
+    else if ( aSig == 0 )
+        return a;
+
+    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;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+
+    if ( aExp == 0x7FF ) {
+        return a;
+    }
+    if (aExp == 0 && aSig == 0)
+        return a;
+
+    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;
+
+    aSig1 = extractFloat128Frac1( a );
+    aSig0 = extractFloat128Frac0( a );
+    aExp = extractFloat128Exp( a );
+    aSign = extractFloat128Sign( a );
+    if ( aExp == 0x7FFF ) {
+        return a;
+    }
+    if ( aExp != 0 )
+        aSig0 |= LIT64( 0x0001000000000000 );
+    else if ( aSig0 == 0 && aSig1 == 0 )
+        return a;
+
+    aExp += n - 1;
+    return normalizeRoundAndPackFloat128( aSign, aExp, aSig0, aSig1
+                                          STATUS_VAR );
+
+}
+#endif