\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`
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
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
(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
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).
projects / ghc-hetmet.git / commitdiff