2 % (c) The AQUA Project, Glasgow University, 1994-1996
5 \section[PrelNum]{Module @PrelNum@}
7 Numeric part of the prelude.
12 {-# OPTIONS -fno-implicit-prelude -#include "cbits/floatExtreme.h" #-}
15 #include "../includes/ieee-flpt.h"
24 import {-# SOURCE #-} PrelErr ( error )
28 import PrelArr ( Array, array, (!) )
29 import PrelIOBase ( unsafePerformIO )
31 import PrelCCall () -- we need the definitions of CCallable and
32 -- CReturnable for the _ccall_s herein.
36 infixl 7 /, %, `quot`, `rem`, `div`, `mod`
40 %*********************************************************
42 \subsection{Standard numeric classes}
44 %*********************************************************
47 class (Num a, Ord a) => Real a where
48 toRational :: a -> Rational
50 class (Real a, Enum a) => Integral a where
51 quot, rem, div, mod :: a -> a -> a
52 quotRem, divMod :: a -> a -> (a,a)
53 toInteger :: a -> Integer
54 toInt :: a -> Int -- partain: Glasgow extension
56 n `quot` d = q where (q,r) = quotRem n d
57 n `rem` d = r where (q,r) = quotRem n d
58 n `div` d = q where (q,r) = divMod n d
59 n `mod` d = r where (q,r) = divMod n d
60 divMod n d = if signum r == negate (signum d) then (q-1, r+d) else qr
61 where qr@(q,r) = quotRem n d
63 class (Num a) => Fractional a where
66 fromRational :: Rational -> a
70 class (Fractional a) => Floating a where
72 exp, log, sqrt :: a -> a
73 (**), logBase :: a -> a -> a
74 sin, cos, tan :: a -> a
75 asin, acos, atan :: a -> a
76 sinh, cosh, tanh :: a -> a
77 asinh, acosh, atanh :: a -> a
79 x ** y = exp (log x * y)
80 logBase x y = log y / log x
83 tanh x = sinh x / cosh x
85 class (Real a, Fractional a) => RealFrac a where
86 properFraction :: (Integral b) => a -> (b,a)
87 truncate, round :: (Integral b) => a -> b
88 ceiling, floor :: (Integral b) => a -> b
90 truncate x = m where (m,_) = properFraction x
92 round x = let (n,r) = properFraction x
93 m = if r < 0 then n - 1 else n + 1
94 in case signum (abs r - 0.5) of
96 0 -> if even n then n else m
99 ceiling x = if r > 0 then n + 1 else n
100 where (n,r) = properFraction x
102 floor x = if r < 0 then n - 1 else n
103 where (n,r) = properFraction x
105 class (RealFrac a, Floating a) => RealFloat a where
106 floatRadix :: a -> Integer
107 floatDigits :: a -> Int
108 floatRange :: a -> (Int,Int)
109 decodeFloat :: a -> (Integer,Int)
110 encodeFloat :: Integer -> Int -> a
112 significand :: a -> a
113 scaleFloat :: Int -> a -> a
114 isNaN, isInfinite, isDenormalized, isNegativeZero, isIEEE
117 exponent x = if m == 0 then 0 else n + floatDigits x
118 where (m,n) = decodeFloat x
120 significand x = encodeFloat m (negate (floatDigits x))
121 where (m,_) = decodeFloat x
123 scaleFloat k x = encodeFloat m (n+k)
124 where (m,n) = decodeFloat x
127 %*********************************************************
129 \subsection{Overloaded numeric functions}
131 %*********************************************************
134 even, odd :: (Integral a) => a -> Bool
135 even n = n `rem` 2 == 0
138 {-# SPECIALISE gcd ::
140 Integer -> Integer -> Integer #-}
141 gcd :: (Integral a) => a -> a -> a
142 gcd 0 0 = error "Prelude.gcd: gcd 0 0 is undefined"
143 gcd x y = gcd' (abs x) (abs y)
145 gcd' x y = gcd' y (x `rem` y)
147 {-# SPECIALISE lcm ::
149 Integer -> Integer -> Integer #-}
150 lcm :: (Integral a) => a -> a -> a
153 lcm x y = abs ((x `quot` (gcd x y)) * y)
155 {-# SPECIALISE (^) ::
156 Integer -> Integer -> Integer,
157 Integer -> Int -> Integer,
158 Int -> Int -> Int #-}
159 (^) :: (Num a, Integral b) => a -> b -> a
161 x ^ n | n > 0 = f x (n-1) x
163 f x n y = g x n where
164 g x n | even n = g (x*x) (n `quot` 2)
165 | otherwise = f x (n-1) (x*y)
166 _ ^ _ = error "Prelude.^: negative exponent"
168 {-# SPECIALISE (^^) ::
169 Double -> Int -> Double,
170 Rational -> Int -> Rational #-}
171 (^^) :: (Fractional a, Integral b) => a -> b -> a
172 x ^^ n = if n >= 0 then x^n else recip (x^(negate n))
174 {-# SPECIALIZE fromIntegral ::
184 Integer -> Double #-}
185 fromIntegral :: (Integral a, Num b) => a -> b
186 fromIntegral = fromInteger . toInteger
188 {-# SPECIALIZE fromRealFrac ::
193 Rational -> Rational,
198 fromRealFrac :: (RealFrac a, Fractional b) => a -> b
199 fromRealFrac = fromRational . toRational
201 atan2 :: (RealFloat a) => a -> a -> a
202 atan2 y x = case (signum y, signum x) of
206 (-1, 0) -> (negate pi)/2
207 ( _, 1) -> atan (y/x)
208 ( _,-1) -> atan (y/x) + pi
209 ( 0, 0) -> error "Prelude.atan2: atan2 of origin"
213 %*********************************************************
215 \subsection{Instances for @Int@}
217 %*********************************************************
220 instance Real Int where
221 toRational x = toInteger x % 1
223 instance Integral Int where
224 a@(I# _) `quotRem` b@(I# _) = (a `quotInt` b, a `remInt` b)
225 -- OK, so I made it a little stricter. Shoot me. (WDP 94/10)
227 -- Following chks for zero divisor are non-standard (WDP)
228 a `quot` b = if b /= 0
230 else error "Integral.Int.quot{PreludeCore}: divide by 0\n"
231 a `rem` b = if b /= 0
233 else error "Integral.Int.rem{PreludeCore}: divide by 0\n"
236 | n > 0 && d < 0 = mk_neg (quotInt (n-d-1) d)
237 | n < 0 && d > 0 = mk_neg (quotInt (n-d+1) d)
238 | otherwise = quotInt n d
241 - the result of (integral) division is
242 defined as being truncated towards
243 negative infinity. (see Sec 6.3.2 of
244 the Haskell 1.4 report.)
246 - in the case of Int, if either nominator or
247 denominator is negative, we adjust the nominator
248 to account for the above property before
249 computing the quotient.
251 - in the case of Int, the adjustment of the
252 nominator runs the risk of overflowing. If
253 we make the assumption that arithmetic is
254 modulo word size, and adjust the final result
262 x `mod` y = if x > 0 && y < 0 || x < 0 && y > 0 then
263 if r/=0 then r+y else 0
268 divMod x@(I# _) y@(I# _) = (x `div` y, x `mod` y)
269 -- Stricter. Sorry if you don't like it. (WDP 94/10)
271 --OLD: even x = eqInt (x `mod` 2) 0
272 --OLD: odd x = neInt (x `mod` 2) 0
274 toInteger (I# n#) = int2Integer# n# -- give back a full-blown Integer
280 %*********************************************************
282 \subsection{Type @Integer@}
284 %*********************************************************
286 These types are used to return from integer primops
289 data Return2GMPs = Return2GMPs Int# Int# ByteArray# Int# Int# ByteArray#
290 data ReturnIntAndGMP = ReturnIntAndGMP Int# Int# Int# ByteArray#
296 instance Eq Integer where
297 (J# a1 s1 d1) == (J# a2 s2 d2)
298 = (cmpInteger# a1 s1 d1 a2 s2 d2) ==# 0#
300 (J# a1 s1 d1) /= (J# a2 s2 d2)
301 = (cmpInteger# a1 s1 d1 a2 s2 d2) /=# 0#
303 instance Ord Integer where
304 (J# a1 s1 d1) <= (J# a2 s2 d2)
305 = (cmpInteger# a1 s1 d1 a2 s2 d2) <=# 0#
307 (J# a1 s1 d1) < (J# a2 s2 d2)
308 = (cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#
310 (J# a1 s1 d1) >= (J# a2 s2 d2)
311 = (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
313 (J# a1 s1 d1) > (J# a2 s2 d2)
314 = (cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#
316 x@(J# a1 s1 d1) `max` y@(J# a2 s2 d2)
317 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#) then x else y
319 x@(J# a1 s1 d1) `min` y@(J# a2 s2 d2)
320 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#) then x else y
322 compare (J# a1 s1 d1) (J# a2 s2 d2)
323 = case cmpInteger# a1 s1 d1 a2 s2 d2 of { res# ->
324 if res# <# 0# then LT else
325 if res# ># 0# then GT else EQ
328 instance Num Integer where
329 (+) (J# a1 s1 d1) (J# a2 s2 d2)
330 = plusInteger# a1 s1 d1 a2 s2 d2
332 (-) (J# a1 s1 d1) (J# a2 s2 d2)
333 = minusInteger# a1 s1 d1 a2 s2 d2
335 negate (J# a s d) = negateInteger# a s d
337 (*) (J# a1 s1 d1) (J# a2 s2 d2)
338 = timesInteger# a1 s1 d1 a2 s2 d2
340 -- ORIG: abs n = if n >= 0 then n else -n
343 = case 0 of { J# a2 s2 d2 ->
344 if (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
346 else negateInteger# a1 s1 d1
349 signum n@(J# a1 s1 d1)
350 = case 0 of { J# a2 s2 d2 ->
352 cmp = cmpInteger# a1 s1 d1 a2 s2 d2
355 else if cmp ==# 0# then 0
361 fromInt (I# n#) = int2Integer# n# -- gives back a full-blown Integer
363 instance Real Integer where
366 instance Integral Integer where
367 quotRem (J# a1 s1 d1) (J# a2 s2 d2)
368 = case (quotRemInteger# a1 s1 d1 a2 s2 d2) of
369 Return2GMPs a3 s3 d3 a4 s4 d4
370 -> (J# a3 s3 d3, J# a4 s4 d4)
372 {- USING THE UNDERLYING "GMP" CODE IS DUBIOUS FOR NOW:
374 divMod (J# a1 s1 d1) (J# a2 s2 d2)
375 = case (divModInteger# a1 s1 d1 a2 s2 d2) of
376 Return2GMPs a3 s3 d3 a4 s4 d4
377 -> (J# a3 s3 d3, J# a4 s4 d4)
380 toInt (J# a s d) = case (integer2Int# a s d) of { n# -> I# n# }
382 -- the rest are identical to the report default methods;
383 -- you get slightly better code if you let the compiler
384 -- see them right here:
385 n `quot` d = if d /= 0 then q else
386 error "Integral.Integer.quot{PreludeCore}: divide by 0\n"
387 where (q,r) = quotRem n d
388 n `rem` d = if d /= 0 then r else
389 error "Integral.Integer.quot{PreludeCore}: divide by 0\n"
390 where (q,r) = quotRem n d
391 n `div` d = q where (q,r) = divMod n d
392 n `mod` d = r where (q,r) = divMod n d
394 divMod n d = case (quotRem n d) of { qr@(q,r) ->
395 if signum r == negate (signum d) then (q - 1, r+d) else qr }
396 -- Case-ified by WDP 94/10
398 instance Enum Integer where
399 toEnum n = toInteger n
401 enumFrom n = n : enumFrom (n + 1)
402 enumFromThen m n = en' m (n - m)
403 where en' m n = m : en' (m + n) n
404 enumFromTo n m = takeWhile (<= m) (enumFrom n)
405 enumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
408 instance Show Integer where
409 showsPrec x = showSignedInteger x
410 showList = showList__ (showsPrec 0)
412 instance Ix Integer where
415 | inRange b i = fromInteger (i - m)
416 | otherwise = error "Integer.index: Index out of range."
417 inRange (m,n) i = m <= i && i <= n
419 integer_0, integer_1, integer_2, integer_m1 :: Integer
420 integer_0 = int2Integer# 0#
421 integer_1 = int2Integer# 1#
422 integer_2 = int2Integer# 2#
423 integer_m1 = int2Integer# (negateInt# 1#)
427 %*********************************************************
429 \subsection{Type @Float@}
431 %*********************************************************
434 instance Eq Float where
435 (F# x) == (F# y) = x `eqFloat#` y
437 instance Ord Float where
438 (F# x) `compare` (F# y) | x `ltFloat#` y = LT
439 | x `eqFloat#` y = EQ
442 (F# x) < (F# y) = x `ltFloat#` y
443 (F# x) <= (F# y) = x `leFloat#` y
444 (F# x) >= (F# y) = x `geFloat#` y
445 (F# x) > (F# y) = x `gtFloat#` y
447 instance Num Float where
448 (+) x y = plusFloat x y
449 (-) x y = minusFloat x y
450 negate x = negateFloat x
451 (*) x y = timesFloat x y
453 | otherwise = negateFloat x
454 signum x | x == 0.0 = 0
456 | otherwise = negate 1
457 fromInteger n = encodeFloat n 0
458 fromInt i = int2Float i
460 instance Real Float where
461 toRational x = (m%1)*(b%1)^^n
462 where (m,n) = decodeFloat x
465 instance Fractional Float where
466 (/) x y = divideFloat x y
467 fromRational x = fromRat x
470 instance Floating Float where
471 pi = 3.141592653589793238
484 (**) x y = powerFloat x y
485 logBase x y = log y / log x
487 asinh x = log (x + sqrt (1.0+x*x))
488 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
489 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
491 instance RealFrac Float where
493 {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
494 {-# SPECIALIZE truncate :: Float -> Int #-}
495 {-# SPECIALIZE round :: Float -> Int #-}
496 {-# SPECIALIZE ceiling :: Float -> Int #-}
497 {-# SPECIALIZE floor :: Float -> Int #-}
499 {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
500 {-# SPECIALIZE truncate :: Float -> Integer #-}
501 {-# SPECIALIZE round :: Float -> Integer #-}
502 {-# SPECIALIZE ceiling :: Float -> Integer #-}
503 {-# SPECIALIZE floor :: Float -> Integer #-}
506 = case (decodeFloat x) of { (m,n) ->
507 let b = floatRadix x in
509 (fromInteger m * fromInteger b ^ n, 0.0)
511 case (quotRem m (b^(negate n))) of { (w,r) ->
512 (fromInteger w, encodeFloat r n)
516 truncate x = case properFraction x of
519 round x = case properFraction x of
521 m = if r < 0.0 then n - 1 else n + 1
522 half_down = abs r - 0.5
524 case (compare half_down 0.0) of
526 EQ -> if even n then n else m
529 ceiling x = case properFraction x of
530 (n,r) -> if r > 0.0 then n + 1 else n
532 floor x = case properFraction x of
533 (n,r) -> if r < 0.0 then n - 1 else n
535 instance RealFloat Float where
536 floatRadix _ = FLT_RADIX -- from float.h
537 floatDigits _ = FLT_MANT_DIG -- ditto
538 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
541 = case decodeFloat# f# of
542 ReturnIntAndGMP exp# a# s# d# ->
543 (J# a# s# d#, I# exp#)
545 encodeFloat (J# a# s# d#) (I# e#)
546 = case encodeFloat# a# s# d# e# of { flt# -> F# flt# }
548 exponent x = case decodeFloat x of
549 (m,n) -> if m == 0 then 0 else n + floatDigits x
551 significand x = case decodeFloat x of
552 (m,_) -> encodeFloat m (negate (floatDigits x))
554 scaleFloat k x = case decodeFloat x of
555 (m,n) -> encodeFloat m (n+k)
557 (0::Int) /= unsafePerformIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
559 (0::Int) /= unsafePerformIO (_ccall_ isFloatInfinite x) {- ditto! -}
561 (0::Int) /= unsafePerformIO (_ccall_ isFloatDenormalized x) -- ..
563 (0::Int) /= unsafePerformIO (_ccall_ isFloatNegativeZero x) -- ...
566 instance Show Float where
567 showsPrec x = showSigned showFloat x
568 showList = showList__ (showsPrec 0)
571 %*********************************************************
573 \subsection{Type @Double@}
575 %*********************************************************
578 instance Eq Double where
579 (D# x) == (D# y) = x ==## y
581 instance Ord Double where
582 (D# x) `compare` (D# y) | x <## y = LT
586 (D# x) < (D# y) = x <## y
587 (D# x) <= (D# y) = x <=## y
588 (D# x) >= (D# y) = x >=## y
589 (D# x) > (D# y) = x >## y
591 instance Num Double where
592 (+) x y = plusDouble x y
593 (-) x y = minusDouble x y
594 negate x = negateDouble x
595 (*) x y = timesDouble x y
597 | otherwise = negateDouble x
598 signum x | x == 0.0 = 0
600 | otherwise = negate 1
601 fromInteger n = encodeFloat n 0
602 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
604 instance Real Double where
605 toRational x = (m%1)*(b%1)^^n
606 where (m,n) = decodeFloat x
609 instance Fractional Double where
610 (/) x y = divideDouble x y
611 fromRational x = fromRat x
614 instance Floating Double where
615 pi = 3.141592653589793238
618 sqrt x = sqrtDouble x
622 asin x = asinDouble x
623 acos x = acosDouble x
624 atan x = atanDouble x
625 sinh x = sinhDouble x
626 cosh x = coshDouble x
627 tanh x = tanhDouble x
628 (**) x y = powerDouble x y
629 logBase x y = log y / log x
631 asinh x = log (x + sqrt (1.0+x*x))
632 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
633 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
635 instance RealFrac Double where
637 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
638 {-# SPECIALIZE truncate :: Double -> Int #-}
639 {-# SPECIALIZE round :: Double -> Int #-}
640 {-# SPECIALIZE ceiling :: Double -> Int #-}
641 {-# SPECIALIZE floor :: Double -> Int #-}
643 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
644 {-# SPECIALIZE truncate :: Double -> Integer #-}
645 {-# SPECIALIZE round :: Double -> Integer #-}
646 {-# SPECIALIZE ceiling :: Double -> Integer #-}
647 {-# SPECIALIZE floor :: Double -> Integer #-}
649 #if defined(__UNBOXED_INSTANCES__)
650 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
651 {-# SPECIALIZE truncate :: Double -> Int# #-}
652 {-# SPECIALIZE round :: Double -> Int# #-}
653 {-# SPECIALIZE ceiling :: Double -> Int# #-}
654 {-# SPECIALIZE floor :: Double -> Int# #-}
658 = case (decodeFloat x) of { (m,n) ->
659 let b = floatRadix x in
661 (fromInteger m * fromInteger b ^ n, 0.0)
663 case (quotRem m (b^(negate n))) of { (w,r) ->
664 (fromInteger w, encodeFloat r n)
668 truncate x = case properFraction x of
671 round x = case properFraction x of
673 m = if r < 0.0 then n - 1 else n + 1
674 half_down = abs r - 0.5
676 case (compare half_down 0.0) of
678 EQ -> if even n then n else m
681 ceiling x = case properFraction x of
682 (n,r) -> if r > 0.0 then n + 1 else n
684 floor x = case properFraction x of
685 (n,r) -> if r < 0.0 then n - 1 else n
687 instance RealFloat Double where
688 floatRadix _ = FLT_RADIX -- from float.h
689 floatDigits _ = DBL_MANT_DIG -- ditto
690 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
693 = case decodeDouble# d# of
694 ReturnIntAndGMP exp# a# s# d# ->
695 (J# a# s# d#, I# exp#)
697 encodeFloat (J# a# s# d#) (I# e#)
698 = case encodeDouble# a# s# d# e# of { dbl# -> D# dbl# }
700 exponent x = case decodeFloat x of
701 (m,n) -> if m == 0 then 0 else n + floatDigits x
703 significand x = case decodeFloat x of
704 (m,_) -> encodeFloat m (negate (floatDigits x))
706 scaleFloat k x = case decodeFloat x of
707 (m,n) -> encodeFloat m (n+k)
709 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
711 (0::Int) /= unsafePerformIO (_ccall_ isDoubleInfinite x) {- ditto -}
713 (0::Int) /= unsafePerformIO (_ccall_ isDoubleDenormalized x) -- ..
715 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNegativeZero x) -- ...
718 instance Show Double where
719 showsPrec x = showSigned showFloat x
720 showList = showList__ (showsPrec 0)
724 %*********************************************************
726 \subsection{Common code for @Float@ and @Double@}
728 %*********************************************************
730 The @Enum@ instances for Floats and Doubles are slightly unusual.
731 The @toEnum@ function truncates numbers to Int. The definitions
732 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
733 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
734 dubious. This example may have either 10 or 11 elements, depending on
735 how 0.1 is represented.
737 NOTE: The instances for Float and Double do not make use of the default
738 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
739 a `non-lossy' conversion to and from Ints. Instead we make use of the
740 1.2 default methods (back in the days when Enum had Ord as a superclass)
741 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
744 instance Enum Float where
745 toEnum = fromIntegral
746 fromEnum = fromInteger . truncate -- may overflow
747 enumFrom = numericEnumFrom
748 enumFromThen = numericEnumFromThen
749 enumFromThenTo = numericEnumFromThenTo
751 instance Enum Double where
752 toEnum = fromIntegral
753 fromEnum = fromInteger . truncate -- may overflow
754 enumFrom = numericEnumFrom
755 enumFromThen = numericEnumFromThen
756 enumFromThenTo = numericEnumFromThenTo
758 numericEnumFrom :: (Real a) => a -> [a]
759 numericEnumFromThen :: (Real a) => a -> a -> [a]
760 numericEnumFromThenTo :: (Real a) => a -> a -> a -> [a]
761 numericEnumFrom = iterate (+1)
762 numericEnumFromThen n m = iterate (+(m-n)) n
763 numericEnumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
764 (numericEnumFromThen n m)
768 %*********************************************************
770 \subsection{The @Ratio@ and @Rational@ types}
772 %*********************************************************
775 data (Eval a, Integral a) => Ratio a = !a :% !a deriving (Eq)
776 type Rational = Ratio Integer
780 {-# SPECIALISE (%) :: Integer -> Integer -> Rational #-}
782 (%) :: (Integral a) => a -> a -> Ratio a
783 numerator, denominator :: (Integral a) => Ratio a -> a
784 approxRational :: (RealFrac a) => a -> a -> Rational
788 \tr{reduce} is a subsidiary function used only in this module .
789 It normalises a ratio by dividing both numerator and denominator by
790 their greatest common divisor.
793 reduce x 0 = error "{Ratio.%}: zero denominator"
794 reduce x y = (x `quot` d) :% (y `quot` d)
799 x % y = reduce (x * signum y) (abs y)
803 denominator (x:%y) = y
807 @approxRational@, applied to two real fractional numbers x and epsilon,
808 returns the simplest rational number within epsilon of x. A rational
809 number n%d in reduced form is said to be simpler than another n'%d' if
810 abs n <= abs n' && d <= d'. Any real interval contains a unique
811 simplest rational; here, for simplicity, we assume a closed rational
812 interval. If such an interval includes at least one whole number, then
813 the simplest rational is the absolutely least whole number. Otherwise,
814 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
815 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
816 the simplest rational between d'%r' and d%r.
819 approxRational x eps = simplest (x-eps) (x+eps)
820 where simplest x y | y < x = simplest y x
822 | x > 0 = simplest' n d n' d'
823 | y < 0 = - simplest' (-n') d' (-n) d
825 where xr = toRational x
832 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
834 | q /= q' = (q+1) :% 1
835 | otherwise = (q*n''+d'') :% n''
836 where (q,r) = quotRem n d
837 (q',r') = quotRem n' d'
838 nd'' = simplest' d' r' d r
840 d'' = denominator nd''
845 instance (Integral a) => Ord (Ratio a) where
846 (x:%y) <= (x':%y') = x * y' <= x' * y
847 (x:%y) < (x':%y') = x * y' < x' * y
849 instance (Integral a) => Num (Ratio a) where
850 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
851 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
852 (x:%y) * (x':%y') = reduce (x * x') (y * y')
853 negate (x:%y) = (-x) :% y
854 abs (x:%y) = abs x :% y
855 signum (x:%y) = signum x :% 1
856 fromInteger x = fromInteger x :% 1
858 instance (Integral a) => Real (Ratio a) where
859 toRational (x:%y) = toInteger x :% toInteger y
861 instance (Integral a) => Fractional (Ratio a) where
862 (x:%y) / (x':%y') = (x*y') % (y*x')
863 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
864 fromRational (x:%y) = fromInteger x :% fromInteger y
866 instance (Integral a) => RealFrac (Ratio a) where
867 properFraction (x:%y) = (fromIntegral q, r:%y)
868 where (q,r) = quotRem x y
870 instance (Integral a) => Enum (Ratio a) where
871 enumFrom = iterate ((+)1)
872 enumFromThen n m = iterate ((+)(m-n)) n
873 toEnum n = fromIntegral n :% 1
874 fromEnum = fromInteger . truncate
876 instance (Integral a) => Show (Ratio a) where
877 showsPrec p (x:%y) = showParen (p > ratio_prec)
878 (shows x . showString " % " . shows y)
880 -- defn. also used by the Read (Ratio a) instance PrelRead.
887 --Exported from std library Numeric, defined here to
888 --avoid mut. rec. between PrelNum and Numeric.
889 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
890 showSigned showPos p x
891 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
892 | otherwise = showPos x
894 showSignedInteger :: Int -> Integer -> ShowS
895 showSignedInteger p n r
896 = -- from HBC version; support code follows
897 if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
899 jtos :: Integer -> String
901 | n < 0 = '-' : jtos' (-n) []
902 | otherwise = jtos' n []
904 jtos' :: Integer -> String -> String
906 | n < 10 = chr (fromInteger (n + ord_0)) : cs
907 | otherwise = jtos' q (chr (toInt r + (ord_0::Int)) : cs)
909 (q,r) = n `quotRem` 10
911 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
913 -- These are the format types. This type is not exported.
915 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
917 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
918 formatRealFloat fmt decs x = s
921 base_i = toInteger base
925 | isInfinite x = (\ str -> if x < 0 then '-':str else str) "Infinity"
926 | x < 0 || isNegativeZero x = '-' : doFmt fmt (floatToDigits base_i (-x))
927 | otherwise = doFmt fmt (floatToDigits base_i x)
930 let ds = map intToDigit is in
933 doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
938 let e' = if e==0 then 0 else e-1 in
941 (d:ds) -> d : '.' : ds ++ "e") ++ show e'
943 let dec' = max dec 1 in
945 [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
948 (ei,is') = roundTo base (dec'+1) is
949 d:ds = map intToDigit (if ei > 0 then init is' else is')
951 d:'.':ds ++ 'e':show (e-1+ei)
954 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
959 f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
960 f n s "" = f (n-1) ('0':s) ""
961 f n s (d:ds) = f (n-1) (d:s) ds
965 let dec' = max dec 1 in
968 (ei,is') = roundTo base (dec' + e) is
969 (ls,rs) = splitAt (e+ei) (map intToDigit is')
971 mk0 ls ++ (if null rs then "" else '.':rs)
974 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
975 d:ds = map intToDigit (if ei > 0 then is' else 0:is')
980 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
991 f n [] = (0, replicate n 0)
992 f 0 (i:_) = (if i>=b2 then 1 else 0, [])
998 if i' == base then (1,0:ds) else (0,i':ds)
1001 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
1002 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
1003 -- This version uses a much slower logarithm estimator. It should be improved.
1005 -- This function returns a list of digits (Ints in [0..base-1]) and an
1007 --floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
1008 floatToDigits _ 0 = ([0], 0)
1009 floatToDigits base x =
1011 (f0, e0) = decodeFloat x
1012 (minExp0, _) = floatRange x
1015 minExp = minExp0 - p -- the real minimum exponent
1016 -- Haskell requires that f be adjusted so denormalized numbers
1017 -- will have an impossibly low exponent. Adjust for this.
1019 let n = minExp - e0 in
1020 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
1024 if f == b^(p-1) then
1025 (f*be*b*2, 2*b, be*b, b)
1029 if e > minExp && f == b^(p-1) then
1030 (f*b*2, b^(-e+1)*2, b, 1)
1032 (f*2, b^(-e)*2, 1, 1)
1036 if b == 2 && base == 10 then
1037 -- logBase 10 2 is slightly bigger than 3/10 so
1038 -- the following will err on the low side. Ignoring
1039 -- the fraction will make it err even more.
1040 -- Haskell promises that p-1 <= logBase b f < p.
1041 (p - 1 + e0) * 3 `div` 10
1043 ceiling ((log (fromInteger (f+1)) + fromInt e * log (fromInteger b)) /
1044 log (fromInteger base))
1048 if r + mUp <= expt base n * s then n else fixup (n+1)
1050 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
1054 gen ds rn sN mUpN mDnN =
1056 (dn, rn') = (rn * base) `divMod` sN
1060 case (rn' < mDnN', rn' + mUpN' > sN) of
1061 (True, False) -> dn : ds
1062 (False, True) -> dn+1 : ds
1063 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
1064 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
1068 gen [] r (s * expt base k) mUp mDn
1070 let bk = expt base (-k) in
1071 gen [] (r * bk) s (mUp * bk) (mDn * bk)
1073 (map toInt (reverse rds), k)
1077 @showRational@ converts a Rational to a string that looks like a
1078 floating point number, but without converting to any floating type
1079 (because of the possible overflow).
1081 From/by Lennart, 94/09/26
1084 showRational :: Int -> Rational -> String
1089 let (r', e) = normalize r
1092 startExpExp = 4 :: Int
1094 -- make sure 1 <= r < 10
1095 normalize :: Rational -> (Rational, Int)
1096 normalize r = if r < 1 then
1097 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1099 norm startExpExp r 0
1100 where norm :: Int -> Rational -> Int -> (Rational, Int)
1101 -- Invariant: r*10^e == original r
1106 in if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1109 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1111 prR :: Int -> Rational -> Int -> String
1112 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
1113 prR n r e | r >= 10 = prR n (r/10) (e+1)
1115 let s = show ((round (r * 10^n))::Integer)
1117 in if e > 0 && e < 8 then
1118 take e s ++ "." ++ drop0 (drop e s)
1119 else if e <= 0 && e > -3 then
1120 "0." ++ take (-e) (repeat '0') ++ drop0 s
1122 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1126 [In response to a request for documentation of how fromRational works,
1127 Joe Fasel writes:] A quite reasonable request! This code was added to
1128 the Prelude just before the 1.2 release, when Lennart, working with an
1129 early version of hbi, noticed that (read . show) was not the identity
1130 for floating-point numbers. (There was a one-bit error about half the
1131 time.) The original version of the conversion function was in fact
1132 simply a floating-point divide, as you suggest above. The new version
1133 is, I grant you, somewhat denser.
1135 Unfortunately, Joe's code doesn't work! Here's an example:
1137 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
1142 1.8217369128763981e-300
1144 Lennart's code follows, and it works...
1147 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
1148 fromRat :: (RealFloat a) => Rational -> a
1152 -- If the exponent of the nearest floating-point number to x
1153 -- is e, then the significand is the integer nearest xb^(-e),
1154 -- where b is the floating-point radix. We start with a good
1155 -- guess for e, and if it is correct, the exponent of the
1156 -- floating-point number we construct will again be e. If
1157 -- not, one more iteration is needed.
1159 f e = if e' == e then y else f e'
1160 where y = encodeFloat (round (x * (1 % b)^^e)) e
1161 (_,e') = decodeFloat y
1164 -- We obtain a trial exponent by doing a floating-point
1165 -- division of x's numerator by its denominator. The
1166 -- result of this division may not itself be the ultimate
1167 -- result, because of an accumulation of three rounding
1170 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
1171 / fromInteger (denominator x))
1174 Now, here's Lennart's code.
1177 {-# SPECIALISE fromRat ::
1179 Rational -> Float #-}
1181 --fromRat :: (RealFloat a) => Rational -> a
1183 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
1184 | x < 0 = - fromRat' (-x) -- first.
1185 | otherwise = fromRat' x
1187 -- Conversion process:
1188 -- Scale the rational number by the RealFloat base until
1189 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
1190 -- Then round the rational to an Integer and encode it with the exponent
1191 -- that we got from the scaling.
1192 -- To speed up the scaling process we compute the log2 of the number to get
1193 -- a first guess of the exponent.
1195 fromRat' :: (RealFloat a) => Rational -> a
1197 where b = floatRadix r
1200 (minExp0, _) = floatRange r
1202 minExp = minExp0 - p -- the real minimum exponent
1204 xMin = toRational (expt b (p-1))
1205 xMax = toRational (expt b p)
1207 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
1208 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
1209 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
1211 r = encodeFloat (round x') p'
1213 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
1214 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
1215 scaleRat b minExp xMin xMax p x
1216 | p <= minExp = (x, p)
1217 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
1218 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
1219 | otherwise = (x, p)
1221 -- Exponentiation with a cache for the most common numbers.
1225 expt :: Integer -> Int -> Integer
1227 | base == 2 && n >= minExpt && n <= maxExpt = expts!n
1228 | otherwise = base^n
1230 expts :: Array Int Integer
1231 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1233 -- Compute the (floor of the) log of i in base b.
1234 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1235 -- be very slow! We are just slightly more clever.
1236 integerLogBase :: Integer -> Integer -> Int
1239 | otherwise = doDiv (i `div` (b^l)) l
1241 -- Try squaring the base first to cut down the number of divisions.
1242 l = 2 * integerLogBase (b*b) i
1244 doDiv :: Integer -> Int -> Int
1245 doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
1249 %*********************************************************
1251 \subsection{Numeric primops}
1253 %*********************************************************
1255 Definitions of the boxed PrimOps; these will be
1256 used in the case of partial applications, etc.
1259 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
1260 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
1261 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
1262 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1263 negateFloat (F# x) = F# (negateFloat# x)
1265 gtFloat (F# x) (F# y) = gtFloat# x y
1266 geFloat (F# x) (F# y) = geFloat# x y
1267 eqFloat (F# x) (F# y) = eqFloat# x y
1268 neFloat (F# x) (F# y) = neFloat# x y
1269 ltFloat (F# x) (F# y) = ltFloat# x y
1270 leFloat (F# x) (F# y) = leFloat# x y
1272 float2Int (F# x) = I# (float2Int# x)
1273 int2Float (I# x) = F# (int2Float# x)
1275 expFloat (F# x) = F# (expFloat# x)
1276 logFloat (F# x) = F# (logFloat# x)
1277 sqrtFloat (F# x) = F# (sqrtFloat# x)
1278 sinFloat (F# x) = F# (sinFloat# x)
1279 cosFloat (F# x) = F# (cosFloat# x)
1280 tanFloat (F# x) = F# (tanFloat# x)
1281 asinFloat (F# x) = F# (asinFloat# x)
1282 acosFloat (F# x) = F# (acosFloat# x)
1283 atanFloat (F# x) = F# (atanFloat# x)
1284 sinhFloat (F# x) = F# (sinhFloat# x)
1285 coshFloat (F# x) = F# (coshFloat# x)
1286 tanhFloat (F# x) = F# (tanhFloat# x)
1288 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
1290 -- definitions of the boxed PrimOps; these will be
1291 -- used in the case of partial applications, etc.
1293 plusDouble (D# x) (D# y) = D# (x +## y)
1294 minusDouble (D# x) (D# y) = D# (x -## y)
1295 timesDouble (D# x) (D# y) = D# (x *## y)
1296 divideDouble (D# x) (D# y) = D# (x /## y)
1297 negateDouble (D# x) = D# (negateDouble# x)
1299 gtDouble (D# x) (D# y) = x >## y
1300 geDouble (D# x) (D# y) = x >=## y
1301 eqDouble (D# x) (D# y) = x ==## y
1302 neDouble (D# x) (D# y) = x /=## y
1303 ltDouble (D# x) (D# y) = x <## y
1304 leDouble (D# x) (D# y) = x <=## y
1306 double2Int (D# x) = I# (double2Int# x)
1307 int2Double (I# x) = D# (int2Double# x)
1308 double2Float (D# x) = F# (double2Float# x)
1309 float2Double (F# x) = D# (float2Double# x)
1311 expDouble (D# x) = D# (expDouble# x)
1312 logDouble (D# x) = D# (logDouble# x)
1313 sqrtDouble (D# x) = D# (sqrtDouble# x)
1314 sinDouble (D# x) = D# (sinDouble# x)
1315 cosDouble (D# x) = D# (cosDouble# x)
1316 tanDouble (D# x) = D# (tanDouble# x)
1317 asinDouble (D# x) = D# (asinDouble# x)
1318 acosDouble (D# x) = D# (acosDouble# x)
1319 atanDouble (D# x) = D# (atanDouble# x)
1320 sinhDouble (D# x) = D# (sinhDouble# x)
1321 coshDouble (D# x) = D# (coshDouble# x)
1322 tanhDouble (D# x) = D# (tanhDouble# x)
1324 powerDouble (D# x) (D# y) = D# (x **## y)