-numerator (x:%y) = x
-
-denominator (x:%y) = y
-\end{code}
-
-
-@approxRational@, applied to two real fractional numbers x and epsilon,
-returns the simplest rational number within epsilon of x. A rational
-number n%d in reduced form is said to be simpler than another n'%d' if
-abs n <= abs n' && d <= d'. Any real interval contains a unique
-simplest rational; here, for simplicity, we assume a closed rational
-interval. If such an interval includes at least one whole number, then
-the simplest rational is the absolutely least whole number. Otherwise,
-the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
-and abs r' < d', and the simplest rational is q%1 + the reciprocal of
-the simplest rational between d'%r' and d%r.
-
-\begin{code}
-approxRational x eps = simplest (x-eps) (x+eps)
- where simplest x y | y < x = simplest y x
- | x == y = xr
- | x > 0 = simplest' n d n' d'
- | y < 0 = - simplest' (-n') d' (-n) d
- | otherwise = 0 :% 1
- 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
- | q /= q' = (q+1) :% 1
- | otherwise = (q*n''+d'') :% n''
- where (q,r) = quotRem n d
- (q',r') = quotRem n' d'
- nd'' = simplest' d' r' d r
- n'' = numerator nd''
- d'' = denominator nd''
-\end{code}
-
-
-\begin{code}
-instance (Integral a) => Ord (Ratio a) where
- (x:%y) <= (x':%y') = x * y' <= x' * y
- (x:%y) < (x':%y') = x * y' < x' * y
-
-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
- signum (x:%y) = signum x :% 1
- fromInteger x = fromInteger x :% 1
-
-instance (Integral a) => Real (Ratio a) where
- toRational (x:%y) = toInteger x :% toInteger y
-
-instance (Integral a) => Fractional (Ratio a) where
- (x:%y) / (x':%y') = (x*y') % (y*x')
- recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
- fromRational (x:%y) = fromInteger x :% fromInteger y
-
-instance (Integral a) => RealFrac (Ratio a) where
- properFraction (x:%y) = (fromIntegral q, r:%y)
- where (q,r) = quotRem x y
-
-instance (Integral a) => Enum (Ratio a) where
- enumFrom = iterate ((+)1)
- enumFromThen n m = iterate ((+)(m-n)) n
- toEnum n = fromIntegral n :% 1
- fromEnum = fromInteger . truncate
-
-ratio_prec :: Int
-ratio_prec = 7
-
-instance (Integral a) => Show (Ratio a) where
- showsPrec p (x:%y) = showParen (p > ratio_prec)
- (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 (toInt (n `rem` 10) + (ord_0::Int)) : 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)