[project @ 1997-10-05 21:33:42 by sof]
[ghc-hetmet.git] / ghc / lib / ghc / PrelNum.lhs
index cadad79..95724b6 100644 (file)
@@ -9,21 +9,27 @@ Numeric part of the prelude.
 It's rather big!
 
 \begin{code}
+{-# OPTIONS -fno-implicit-prelude -#include "cbits/floatExtreme.h" #-}
 {-# OPTIONS -H20m #-}
+
 #include "../includes/ieee-flpt.h"
+
 \end{code}
 
 \begin{code}
-{-# OPTIONS -fno-implicit-prelude #-}
-
 module PrelNum where
 
+import PrelBase
+import GHC
 import {-# SOURCE #-}  IOBase  ( error )
 import PrelList
-import PrelBase
+
 import ArrBase ( Array, array, (!) )
+import UnsafeST ( unsafePerformPrimIO )
 import Ix      ( Ix(..) )
-import GHC
+import Foreign ()              -- This import tells the dependency analyser to compile Foreign first.
+                               -- There's an implicit dependency on Foreign because the ccalls in
+                               -- PrelNum implicitly mention CCallable.
 
 infixr 8  ^, ^^, **
 infixl 7  /, %, `quot`, `rem`, `div`, `mod`
@@ -329,6 +335,8 @@ instance  Integral Integer where
                   -- Case-ified by WDP 94/10
 
 instance  Enum Integer  where
+    toEnum n            =  toInteger n
+    fromEnum n          =  toInt n
     enumFrom n           =  n : enumFrom (n + 1)
     enumFromThen m n     =  en' m (n - m)
                            where en' m n = m : en' (m + n) n
@@ -395,7 +403,7 @@ instance  Real Float  where
 
 instance  Fractional Float  where
     (/) x y            =  divideFloat x y
-    fromRational x     =  fromRational__ x
+    fromRational x     =  fromRat x
     recip x            =  1.0 / x
 
 instance  Floating Float  where
@@ -484,6 +492,15 @@ instance  RealFloat Float  where
 
     scaleFloat k x     = case decodeFloat x of
                            (m,n) -> encodeFloat m (n+k)
+    isNaN x = 
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
+    isInfinite x =
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatInfinite x) {- ditto! -}
+    isDenormalized x =
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatDenormalized x) -- ..
+    isNegativeZero x =
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatNegativeZero x) -- ...
+    isIEEE x    = True
 
 instance  Show Float  where
     showsPrec   x = showSigned showFloat x
@@ -530,7 +547,7 @@ instance  Real Double  where
 
 instance  Fractional Double  where
     (/) x y            =  divideDouble x y
-    fromRational x     =  fromRational__ x
+    fromRational x     =  fromRat x
     recip x            =  1.0 / x
 
 instance  Floating Double  where
@@ -627,6 +644,15 @@ instance  RealFloat Double  where
 
     scaleFloat k x     = case decodeFloat x of
                            (m,n) -> encodeFloat m (n+k)
+    isNaN x = 
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
+    isInfinite x =
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleInfinite x) {- ditto -}
+    isDenormalized x =
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleDenormalized x) -- ..
+    isNegativeZero x =
+      (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleNegativeZero x) -- ...
+    isIEEE x    = True
 
 instance  Show Double  where
     showsPrec   x = showSigned showFloat x
@@ -640,30 +666,41 @@ instance  Show Double  where
 %*                                                     *
 %*********************************************************
 
-The Enum instances for Floats and Doubles are slightly unusual.
-The `toEnum' function truncates numbers to Int.  The definitions
-of enumFrom and enumFromThen allow floats to be used in arithmetic
+The @Enum@ instances for Floats and Doubles are slightly unusual.
+The @toEnum@ function truncates numbers to Int.  The definitions
+of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
 series: [0,0.1 .. 1.0].  However, roundoff errors make these somewhat
 dubious.  This example may have either 10 or 11 elements, depending on
 how 0.1 is represented.
 
+NOTE: The instances for Float and Double do not make use of the default
+methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
+a `non-lossy' conversion to and from Ints. Instead we make use of the 
+1.2 default methods (back in the days when Enum had Ord as a superclass)
+for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
+
 \begin{code}
 instance  Enum Float  where
-    toEnum              =  fromIntegral
-    fromEnum            =  fromInteger . truncate   -- may overflow
-    enumFrom           =  numericEnumFrom
-    enumFromThen       =  numericEnumFromThen
+    toEnum         =  fromIntegral
+    fromEnum       =  fromInteger . truncate   -- may overflow
+    enumFrom      =  numericEnumFrom
+    enumFromThen   =  numericEnumFromThen
+    enumFromThenTo =  numericEnumFromThenTo
 
 instance  Enum Double  where
-    toEnum              =  fromIntegral
-    fromEnum            =  fromInteger . truncate   -- may overflow
-    enumFrom           =  numericEnumFrom
-    enumFromThen       =  numericEnumFromThen
+    toEnum         =  fromIntegral
+    fromEnum       =  fromInteger . truncate   -- may overflow
+    enumFrom      =  numericEnumFrom
+    enumFromThen   =  numericEnumFromThen
+    enumFromThenTo =  numericEnumFromThenTo
 
 numericEnumFrom                :: (Real a) => a -> [a]
 numericEnumFromThen    :: (Real a) => a -> a -> [a]
+numericEnumFromThenTo   :: (Real a) => a -> a -> a -> [a]
 numericEnumFrom                =  iterate (+1)
 numericEnumFromThen n m        =  iterate (+(m-n)) n
+numericEnumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
+                                     (numericEnumFromThen n m)
 \end{code}
 
 
@@ -690,7 +727,7 @@ It normalises a ratio by dividing both numerator and denominator by
 their greatest common divisor.
 
 \begin{code}
-reduce _ 0             =  error "{Ratio.%}: zero denominator"
+reduce x 0             =  error "{Ratio.%}: zero denominator"
 reduce x y             =  (x `quot` d) :% (y `quot` d)
                           where d = gcd x y
 \end{code}
@@ -722,8 +759,12 @@ approxRational x eps       =  simplest (x-eps) (x+eps)
                           | x > 0      =  simplest' n d n' d'
                           | y < 0      =  - simplest' (-n') d' (-n) d
                           | otherwise  =  0 :% 1
-                                       where xr@(n:%d) = toRational x
-                                             (n':%d')  = toRational y
+                                       where xr  = toRational x
+                                             n   = numerator xr
+                                             d   = denominator xr
+                                             nd' = toRational y
+                                             n'  = numerator nd'
+                                             d'  = denominator nd'
 
              simplest' n d n' d'       -- assumes 0 < n%d < n'%d'
                        | r == 0     =  q :% 1
@@ -731,7 +772,9 @@ approxRational x eps        =  simplest (x-eps) (x+eps)
                        | otherwise  =  (q*n''+d'') :% n''
                                     where (q,r)      =  quotRem n d
                                           (q',r')    =  quotRem n' d'
-                                          (n'':%d'') =  simplest' d' r' d r
+                                          nd''       =  simplest' d' r' d r
+                                          n''        =  numerator nd''
+                                          d''        =  denominator nd''
 \end{code}
 
 
@@ -742,6 +785,7 @@ instance  (Integral a)      => Ord (Ratio a)  where
 
 instance  (Integral a) => Num (Ratio a)  where
     (x:%y) + (x':%y')  =  reduce (x*y' + x'*y) (y*y')
+    (x:%y) - (x':%y')  =  reduce (x*y' - x'*y) (y*y')
     (x:%y) * (x':%y')  =  reduce (x * x') (y * y')
     negate (x:%y)      =  (-x) :% y
     abs (x:%y)         =  abs x :% y
@@ -774,6 +818,252 @@ instance  (Integral a)  => Show (Ratio a)  where
                               (shows x . showString " % " . shows y)
 \end{code}
 
+\begin{code}
+--Exported from std library Numeric, defined here to
+--avoid mut. rec. between PrelNum and Numeric.
+showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
+showSigned showPos p x = if x < 0 then showParen (p > 6)
+                                                (showChar '-' . showPos (-x))
+                                 else showPos x
+
+showSignedInteger :: Int -> Integer -> ShowS
+showSignedInteger p n r
+  = -- from HBC version; support code follows
+    if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
+
+jtos :: Integer -> String
+jtos n 
+  = if n < 0 then
+        '-' : jtos' (-n) []
+    else 
+       jtos' n []
+
+jtos' :: Integer -> String -> String
+jtos' n cs
+  = if n < 10 then
+       chr (fromInteger (n + ord_0)) : cs
+    else 
+       jtos' (n `quot` 10) (chr (fromInteger (n `rem` 10 + ord_0)) : cs)
+
+showFloat x  =  showString (formatRealFloat FFGeneric Nothing x)
+
+-- These are the format types.  This type is not exported.
+
+data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
+
+formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
+formatRealFloat fmt decs x = s
+ where 
+  base = 10
+  s = if isNaN x 
+      then "NaN"
+      else 
+       if isInfinite x then
+          if x < 0 then "-Infinity" else "Infinity"
+       else
+          if x < 0 || isNegativeZero x then
+            '-':doFmt fmt (floatToDigits (toInteger base) (-x))
+          else
+           doFmt fmt (floatToDigits (toInteger base) x)
+
+  doFmt fmt (is, e) =
+    let ds = map intToDigit is in
+    case fmt of
+     FFGeneric ->
+      doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
+           (is,e)
+     FFExponent ->
+      case decs of
+       Nothing ->
+        let e' = if e==0 then 0 else e-1 in
+       (case ds of
+          [d]    -> d : ".0e"
+         (d:ds) -> d : '.' : ds ++ "e") ++ show e'
+       Just dec ->
+        let dec' = max dec 1 in
+        case is of
+         [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
+         _ ->
+          let
+          (ei,is') = roundTo base (dec'+1) is
+          d:ds = map intToDigit (if ei > 0 then init is' else is')
+          in
+         d:'.':ds ++ 'e':show (e-1+ei)
+     FFFixed ->
+      let
+       mk0 ls = case ls of { "" -> "0" ; _ -> ls}
+      in
+      case decs of
+       Nothing ->
+         let
+         f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
+         f n s "" = f (n-1) ('0':s) ""
+         f n s (d:ds) = f (n-1) (d:s) ds
+        in
+        f e "" ds
+       Just dec ->
+        let dec' = max dec 1 in
+       if e >= 0 then
+        let
+         (ei,is') = roundTo base (dec' + e) is
+         (ls,rs)  = splitAt (e+ei) (map intToDigit is')
+        in
+        mk0 ls ++ (if null rs then "" else '.':rs)
+       else
+        let
+         (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
+         d:ds = map intToDigit (if ei > 0 then is' else 0:is')
+        in
+        d : '.' : ds
+        
+
+roundTo :: Int -> Int -> [Int] -> (Int,[Int])
+roundTo base d is =
+ let
+  v = f d is
+ in
+ case v of
+  (0,is) -> v
+  (1,is) -> (1, 1:is)
+ where
+  b2 = base `div` 2
+
+  f n [] = (0, replicate n 0)
+  f 0 (i:_) = (if i>=b2 then 1 else 0, [])
+  f d (i:is) =
+    let 
+     (c,ds) = f (d-1) is
+     i' = c + i
+    in
+    if i' == base then (1,0:ds) else (0,i':ds)
+
+--
+-- Based on "Printing Floating-Point Numbers Quickly and Accurately"
+-- by R.G. Burger and R.K. Dybvig in PLDI 96.
+-- This version uses a much slower logarithm estimator. It should be improved.
+
+-- This function returns a list of digits (Ints in [0..base-1]) and an
+-- exponent.
+--floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
+floatToDigits _ 0 = ([0], 0)
+floatToDigits base x =
+ let 
+  (f0, e0) = decodeFloat x
+  (minExp0, _) = floatRange x
+  p = floatDigits x
+  b = floatRadix x
+  minExp = minExp0 - p -- the real minimum exponent
+  -- Haskell requires that f be adjusted so denormalized numbers
+  -- will have an impossibly low exponent.  Adjust for this.
+  (f, e) = 
+   let n = minExp - e0 in
+   if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
+  (r, s, mUp, mDn) =
+   if e >= 0 then
+    let be = b^ e in
+    if f == b^(p-1) then
+      (f*be*b*2, 2*b, be*b, b)
+    else
+      (f*be*2, 2, be, be)
+   else
+    if e > minExp && f == b^(p-1) then
+      (f*b*2, b^(-e+1)*2, b, 1)
+    else
+      (f*2, b^(-e)*2, 1, 1)
+  k =
+   let 
+    k0 =
+     if b == 2 && base == 10 then
+        -- logBase 10 2 is slightly bigger than 3/10 so
+       -- the following will err on the low side.  Ignoring
+       -- the fraction will make it err even more.
+       -- Haskell promises that p-1 <= logBase b f < p.
+       (p - 1 + e0) * 3 `div` 10
+     else
+        ceiling ((log (fromInteger (f+1)) +
+                fromInt e * log (fromInteger b)) /
+                 fromInt e * log (fromInteger b))
+
+    fixup n =
+      if n >= 0 then
+        if r + mUp <= expt base n * s then n else fixup (n+1)
+      else
+        if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
+   in
+   fixup k0
+
+  gen ds rn sN mUpN mDnN =
+   let
+    (dn, rn') = (rn * base) `divMod` sN
+    mUpN' = mUpN * base
+    mDnN' = mDnN * base
+   in
+   case (rn' < mDnN', rn' + mUpN' > sN) of
+    (True,  False) -> dn : ds
+    (False, True)  -> dn+1 : ds
+    (True,  True)  -> if rn' * 2 < sN then dn : ds else dn+1 : ds
+    (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
+  
+  rds = 
+   if k >= 0 then
+      gen [] r (s * expt base k) mUp mDn
+   else
+     let bk = expt base (-k) in
+     gen [] (r * bk) s (mUp * bk) (mDn * bk)
+ in
+ (map toInt (reverse rds), k)
+
+\end{code}
+
+@showRational@ converts a Rational to a string that looks like a
+floating point number, but without converting to any floating type
+(because of the possible overflow).
+
+From/by Lennart, 94/09/26
+
+\begin{code}
+showRational :: Int -> Rational -> String
+showRational n r =
+    if r == 0 then
+       "0.0"
+    else
+       let (r', e) = normalize r
+       in  prR n r' e
+
+startExpExp = 4 :: Int
+
+-- make sure 1 <= r < 10
+normalize :: Rational -> (Rational, Int)
+normalize r = if r < 1 then
+                 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
+             else
+                 norm startExpExp r 0
+       where norm :: Int -> Rational -> Int -> (Rational, Int)
+             -- Invariant: r*10^e == original r
+             norm 0  r e = (r, e)
+             norm ee r e =
+               let n = 10^ee
+                   tn = 10^n
+               in  if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
+
+drop0 "" = ""
+drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
+
+prR :: Int -> Rational -> Int -> String
+prR n r e | r <  1  = prR n (r*10) (e-1)               -- final adjustment
+prR n r e | r >= 10 = prR n (r/10) (e+1)
+prR n r e0 =
+       let s = show ((round (r * 10^n))::Integer)
+           e = e0+1
+       in  if e > 0 && e < 8 then
+               take e s ++ "." ++ drop0 (drop e s)
+           else if e <= 0 && e > -3 then
+               "0." ++ take (-e) (repeat '0') ++ drop0 s
+           else
+               head s : "."++ drop0 (tail s) ++ "e" ++ show e0
+\end{code}
+
+
 [In response to a request for documentation of how fromRational works,
 Joe Fasel writes:] A quite reasonable request!  This code was added to
 the Prelude just before the 1.2 release, when Lennart, working with an
@@ -796,8 +1086,8 @@ Lennart's code follows, and it works...
 
 \begin{pseudocode}
 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
-fromRational__ :: (RealFloat a) => Rational -> a
-fromRational__ x = x'
+fromRat :: (RealFloat a) => Rational -> a
+fromRat x = x'
        where x' = f e
 
 --             If the exponent of the nearest floating-point number to x 
@@ -825,8 +1115,8 @@ fromRational__ x = x'
 Now, here's Lennart's code.
 
 \begin{code}
-fromRational__ :: (RealFloat a) => Rational -> a
-fromRational__ x = 
+--fromRat :: (RealFloat a) => Rational -> a
+fromRat x = 
     if x == 0 then encodeFloat 0 0             -- Handle exceptional cases
     else if x < 0 then - fromRat' (-x)         -- first.
     else fromRat' x
@@ -891,129 +1181,6 @@ integerLogBase b i =
        in  doDiv (i `div` (b^l)) l
 \end{code}
 
-%*********************************************************
-%*                                                     *
-\subsection{Showing numbers}
-%*                                                     *
-%*********************************************************
-
-\begin{code}
-showInteger n r
-  = case quotRem n 10 of                    { (n', d) ->
-    case (chr (ord_0 + fromIntegral d)) of { C# c# -> -- stricter than necessary
-    let
-       r' = C# c# : r
-    in
-    if n' == 0 then r' else showInteger n' r'
-    }}
-
-showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
-showSigned showPos p x = if x < 0 then showParen (p > 6)
-                                                (showChar '-' . showPos (-x))
-                                 else showPos x
-
-showSignedInteger :: Int -> Integer -> ShowS
-showSignedInteger p n r
-  = -- from HBC version; support code follows
-    if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
-
-jtos :: Integer -> String
-jtos n 
-  = if n < 0 then
-        '-' : jtos' (-n) []
-    else 
-       jtos' n []
-
-jtos' :: Integer -> String -> String
-jtos' n cs
-  = if n < 10 then
-       chr (fromInteger (n + ord_0)) : cs
-    else 
-       jtos' (n `quot` 10) (chr (fromInteger (n `rem` 10 + ord_0)) : cs)
-\end{code}
-
-The functions showFloat below uses rational arithmetic
-to insure correct conversion between the floating-point radix and
-decimal.  It is often possible to use a higher-precision floating-
-point type to obtain the same results.
-
-\begin{code}
-{-# GENERATE_SPECS showFloat a{Double#,Double} #-}
-showFloat:: (RealFloat a) => a -> ShowS
-showFloat x =
-    if x == 0 then showString ("0." ++ take (m-1) zeros)
-             else if e >= m-1 || e < 0 then showSci else showFix
-    where
-    showFix    = showString whole . showChar '.' . showString frac
-                 where (whole,frac) = splitAt (e+1) (show sig)
-    showSci    = showChar d . showChar '.' . showString frac
-                     . showChar 'e' . shows e
-                 where (d:frac) = show sig
-    (m, sig, e) = if b == 10 then (w,          s,   n+w-1)
-                            else (m', sig', e'   )
-    m'         = ceiling
-                     ((fromInt w * log (fromInteger b)) / log 10 :: Double)
-                 + 1
-    (sig', e') = if      sig1 >= 10^m'     then (round (t/10), e1+1)
-                 else if sig1 <  10^(m'-1) then (round (t*10), e1-1)
-                                           else (sig1,          e1  )
-    sig1       = round t
-    t          = s%1 * (b%1)^^n * 10^^(m'-e1-1)
-    e1         = floor (logBase 10 x)
-    (s, n)     = decodeFloat x
-    b          = floatRadix x
-    w          = floatDigits x
-zeros = repeat '0'
-\end{code}
-
-@showRational@ converts a Rational to a string that looks like a
-floating point number, but without converting to any floating type
-(because of the possible overflow).
-
-From/by Lennart, 94/09/26
-
-\begin{code}
-showRational :: Int -> Rational -> String
-showRational n r =
-    if r == 0 then
-       "0.0"
-    else
-       let (r', e) = normalize r
-       in  prR n r' e
-
-startExpExp = 4 :: Int
-
--- make sure 1 <= r < 10
-normalize :: Rational -> (Rational, Int)
-normalize r = if r < 1 then
-                 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
-             else
-                 norm startExpExp r 0
-       where norm :: Int -> Rational -> Int -> (Rational, Int)
-             -- Invariant: r*10^e == original r
-             norm 0  r e = (r, e)
-             norm ee r e =
-               let n = 10^ee
-                   tn = 10^n
-               in  if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
-
-drop0 "" = ""
-drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
-
-prR :: Int -> Rational -> Int -> String
-prR n r e | r <  1  = prR n (r*10) (e-1)               -- final adjustment
-prR n r e | r >= 10 = prR n (r/10) (e+1)
-prR n r e0 =
-       let s = show ((round (r * 10^n))::Integer)
-           e = e0+1
-       in  if e > 0 && e < 8 then
-               take e s ++ "." ++ drop0 (drop e s)
-           else if e <= 0 && e > -3 then
-               "0." ++ take (-e) (repeat '0') ++ drop0 s
-           else
-               head s : "."++ drop0 (tail s) ++ "e" ++ show e0
-\end{code}
 
 %*********************************************************
 %*                                                     *