import {-# SOURCE #-} IOBase ( error )
import PrelList
import PrelBase
+import ArrBase ( Array, array, (!) )
+import Ix ( Ix(..) )
import GHC
infixr 8 ^, ^^, **
showsPrec x = showSignedInteger x
showList = showList__ (showsPrec 0)
+instance Ix Integer where
+ range (m,n) = [m..n]
+ index b@(m,n) i
+ | inRange b i = fromInteger (i - m)
+ | otherwise = error "Integer.index: Index out of range."
+ inRange (m,n) i = m <= i && i <= n
+
integer_0, integer_1, integer_2, integer_m1 :: Integer
-integer_0 = 0; integer_1 = 1; integer_2 = 2; integer_m1 = -1
+integer_0 = int2Integer# 0#
+integer_1 = int2Integer# 1#
+integer_2 = int2Integer# 2#
+integer_m1 = int2Integer# (negateInt# 1#)
\end{code}
(F# x) < (F# y) = x `ltFloat#` y
(F# x) <= (F# y) = x `leFloat#` y
(F# x) >= (F# y) = x `geFloat#` y
- (F# x) > (F# y) = x `geFloat#` y
+ (F# x) > (F# y) = x `gtFloat#` y
instance Num Float where
(+) x y = plusFloat x y
%*********************************************************
\begin{code}
-data (Integral a) => Ratio a = a :% a deriving (Eq)
+data (Integral a) => Ratio a = !a :% !a deriving (Eq)
type Rational = Ratio Integer
\end{code}
numerator, denominator :: (Integral a) => Ratio a -> a
approxRational :: (RealFrac a) => a -> a -> Rational
+\end{code}
+
+\tr{reduce} is a subsidiary function used only in this module .
+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 y = (x `quot` d) :% (y `quot` d)
where d = gcd x y
+\end{code}
+\begin{code}
x % y = reduce (x * signum y) (abs y)
numerator (x:%y) = x
(shows x . showString " % " . shows y)
\end{code}
-{-
-[In response to a request by simonpj, Joe Fasel writes:]
+[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.
-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:
-How's this?
+main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
-Joe
--}
+This program prints
+ 0.0000000000000000
+instead of
+ 1.8217369128763981e-300
-\begin{code}
+Lennart's code follows, and it works...
+
+\begin{pseudocode}
{-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
fromRational__ :: (RealFloat a) => Rational -> a
fromRational__ x = x'
(s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
/ fromInteger (denominator x))
-\end{code}
+\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}
%*********************************************************
%* *
projects / ghc-hetmet.git / blobdiff