FIX #2271
authorDaniel Fischer <daniel.is.fischer@web.de>
Mon, 18 Oct 2010 21:03:37 +0000 (21:03 +0000)
committerDaniel Fischer <daniel.is.fischer@web.de>
Mon, 18 Oct 2010 21:03:37 +0000 (21:03 +0000)
Faster rounding functions for Double and float with Int or Integer results.
Fixes #2271.
Since some glibc's have buggy rintf or rint functions and the behaviour of
these functions depends on the setting of the rounding mode, we provide our
own implementations which always round ties to even.

Also added rewrite rules and removed trailing whitespace.

GHC/Float.lhs
GHC/Float/RealFracMethods.hs [new file with mode: 0644]
base.cabal
cbits/primFloat.c

index 02aba8c..a9230c2 100644 (file)
@@ -27,7 +27,8 @@
 #include "ieee-flpt.h"
 
 -- #hide
-module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double# )
+module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double#
+                , double2Int, int2Double, float2Int, int2Float )
     where
 
 import Data.Maybe
@@ -40,6 +41,7 @@ import GHC.Show
 import GHC.Num
 import GHC.Real
 import GHC.Arr
+import GHC.Float.RealFracMethods
 
 infixr 8  **
 \end{code}
@@ -197,19 +199,25 @@ instance  Fractional Float  where
     fromRational x      =  fromRat x
     recip x             =  1.0 / x
 
-{-# RULES "truncate/Float->Int" truncate = float2Int #-}
+-- RULES for Integer and Int
+{-# RULES
+"properFraction/Float->Integer"     properFraction = properFractionFloatInteger
+"truncate/Float->Integer"           truncate = truncateFloatInteger
+"floor/Float->Integer"              floor = floorFloatInteger
+"ceiling/Float->Integer"            ceiling = ceilingFloatInteger
+"round/Float->Integer"              round = roundFloatInteger
+"properFraction/Float->Int"         properFraction = properFractionFloatInt
+"truncate/Float->Int"               truncate = float2Int
+"floor/Float->Int"                  floor = floorFloatInt
+"ceiling/Float->Int"                ceiling = ceilingFloatInt
+"round/Float->Int"                  round = roundFloatInt
+  #-}
 instance  RealFrac Float  where
 
-    {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
-    {-# SPECIALIZE round    :: Float -> Int #-}
-
-    {-# SPECIALIZE properFraction :: Float  -> (Integer, Float) #-}
-    {-# SPECIALIZE round    :: Float -> Integer #-}
-
         -- ceiling, floor, and truncate are all small
-    {-# INLINE ceiling #-}
-    {-# INLINE floor #-}
-    {-# INLINE truncate #-}
+    {-# INLINE [1] ceiling #-}
+    {-# INLINE [1] floor #-}
+    {-# INLINE [1] truncate #-}
 
 -- We assume that FLT_RADIX is 2 so that we can use more efficient code
 #if FLT_RADIX != 2
@@ -352,27 +360,32 @@ instance  Floating Double  where
     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 
-{-# RULES "truncate/Double->Int" truncate = double2Int #-}
+-- RULES for Integer and Int
+{-# RULES
+"properFraction/Double->Integer"    properFraction = properFractionDoubleInteger
+"truncate/Double->Integer"          truncate = truncateDoubleInteger
+"floor/Double->Integer"             floor = floorDoubleInteger
+"ceiling/Double->Integer"           ceiling = ceilingDoubleInteger
+"round/Double->Integer"             round = roundDoubleInteger
+"properFraction/Double->Int"        properFraction = properFractionDoubleInt
+"truncate/Double->Int"              truncate = double2Int
+"floor/Double->Int"                 floor = floorDoubleInt
+"ceiling/Double->Int"               ceiling = ceilingDoubleInt
+"round/Double->Int"                 round = roundDoubleInt
+  #-}
 instance  RealFrac Double  where
 
-    {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
-    {-# SPECIALIZE round    :: Double -> Int #-}
-
-    {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
-    {-# SPECIALIZE round    :: Double -> Integer #-}
-
         -- ceiling, floor, and truncate are all small
-    {-# INLINE ceiling #-}
-    {-# INLINE floor #-}
-    {-# INLINE truncate #-}
+    {-# INLINE [1] ceiling #-}
+    {-# INLINE [1] floor #-}
+    {-# INLINE [1] truncate #-}
 
     properFraction x
       = case (decodeFloat x)      of { (m,n) ->
-        let  b = floatRadix x     in
         if n >= 0 then
-            (fromInteger m * fromInteger b ^ n, 0.0)
+            (fromInteger m * 2 ^ n, 0.0)
         else
-            case (quotRem m (b^(negate n))) of { (w,r) ->
+            case (quotRem m (2^(negate n))) of { (w,r) ->
             (fromInteger w, encodeFloat r n)
             }
         }
@@ -851,12 +864,6 @@ neFloat     (F# x) (F# y) = neFloat# x y
 ltFloat     (F# x) (F# y) = ltFloat# x y
 leFloat     (F# x) (F# y) = leFloat# x y
 
-float2Int :: Float -> Int
-float2Int   (F# x) = I# (float2Int# x)
-
-int2Float :: Int -> Float
-int2Float   (I# x) = F# (int2Float# x)
-
 expFloat, logFloat, sqrtFloat :: Float -> Float
 sinFloat, cosFloat, tanFloat  :: Float -> Float
 asinFloat, acosFloat, atanFloat  :: Float -> Float
@@ -897,12 +904,6 @@ neDouble    (D# x) (D# y) = x /=## y
 ltDouble    (D# x) (D# y) = x <## y
 leDouble    (D# x) (D# y) = x <=## y
 
-double2Int :: Double -> Int
-double2Int   (D# x) = I# (double2Int#   x)
-
-int2Double :: Int -> Double
-int2Double   (I# x) = D# (int2Double#   x)
-
 double2Float :: Double -> Float
 double2Float (D# x) = F# (double2Float# x)
 
diff --git a/GHC/Float/RealFracMethods.hs b/GHC/Float/RealFracMethods.hs
new file mode 100644 (file)
index 0000000..1728e75
--- /dev/null
@@ -0,0 +1,341 @@
+{-# LANGUAGE CPP, MagicHash, UnboxedTuples, ForeignFunctionInterface,
+    NoImplicitPrelude #-}
+{-# OPTIONS_HADDOCK hide #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  GHC.Float.RealFracMethods
+-- Copyright   :  (c) Daniel Fischer 2010
+-- License     :  see libraries/base/LICENSE
+--
+-- Maintainer  :  cvs-ghc@haskell.org
+-- Stability   :  internal
+-- Portability :  non-portable (GHC Extensions)
+--
+-- Methods for the RealFrac instances for 'Float' and 'Double',
+-- with specialised versions for 'Int'.
+--
+-- Moved to their own module to not bloat GHC.Float further.
+--
+-----------------------------------------------------------------------------
+
+#include "MachDeps.h"
+
+-- #hide
+module GHC.Float.RealFracMethods
+    ( -- * Double methods
+      -- ** Integer results
+      properFractionDoubleInteger
+    , truncateDoubleInteger
+    , floorDoubleInteger
+    , ceilingDoubleInteger
+    , roundDoubleInteger
+      -- ** Int results
+    , properFractionDoubleInt
+    , floorDoubleInt
+    , ceilingDoubleInt
+    , roundDoubleInt
+      -- * Double/Int conversions, wrapped primops
+    , double2Int
+    , int2Double
+      -- * Float methods
+      -- ** Integer results
+    , properFractionFloatInteger
+    , truncateFloatInteger
+    , floorFloatInteger
+    , ceilingFloatInteger
+    , roundFloatInteger
+      -- ** Int results
+    , properFractionFloatInt
+    , floorFloatInt
+    , ceilingFloatInt
+    , roundFloatInt
+      -- * Float/Int conversions, wrapped primops
+    , float2Int
+    , int2Float
+    ) where
+
+import GHC.Integer
+
+import GHC.Base
+import GHC.Num
+
+#if WORD_SIZE_IN_BITS < 64
+
+import GHC.IntWord64
+
+#define TO64 integerToInt64
+#define FROM64 int64ToInteger
+#define MINUS64 minusInt64#
+#define NEGATE64 negateInt64#
+
+#else
+
+#define TO64 toInt#
+#define FROM64 smallInteger
+#define MINUS64 ( -# )
+#define NEGATE64 negateInt#
+
+uncheckedIShiftRA64# :: Int# -> Int# -> Int#
+uncheckedIShiftRA64# = uncheckedIShiftRA#
+
+uncheckedIShiftL64# :: Int# -> Int# -> Int#
+uncheckedIShiftL64# = uncheckedIShiftL#
+
+#endif
+
+default ()
+
+------------------------------------------------------------------------------
+--                              Float Methods                               --
+------------------------------------------------------------------------------
+
+-- Special Functions for Int, nice, easy and fast.
+-- They should be small enough to be inlined automatically.
+
+-- We have to test for ±0.0 to avoid returning -0.0 in the second
+-- component of the pair. Unfortunately the branching costs a lot
+-- of performance.
+properFractionFloatInt :: Float -> (Int, Float)
+properFractionFloatInt (F# x) =
+    if x `eqFloat#` 0.0#
+        then (I# 0#, F# 0.0#)
+        else case float2Int# x of
+                n -> (I# n, F# (x `minusFloat#` int2Float# n))
+
+-- truncateFloatInt = float2Int
+
+floorFloatInt :: Float -> Int
+floorFloatInt (F# x) =
+    case float2Int# x of
+      n | x `ltFloat#` int2Float# n -> I# (n -# 1#)
+        | otherwise                 -> I# n
+
+ceilingFloatInt :: Float -> Int
+ceilingFloatInt (F# x) =
+    case float2Int# x of
+      n | int2Float# n `ltFloat#` x  -> I# (n +# 1#)
+        | otherwise                 -> I# n
+
+roundFloatInt :: Float -> Int
+roundFloatInt x = float2Int (c_rintFloat x)
+
+-- Functions with Integer results
+
+-- With the new code generator in GHC 7, the explicit bit-fiddling is
+-- slower than the old code for values of small modulus, but when the
+-- 'Int' range is left, the bit-fiddling quickly wins big, so we use that.
+-- If the methods are called on smallish values, hopefully people go
+-- through Int and not larger types.
+
+-- Note: For negative exponents, we must check the validity of the shift
+-- distance for the right shifts of the mantissa.
+
+{-# INLINE properFractionFloatInteger #-}
+properFractionFloatInteger :: Float -> (Integer, Float)
+properFractionFloatInteger v@(F# x) =
+    case decodeFloat_Int# x of
+      (# m, e #)
+        | e <# 0#   ->
+          case negateInt# e of
+            s | s ># 23#    -> (0, v)
+              | m <# 0#     ->
+                case negateInt# (negateInt# m `uncheckedIShiftRA#` s) of
+                  k -> (smallInteger k,
+                            case m -# (k `uncheckedIShiftL#` s) of
+                              r -> F# (encodeFloatInteger (smallInteger r) e))
+              | otherwise           ->
+                case m `uncheckedIShiftRL#` s of
+                  k -> (smallInteger k,
+                            case m -# (k `uncheckedIShiftL#` s) of
+                              r -> F# (encodeFloatInteger (smallInteger r) e))
+        | otherwise -> (shiftLInteger (smallInteger m) e, F# 0.0#)
+
+{-# INLINE truncateFloatInteger #-}
+truncateFloatInteger :: Float -> Integer
+truncateFloatInteger x =
+    case properFractionFloatInteger x of
+      (n, _) -> n
+
+-- floor is easier for negative numbers than truncate, so this gets its
+-- own implementation, it's a little faster.
+{-# INLINE floorFloatInteger #-}
+floorFloatInteger :: Float -> Integer
+floorFloatInteger (F# x) =
+    case decodeFloat_Int# x of
+      (# m, e #)
+        | e <# 0#   ->
+          case negateInt# e of
+            s | s ># 23#    -> if m <# 0# then (-1) else 0
+              | otherwise   -> smallInteger (m `uncheckedIShiftRA#` s)
+        | otherwise -> shiftLInteger (smallInteger m) e
+
+-- ceiling x = -floor (-x)
+-- If giving this its own implementation is faster at all,
+-- it's only marginally so, hence we keep it short.
+{-# INLINE ceilingFloatInteger #-}
+ceilingFloatInteger :: Float -> Integer
+ceilingFloatInteger (F# x) =
+    negateInteger (floorFloatInteger (F# (negateFloat# x)))
+
+{-# INLINE roundFloatInteger #-}
+roundFloatInteger :: Float -> Integer
+roundFloatInteger x = float2Integer (c_rintFloat x)
+
+------------------------------------------------------------------------------
+--                              Double Methods                              --
+------------------------------------------------------------------------------
+
+-- Special Functions for Int, nice, easy and fast.
+-- They should be small enough to be inlined automatically.
+
+-- We have to test for ±0.0 to avoid returning -0.0 in the second
+-- component of the pair. Unfortunately the branching costs a lot
+-- of performance.
+properFractionDoubleInt :: Double -> (Int, Double)
+properFractionDoubleInt (D# x) =
+    if x ==## 0.0##
+        then (I# 0#, D# 0.0##)
+        else case double2Int# x of
+                n -> (I# n, D# (x -## int2Double# n))
+
+-- truncateDoubleInt = double2Int
+
+floorDoubleInt :: Double -> Int
+floorDoubleInt (D# x) =
+    case double2Int# x of
+      n | x <## int2Double# n   -> I# (n -# 1#)
+        | otherwise             -> I# n
+
+ceilingDoubleInt :: Double -> Int
+ceilingDoubleInt (D# x) =
+    case double2Int# x of
+      n | int2Double# n <## x   -> I# (n +# 1#)
+        | otherwise             -> I# n
+
+roundDoubleInt :: Double -> Int
+roundDoubleInt x = double2Int (c_rintDouble x)
+
+-- Functions with Integer results
+
+-- The new Code generator isn't quite as good for the old 'Double' code
+-- as for the 'Float' code, so for 'Double' the bit-fiddling also wins
+-- when the values have small modulus.
+
+-- When the exponent is negative, all mantissae have less than 64 bits
+-- and the right shifting of sized types is much faster than that of
+-- 'Integer's, especially when we can
+
+-- Note: For negative exponents, we must check the validity of the shift
+-- distance for the right shifts of the mantissa.
+
+{-# INLINE properFractionDoubleInteger #-}
+properFractionDoubleInteger :: Double -> (Integer, Double)
+properFractionDoubleInteger v@(D# x) =
+    case decodeDoubleInteger x of
+      (# m, e #)
+        | e <# 0#   ->
+          case negateInt# e of
+            s | s ># 52#    -> (0, v)
+              | m < 0       ->
+                case TO64 (negateInteger m) of
+                  n ->
+                    case n `uncheckedIShiftRA64#` s of
+                      k ->
+                        (FROM64 (NEGATE64 k),
+                          case MINUS64 n (k `uncheckedIShiftL64#` s) of
+                            r ->
+                              D# (encodeDoubleInteger (FROM64 (NEGATE64 r)) e))
+              | otherwise           ->
+                case TO64 m of
+                  n ->
+                    case n `uncheckedIShiftRA64#` s of
+                      k -> (FROM64 k,
+                            case MINUS64 n (k `uncheckedIShiftL64#` s) of
+                              r -> D# (encodeDoubleInteger (FROM64 r) e))
+        | otherwise -> (shiftLInteger m e, D# 0.0##)
+
+{-# INLINE truncateDoubleInteger #-}
+truncateDoubleInteger :: Double -> Integer
+truncateDoubleInteger x =
+    case properFractionDoubleInteger x of
+      (n, _) -> n
+
+-- floor is easier for negative numbers than truncate, so this gets its
+-- own implementation, it's a little faster.
+{-# INLINE floorDoubleInteger #-}
+floorDoubleInteger :: Double -> Integer
+floorDoubleInteger (D# x) =
+    case decodeDoubleInteger x of
+      (# m, e #)
+        | e <# 0#   ->
+          case negateInt# e of
+            s | s ># 52#    -> if m < 0 then (-1) else 0
+              | otherwise   ->
+                case TO64 m of
+                  n -> FROM64 (n `uncheckedIShiftRA64#` s)
+        | otherwise -> shiftLInteger m e
+
+{-# INLINE ceilingDoubleInteger #-}
+ceilingDoubleInteger :: Double -> Integer
+ceilingDoubleInteger (D# x) =
+    negateInteger (floorDoubleInteger (D# (negateDouble# x)))
+
+{-# INLINE roundDoubleInteger #-}
+roundDoubleInteger :: Double -> Integer
+roundDoubleInteger x = double2Integer (c_rintDouble x)
+
+-- Wrappers around double2Int#, int2Double#, float2Int# and int2Float#,
+-- we need them here, so we move them from GHC.Float and re-export them
+-- explicitly from there.
+
+double2Int :: Double -> Int
+double2Int (D# x) = I# (double2Int# x)
+
+int2Double :: Int -> Double
+int2Double (I# i) = D# (int2Double# i)
+
+float2Int :: Float -> Int
+float2Int (F# x) = I# (float2Int# x)
+
+int2Float :: Int -> Float
+int2Float (I# i) = F# (int2Float# i)
+
+-- Quicker conversions from 'Double' and 'Float' to 'Integer',
+-- assuming the floating point value is integral.
+--
+-- Note: Since the value is integral, the exponent can't be less than
+-- (-TYP_MANT_DIG), so we need not check the validity of the shift
+-- distance for the right shfts here.
+
+{-# INLINE double2Integer #-}
+double2Integer :: Double -> Integer
+double2Integer (D# x) =
+    case decodeDoubleInteger x of
+      (# m, e #)
+        | e <# 0#   ->
+          case TO64 m of
+            n -> FROM64 (n `uncheckedIShiftRA64#` negateInt# e)
+        | otherwise -> shiftLInteger m e
+
+{-# INLINE float2Integer #-}
+float2Integer :: Float -> Integer
+float2Integer (F# x) =
+    case decodeFloat_Int# x of
+      (# m, e #)
+        | e <# 0#   -> smallInteger (m `uncheckedIShiftRA#` negateInt# e)
+        | otherwise -> shiftLInteger (smallInteger m) e
+
+-- Foreign imports, the rounding is done faster in C when the value
+-- isn't integral, so we call out for rounding. For values of large
+-- modulus, calling out to C is slower than staying in Haskell, but
+-- presumably 'round' is mostly called for values with smaller modulus,
+-- when calling out to C is a major win.
+-- For all other functions, calling out to C gives at most a marginal
+-- speedup for values of small modulus and is much slower than staying
+-- in Haskell for values of large modulus, so those are done in Haskell.
+
+foreign import ccall unsafe "rintDouble"
+    c_rintDouble :: Double -> Double
+
+foreign import ccall unsafe "rintFloat"
+    c_rintFloat :: Float -> Float
index c4eadf7..d863d7a 100644 (file)
@@ -52,6 +52,7 @@ Library {
             GHC.Exception,
             GHC.Exts,
             GHC.Float,
+            GHC.Float.RealFracMethods,
             GHC.ForeignPtr,
             GHC.MVar,
             GHC.IO,
index 3fa39d3..a8f4803 100644 (file)
@@ -44,7 +44,7 @@ union stg_ieee754_flt
 };
 
 /*
+
  To recap, here's the representation of a double precision
  IEEE floating point number:
 
@@ -106,7 +106,7 @@ union stg_ieee754_dbl
 
 /*
  * Predicates for testing for extreme IEEE fp values.
- */ 
+ */
 
 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
 #ifdef IEEE_FLOATING_POINT
@@ -115,7 +115,7 @@ HsInt
 isDoubleNaN(HsDouble d)
 {
   union stg_ieee754_dbl u;
-  
+
   u.d = d;
 
   return (
@@ -141,7 +141,7 @@ isDoubleInfinite(HsDouble d)
 }
 
 HsInt
-isDoubleDenormalized(HsDouble d) 
+isDoubleDenormalized(HsDouble d)
 {
     union stg_ieee754_dbl u;
 
@@ -154,16 +154,16 @@ isDoubleDenormalized(HsDouble d)
         - (don't care about setting of sign bit.)
 
     */
-    return (  
+    return (
        u.ieee.exponent  == 0 &&
        (u.ieee.mantissa0 != 0 ||
         u.ieee.mantissa1 != 0)
       );
-        
+
 }
 
 HsInt
-isDoubleNegativeZero(HsDouble d) 
+isDoubleNegativeZero(HsDouble d)
 {
     union stg_ieee754_dbl u;
 
@@ -207,7 +207,7 @@ isFloatInfinite(HsFloat f)
 {
     union stg_ieee754_flt u;
     u.f = f;
-  
+
     /* A float is Inf iff exponent is max (all ones),
        and mantissa is min(all zeros.) */
     return (
@@ -234,7 +234,7 @@ isFloatDenormalized(HsFloat f)
 }
 
 HsInt
-isFloatNegativeZero(HsFloat f) 
+isFloatNegativeZero(HsFloat f)
 {
     union stg_ieee754_flt u;
     u.f = f;
@@ -246,6 +246,184 @@ isFloatNegativeZero(HsFloat f)
        u.ieee.mantissa == 0);
 }
 
+/*
+ There are glibc versions around with buggy rintf or rint, hence we
+ provide our own. We always round ties to even, so we can be simpler.
+*/
+
+#define FLT_HIDDEN 0x800000
+#define FLT_POWER2 0x1000000
+
+HsFloat
+rintFloat(HsFloat f)
+{
+    union stg_ieee754_flt u;
+    u.f = f;
+    /* if real exponent > 22, it's already integral, infinite or nan */
+    if (u.ieee.exponent > 149)  /* 22 + 127 */
+    {
+        return u.f;
+    }
+    if (u.ieee.exponent < 126)  /* (-1) + 127, abs(f) < 0.5 */
+    {
+        /* only used for rounding to Integral a, so don't care about -0.0 */
+        return 0.0;
+    }
+    /* 0.5 <= abs(f) < 2^23 */
+    unsigned int half, mask, mant, frac;
+    half = 1 << (149 - u.ieee.exponent);    /* bit for 0.5 */
+    mask = 2*half - 1;                      /* fraction bits */
+    mant = u.ieee.mantissa | FLT_HIDDEN;    /* add hidden bit */
+    frac = mant & mask;                     /* get fraction */
+    mant ^= frac;                           /* truncate mantissa */
+    if ((frac < half) || ((frac == half) && ((mant & (2*half)) == 0)))
+    {
+        /* this means we have to truncate */
+        if (mant == 0)
+        {
+            /* f == ±0.5, return 0.0 */
+            return 0.0;
+        }
+        else
+        {
+            /* remove hidden bit and set mantissa */
+            u.ieee.mantissa = mant ^ FLT_HIDDEN;
+            return u.f;
+        }
+    }
+    else
+    {
+        /* round away from zero, increment mantissa */
+        mant += 2*half;
+        if (mant == FLT_POWER2)
+        {
+            /* next power of 2, increase exponent an set mantissa to 0 */
+            u.ieee.mantissa = 0;
+            u.ieee.exponent += 1;
+            return u.f;
+        }
+        else
+        {
+            /* remove hidden bit and set mantissa */
+            u.ieee.mantissa = mant ^ FLT_HIDDEN;
+            return u.f;
+        }
+    }
+}
+
+#define DBL_HIDDEN 0x100000
+#define DBL_POWER2 0x200000
+#define LTOP_BIT 0x80000000
+
+HsDouble
+rintDouble(HsDouble d)
+{
+    union stg_ieee754_dbl u;
+    u.d = d;
+    /* if real exponent > 51, it's already integral, infinite or nan */
+    if (u.ieee.exponent > 1074) /* 51 + 1023 */
+    {
+        return u.d;
+    }
+    if (u.ieee.exponent < 1022)  /* (-1) + 1023, abs(d) < 0.5 */
+    {
+        /* only used for rounding to Integral a, so don't care about -0.0 */
+        return 0.0;
+    }
+    unsigned int half, mask, mant, frac;
+    if (u.ieee.exponent < 1043) /* 20 + 1023, real exponent < 20 */
+    {
+        /* the fractional part meets the higher part of the mantissa */
+        half = 1 << (1042 - u.ieee.exponent);   /* bit for 0.5 */
+        mask = 2*half - 1;                      /* fraction bits */
+        mant = u.ieee.mantissa0 | DBL_HIDDEN;   /* add hidden bit */
+        frac = mant & mask;                     /* get fraction */
+        mant ^= frac;                           /* truncate mantissa */
+        if ((frac < half) ||
+            ((frac == half) && (u.ieee.mantissa1 == 0)  /* a tie */
+                && ((mant & (2*half)) == 0)))
+        {
+            /* truncate */
+            if (mant == 0)
+            {
+                /* d = ±0.5, return 0.0 */
+                return 0.0;
+            }
+            /* remove hidden bit and set mantissa */
+            u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
+            u.ieee.mantissa1 = 0;
+            return u.d;
+        }
+        else    /* round away from zero */
+        {
+            /* zero low mantissa bits */
+            u.ieee.mantissa1 = 0;
+            /* increment integer part of mantissa */
+            mant += 2*half;
+            if (mant == DBL_POWER2)
+            {
+                /* power of 2, increment exponent and zero mantissa */
+                u.ieee.mantissa0 = 0;
+                u.ieee.exponent += 1;
+                return u.d;
+            }
+            /* remove hidden bit */
+            u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
+            return u.d;
+        }
+    }
+    else
+    {
+        /* 20 <= real exponent < 52, fractional part entirely in mantissa1 */
+        half = 1 << (1074 - u.ieee.exponent);   /* bit for 0.5 */
+        mask = 2*half - 1;                      /* fraction bits */
+        mant = u.ieee.mantissa1;                /* no hidden bit here */
+        frac = mant & mask;                     /* get fraction */
+        mant ^= frac;                           /* truncate mantissa */
+        if ((frac < half) ||
+            ((frac == half) &&                  /* tie */
+            (((half == LTOP_BIT) ? (u.ieee.mantissa0 & 1)  /* yuck */
+                                : (mant & (2*half)))
+                                        == 0)))
+        {
+            /* truncate */
+            u.ieee.mantissa1 = mant;
+            return u.d;
+        }
+        else
+        {
+            /* round away from zero */
+            /* increment mantissa */
+            mant += 2*half;
+            u.ieee.mantissa1 = mant;
+            if (mant == 0)
+            {
+                /* low part of mantissa overflowed */
+                /* increment high part of mantissa */
+                mant = u.ieee.mantissa0 + 1;
+                if (mant == DBL_HIDDEN)
+                {
+                    /* hit power of 2 */
+                    /* zero mantissa */
+                    u.ieee.mantissa0 = 0;
+                    /* and increment exponent */
+                    u.ieee.exponent += 1;
+                    return u.d;
+                }
+                else
+                {
+                    u.ieee.mantissa0 = mant;
+                    return u.d;
+                }
+            }
+            else
+            {
+                return u.d;
+            }
+        }
+    }
+}
+
 #else /* ! IEEE_FLOATING_POINT */
 
 /* Dummy definitions of predicates - they all return false */
@@ -258,4 +436,77 @@ HsInt isFloatInfinite(f) HsFloat f; { return 0; }
 HsInt isFloatDenormalized(f) HsFloat f; { return 0; }
 HsInt isFloatNegativeZero(f) HsFloat f; { return 0; }
 
+
+/* For exotic floating point formats, we can't do much */
+/* We suppose the format has not too many bits */
+/* I hope nobody tries to build GHC where this is wrong */
+
+#define FLT_UPP 536870912.0
+
+HsFloat
+rintFloat(HsFloat f)
+{
+    if ((f > FLT_UPP) || (f < (-FLT_UPP)))
+    {
+        return f;
+    }
+    else
+    {
+        int i = (int)f;
+        float g = i;
+        float d = f - g;
+        if (d > 0.5)
+        {
+            return g + 1.0;
+        }
+        if (d == 0.5)
+        {
+            return (i & 1) ? (g + 1.0) : g;
+        }
+        if (d == -0.5)
+        {
+            return (i & 1) ? (g - 1.0) : g;
+        }
+        if (d < -0.5)
+        {
+            return g - 1.0;
+        }
+        return g;
+    }
+}
+
+#define DBL_UPP 2305843009213693952.0
+
+HsDouble
+rintDouble(HsDouble d)
+{
+    if ((d > DBL_UPP) || (d < (-DBL_UPP)))
+    {
+        return d;
+    }
+    else
+    {
+        HsInt64 i = (HsInt64)d;
+        double e = i;
+        double r = d - e;
+        if (r > 0.5)
+        {
+            return e + 1.0;
+        }
+        if (r == 0.5)
+        {
+            return (i & 1) ? (e + 1.0) : e;
+        }
+        if (r == -0.5)
+        {
+            return (i & 1) ? (e - 1.0) : e;
+        }
+        if (r < -0.5)
+        {
+            return e - 1.0;
+        }
+        return e;
+    }
+}
+
 #endif /* ! IEEE_FLOATING_POINT */