From 6cad3aa264eec410bb577604f512b681da1dada7 Mon Sep 17 00:00:00 2001 From: sof Date: Sun, 13 Apr 1997 02:35:07 +0000 Subject: [PATCH] [project @ 1997-04-13 02:35:07 by sof] Moved code over to Numeric --- ghc/lib/ghc/PrelNum.lhs | 189 ++--------------------------------------------- 1 file changed, 6 insertions(+), 183 deletions(-) diff --git a/ghc/lib/ghc/PrelNum.lhs b/ghc/lib/ghc/PrelNum.lhs index 3d9ff7a..2cd1078 100644 --- a/ghc/lib/ghc/PrelNum.lhs +++ b/ghc/lib/ghc/PrelNum.lhs @@ -19,13 +19,15 @@ It's rather big! \begin{code} module PrelNum where +import PrelBase +import GHC import {-# SOURCE #-} IOBase ( error ) import PrelList -import PrelBase + import ArrBase ( Array, array, (!) ) import STBase ( unsafePerformPrimIO ) import Ix ( Ix(..) ) -import GHC +import Numeric infixr 8 ^, ^^, ** infixl 7 /, %, `quot`, `rem`, `div`, `mod` @@ -397,7 +399,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 @@ -541,7 +543,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 @@ -801,144 +803,7 @@ instance (Integral a) => Show (Ratio a) where (shows x . showString " % " . shows y) \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 -early version of hbi, noticed that (read . show) was not the identity -for floating-point numbers. (There was a one-bit error about half the -time.) The original version of the conversion function was in fact -simply a floating-point divide, as you suggest above. The new version -is, I grant you, somewhat denser. - -Unfortunately, Joe's code doesn't work! Here's an example: - -main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n") - -This program prints - 0.0000000000000000 -instead of - 1.8217369128763981e-300 - -Lennart's code follows, and it works... - -\begin{pseudocode} -{-# GENERATE_SPECS fromRational__ a{Double#,Double} #-} -fromRational__ :: (RealFloat a) => Rational -> a -fromRational__ x = x' - where x' = f e - --- If the exponent of the nearest floating-point number to x --- is e, then the significand is the integer nearest xb^(-e), --- where b is the floating-point radix. We start with a good --- guess for e, and if it is correct, the exponent of the --- floating-point number we construct will again be e. If --- not, one more iteration is needed. - - f e = if e' == e then y else f e' - where y = encodeFloat (round (x * (1 % b)^^e)) e - (_,e') = decodeFloat y - b = floatRadix x' - --- We obtain a trial exponent by doing a floating-point --- division of x's numerator by its denominator. The --- result of this division may not itself be the ultimate --- result, because of an accumulation of three rounding --- errors. - - (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x' - / fromInteger (denominator x)) -\end{pseudocode} - -Now, here's Lennart's code. - \begin{code} -fromRational__ :: (RealFloat a) => Rational -> a -fromRational__ x = - if x == 0 then encodeFloat 0 0 -- Handle exceptional cases - else if x < 0 then - fromRat' (-x) -- first. - else fromRat' x - --- Conversion process: --- Scale the rational number by the RealFloat base until --- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat). --- Then round the rational to an Integer and encode it with the exponent --- that we got from the scaling. --- To speed up the scaling process we compute the log2 of the number to get --- a first guess of the exponent. - -fromRat' :: (RealFloat a) => Rational -> a -fromRat' x = r - where b = floatRadix r - p = floatDigits r - (minExp0, _) = floatRange r - minExp = minExp0 - p -- the real minimum exponent - xMin = toRational (expt b (p-1)) - xMax = toRational (expt b p) - p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp - f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1 - (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f) - r = encodeFloat (round x') p' - --- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp. -scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int) -scaleRat b minExp xMin xMax p x = - if p <= minExp then - (x, p) - else if x >= xMax then - scaleRat b minExp xMin xMax (p+1) (x/b) - else if x < xMin then - scaleRat b minExp xMin xMax (p-1) (x*b) - else - (x, p) - --- Exponentiation with a cache for the most common numbers. -minExpt = 0::Int -maxExpt = 1100::Int -expt :: Integer -> Int -> Integer -expt base n = - if base == 2 && n >= minExpt && n <= maxExpt then - expts!n - else - base^n -expts :: Array Int Integer -expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]] - --- 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. -integerLogBase :: Integer -> Integer -> Int -integerLogBase b i = - if i < b then - 0 - else - -- Try squaring the base first to cut down the number of divisions. - let l = 2 * integerLogBase (b*b) i - doDiv :: Integer -> Int -> Int - doDiv i l = if i < b then l else doDiv (i `div` b) (l+1) - 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 @@ -959,48 +824,6 @@ jtos' n cs 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 isNaN x then showString "NaN" - else if isInfinite x then - ( if x < 0 then showString "-" else id) . showString "Infinite" - else if x < 0 || isNegativeZero x then showChar '-' . showFloat' (-x) else showFloat' x - -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). -- 1.7.10.4