add GHC.HetMet.{hetmet_kappa,hetmet_kappa_app}
[ghc-base.git] / GHC / Float.lhs
index aa2a534..1c6fd5f 100644 (file)
@@ -1,9 +1,15 @@
 \begin{code}
-{-# OPTIONS_GHC -XNoImplicitPrelude #-}
+{-# LANGUAGE CPP
+           , NoImplicitPrelude
+           , MagicHash
+           , UnboxedTuples
+           , ForeignFunctionInterface
+  #-}
 -- We believe we could deorphan this module, by moving lots of things
 -- around, but we haven't got there yet:
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 {-# OPTIONS_HADDOCK hide #-}
+
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  GHC.Float
@@ -21,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
@@ -34,6 +41,10 @@ import GHC.Show
 import GHC.Num
 import GHC.Real
 import GHC.Arr
+import GHC.Float.RealFracMethods
+import GHC.Float.ConversionUtils
+import GHC.Integer.Logarithms ( integerLogBase# )
+import GHC.Integer.Logarithms.Internals
 
 infixr 8  **
 \end{code}
@@ -182,28 +193,51 @@ instance  Num Float  where
     fromInteger i = F# (floatFromInteger i)
 
 instance  Real Float  where
-    toRational x        =  (m%1)*(b%1)^^n
-                           where (m,n) = decodeFloat x
-                                 b     = floatRadix  x
+    toRational (F# x#)  =
+        case decodeFloat_Int# x# of
+          (# m#, e# #)
+            | e# >=# 0#                                 ->
+                    (smallInteger m# `shiftLInteger` e#) :% 1
+            | (int2Word# m# `and#` 1##) `eqWord#` 0##   ->
+                    case elimZerosInt# m# (negateInt# e#) of
+                      (# n, d# #) -> n :% shiftLInteger 1 d#
+            | otherwise                                 ->
+                    smallInteger m# :% shiftLInteger 1 (negateInt# e#)
 
 instance  Fractional Float  where
     (/) x y             =  divideFloat x y
-    fromRational x      =  fromRat x
+    fromRational (n:%0)
+        | n == 0        = 0/0
+        | n < 0         = (-1)/0
+        | otherwise     = 1/0
+    fromRational (n:%d)
+        | n == 0        = encodeFloat 0 0
+        | n < 0         = -(fromRat'' minEx mantDigs (-n) d)
+        | otherwise     = fromRat'' minEx mantDigs n d
+          where
+            minEx       = FLT_MIN_EXP
+            mantDigs    = FLT_MANT_DIG
     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
@@ -316,13 +350,30 @@ instance  Num Double  where
 
 
 instance  Real Double  where
-    toRational x        =  (m%1)*(b%1)^^n
-                           where (m,n) = decodeFloat x
-                                 b     = floatRadix  x
+    toRational (D# x#)  =
+        case decodeDoubleInteger x# of
+          (# m, e# #)
+            | e# >=# 0#                                         ->
+                shiftLInteger m e# :% 1
+            | (int2Word# (toInt# m) `and#` 1##) `eqWord#` 0##   ->
+                case elimZerosInteger m (negateInt# e#) of
+                    (# n, d# #) ->  n :% shiftLInteger 1 d#
+            | otherwise                                         ->
+                m :% shiftLInteger 1 (negateInt# e#)
 
 instance  Fractional Double  where
     (/) x y             =  divideDouble x y
-    fromRational x      =  fromRat x
+    fromRational (n:%0)
+        | n == 0        = 0/0
+        | n < 0         = (-1)/0
+        | otherwise     = 1/0
+    fromRational (n:%d)
+        | n == 0        = encodeFloat 0 0
+        | n < 0         = -(fromRat'' minEx mantDigs (-n) d)
+        | otherwise     = fromRat'' minEx mantDigs n d
+          where
+            minEx       = DBL_MIN_EXP
+            mantDigs    = DBL_MANT_DIG
     recip x             =  1.0 / x
 
 instance  Floating Double  where
@@ -346,27 +397,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)
             }
         }
@@ -600,7 +656,7 @@ floatToDigits base x =
    if e >= 0 then
     let be = expt b e in
     if f == expt b (p-1) then
-      (f*be*b*2, 2*b, be*b, b)
+      (f*be*b*2, 2*b, be*b, be)     -- according to Burger and Dybvig
     else
       (f*be*2, 2, be, be)
    else
@@ -732,8 +788,10 @@ Now, here's Lennart's code (which works)
 
 \begin{code}
 -- | Converts a 'Rational' value into any type in class 'RealFloat'.
-{-# SPECIALISE fromRat :: Rational -> Double,
-                          Rational -> Float #-}
+{-# RULES
+"fromRat/Float"     fromRat = (fromRational :: Rational -> Float)
+"fromRat/Double"    fromRat = (fromRational :: Rational -> Double)
+  #-}
 fromRat :: (RealFloat a) => Rational -> a
 
 -- Deal with special cases first, delegating the real work to fromRat'
@@ -785,27 +843,106 @@ expt base n =
     if base == 2 && n >= minExpt && n <= maxExpt then
         expts!n
     else
-        base^n
+        if base == 10 && n <= maxExpt10 then
+            expts10!n
+        else
+            base^n
 
 expts :: Array Int Integer
 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
 
+maxExpt10 :: Int
+maxExpt10 = 324
+
+expts10 :: Array Int Integer
+expts10 = array (minExpt,maxExpt10) [(n,10^n) | n <- [minExpt .. maxExpt10]]
+
 -- Compute the (floor of the) log of i in base b.
 -- Simplest way would be just divide i by b until it's smaller then b, but that would
--- be very slow!  We are just slightly more clever.
+-- be very slow!  We are just slightly more clever, except for base 2, where
+-- we take advantage of the representation of Integers.
+-- The general case could be improved by a lookup table for
+-- approximating the result by integerLog2 i / integerLog2 b.
 integerLogBase :: Integer -> Integer -> Int
 integerLogBase b i
    | i < b     = 0
-   | otherwise = doDiv (i `div` (b^l)) l
-       where
-        -- Try squaring the base first to cut down the number of divisions.
-         l = 2 * integerLogBase (b*b) i
+   | b == 2    = I# (integerLog2# i)
+   | otherwise = I# (integerLogBase# b i)
+
+\end{code}
+
+Unfortunately, the old conversion code was awfully slow due to
+a) a slow integer logarithm
+b) repeated calculation of gcd's
+
+For the case of Rational's coming from a Float or Double via toRational,
+we can exploit the fact that the denominator is a power of two, which for
+these brings a huge speedup since we need only shift and add instead
+of division.
 
-         doDiv :: Integer -> Int -> Int
-         doDiv x y
-            | x < b     = y
-            | otherwise = doDiv (x `div` b) (y+1)
+The below is an adaption of fromRat' for the conversion to
+Float or Double exploiting the know floatRadix and avoiding
+divisions as much as possible.
 
+\begin{code}
+{-# SPECIALISE fromRat'' :: Int -> Int -> Integer -> Integer -> Float,
+                            Int -> Int -> Integer -> Integer -> Double #-}
+fromRat'' :: RealFloat a => Int -> Int -> Integer -> Integer -> a
+fromRat'' minEx@(I# me#) mantDigs@(I# md#) n d =
+    case integerLog2IsPowerOf2# d of
+      (# ld#, pw# #)
+        | pw# ==# 0# ->
+          case integerLog2# n of
+            ln# | ln# ># (ld# +# me#) ->
+                  if ln# <# md#
+                    then encodeFloat (n `shiftL` (I# (md# -# 1# -# ln#)))
+                                        (I# (ln# +# 1# -# ld# -# md#))
+                    else let n'  = n `shiftR` (I# (ln# +# 1# -# md#))
+                             n'' = case roundingMode# n (ln# -# md#) of
+                                    0# -> n'
+                                    2# -> n' + 1
+                                    _  -> case fromInteger n' .&. (1 :: Int) of
+                                            0 -> n'
+                                            _ -> n' + 1
+                         in encodeFloat n'' (I# (ln# -# ld# +# 1# -# md#))
+                | otherwise ->
+                  case ld# +# (me# -# md#) of
+                    ld'# | ld'# ># (ln# +# 1#)  -> encodeFloat 0 0
+                         | ld'# ==# (ln# +# 1#) ->
+                           case integerLog2IsPowerOf2# n of
+                            (# _, 0# #) -> encodeFloat 0 0
+                            (# _, _ #)  -> encodeFloat 1 (minEx - mantDigs)
+                         | ld'# <=# 0#  ->
+                           encodeFloat n (I# ((me# -# md#) -# ld'#))
+                         | otherwise    ->
+                           let n' = n `shiftR` (I# ld'#)
+                           in case roundingMode# n (ld'# -# 1#) of
+                                0# -> encodeFloat n' (minEx - mantDigs)
+                                1# -> if fromInteger n' .&. (1 :: Int) == 0
+                                        then encodeFloat n' (minEx-mantDigs)
+                                        else encodeFloat (n' + 1) (minEx-mantDigs)
+                                _  -> encodeFloat (n' + 1) (minEx-mantDigs)
+        | otherwise ->
+          let ln = I# (integerLog2# n)
+              ld = I# ld#
+              p0 = max minEx (ln - ld)
+              (n', d')
+                | p0 < mantDigs = (n `shiftL` (mantDigs - p0), d)
+                | p0 == mantDigs = (n, d)
+                | otherwise     = (n, d `shiftL` (p0 - mantDigs))
+              scale p a b
+                | p <= minEx-mantDigs = (p,a,b)
+                | a < (b `shiftL` (mantDigs-1)) = (p-1, a `shiftL` 1, b)
+                | (b `shiftL` mantDigs) <= a = (p+1, a, b `shiftL` 1)
+                | otherwise = (p, a, b)
+              (p', n'', d'') = scale (p0-mantDigs) n' d'
+              rdq = case n'' `quotRem` d'' of
+                     (q,r) -> case compare (r `shiftL` 1) d'' of
+                                LT -> q
+                                EQ -> if fromInteger q .&. (1 :: Int) == 0
+                                        then q else q+1
+                                GT -> q+1
+          in  encodeFloat rdq p'
 \end{code}
 
 
@@ -836,12 +973,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
@@ -882,12 +1013,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)