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"
235 x `div` y = if x > 0 && y < 0 then quotInt (x-y-1) y
236 else if x < 0 && y > 0 then quotInt (x-y+1) y
238 x `mod` y = if x > 0 && y < 0 || x < 0 && y > 0 then
239 if r/=0 then r+y else 0
244 divMod x@(I# _) y@(I# _) = (x `div` y, x `mod` y)
245 -- Stricter. Sorry if you don't like it. (WDP 94/10)
247 --OLD: even x = eqInt (x `mod` 2) 0
248 --OLD: odd x = neInt (x `mod` 2) 0
250 toInteger (I# n#) = int2Integer# n# -- give back a full-blown Integer
256 %*********************************************************
258 \subsection{Type @Integer@}
260 %*********************************************************
262 These types are used to return from integer primops
265 data Return2GMPs = Return2GMPs Int# Int# ByteArray# Int# Int# ByteArray#
266 data ReturnIntAndGMP = ReturnIntAndGMP Int# Int# Int# ByteArray#
272 instance Eq Integer where
273 (J# a1 s1 d1) == (J# a2 s2 d2)
274 = (cmpInteger# a1 s1 d1 a2 s2 d2) ==# 0#
276 (J# a1 s1 d1) /= (J# a2 s2 d2)
277 = (cmpInteger# a1 s1 d1 a2 s2 d2) /=# 0#
279 instance Ord Integer where
280 (J# a1 s1 d1) <= (J# a2 s2 d2)
281 = (cmpInteger# a1 s1 d1 a2 s2 d2) <=# 0#
283 (J# a1 s1 d1) < (J# a2 s2 d2)
284 = (cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#
286 (J# a1 s1 d1) >= (J# a2 s2 d2)
287 = (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
289 (J# a1 s1 d1) > (J# a2 s2 d2)
290 = (cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#
292 x@(J# a1 s1 d1) `max` y@(J# a2 s2 d2)
293 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#) then x else y
295 x@(J# a1 s1 d1) `min` y@(J# a2 s2 d2)
296 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#) then x else y
298 compare (J# a1 s1 d1) (J# a2 s2 d2)
299 = case cmpInteger# a1 s1 d1 a2 s2 d2 of { res# ->
300 if res# <# 0# then LT else
301 if res# ># 0# then GT else EQ
304 instance Num Integer where
305 (+) (J# a1 s1 d1) (J# a2 s2 d2)
306 = plusInteger# a1 s1 d1 a2 s2 d2
308 (-) (J# a1 s1 d1) (J# a2 s2 d2)
309 = minusInteger# a1 s1 d1 a2 s2 d2
311 negate (J# a s d) = negateInteger# a s d
313 (*) (J# a1 s1 d1) (J# a2 s2 d2)
314 = timesInteger# a1 s1 d1 a2 s2 d2
316 -- ORIG: abs n = if n >= 0 then n else -n
319 = case 0 of { J# a2 s2 d2 ->
320 if (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
322 else negateInteger# a1 s1 d1
325 signum n@(J# a1 s1 d1)
326 = case 0 of { J# a2 s2 d2 ->
328 cmp = cmpInteger# a1 s1 d1 a2 s2 d2
331 else if cmp ==# 0# then 0
337 fromInt (I# n#) = int2Integer# n# -- gives back a full-blown Integer
339 instance Real Integer where
342 instance Integral Integer where
343 quotRem (J# a1 s1 d1) (J# a2 s2 d2)
344 = case (quotRemInteger# a1 s1 d1 a2 s2 d2) of
345 Return2GMPs a3 s3 d3 a4 s4 d4
346 -> (J# a3 s3 d3, J# a4 s4 d4)
348 {- USING THE UNDERLYING "GMP" CODE IS DUBIOUS FOR NOW:
350 divMod (J# a1 s1 d1) (J# a2 s2 d2)
351 = case (divModInteger# a1 s1 d1 a2 s2 d2) of
352 Return2GMPs a3 s3 d3 a4 s4 d4
353 -> (J# a3 s3 d3, J# a4 s4 d4)
356 toInt (J# a s d) = case (integer2Int# a s d) of { n# -> I# n# }
358 -- the rest are identical to the report default methods;
359 -- you get slightly better code if you let the compiler
360 -- see them right here:
361 n `quot` d = if d /= 0 then q else
362 error "Integral.Integer.quot{PreludeCore}: divide by 0\n"
363 where (q,r) = quotRem n d
364 n `rem` d = if d /= 0 then r else
365 error "Integral.Integer.quot{PreludeCore}: divide by 0\n"
366 where (q,r) = quotRem n d
367 n `div` d = q where (q,r) = divMod n d
368 n `mod` d = r where (q,r) = divMod n d
370 divMod n d = case (quotRem n d) of { qr@(q,r) ->
371 if signum r == negate (signum d) then (q - 1, r+d) else qr }
372 -- Case-ified by WDP 94/10
374 instance Enum Integer where
375 toEnum n = toInteger n
377 enumFrom n = n : enumFrom (n + 1)
378 enumFromThen m n = en' m (n - m)
379 where en' m n = m : en' (m + n) n
380 enumFromTo n m = takeWhile (<= m) (enumFrom n)
381 enumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
384 instance Show Integer where
385 showsPrec x = showSignedInteger x
386 showList = showList__ (showsPrec 0)
388 instance Ix Integer where
391 | inRange b i = fromInteger (i - m)
392 | otherwise = error "Integer.index: Index out of range."
393 inRange (m,n) i = m <= i && i <= n
395 integer_0, integer_1, integer_2, integer_m1 :: Integer
396 integer_0 = int2Integer# 0#
397 integer_1 = int2Integer# 1#
398 integer_2 = int2Integer# 2#
399 integer_m1 = int2Integer# (negateInt# 1#)
403 %*********************************************************
405 \subsection{Type @Float@}
407 %*********************************************************
410 instance Eq Float where
411 (F# x) == (F# y) = x `eqFloat#` y
413 instance Ord Float where
414 (F# x) `compare` (F# y) | x `ltFloat#` y = LT
415 | x `eqFloat#` y = EQ
418 (F# x) < (F# y) = x `ltFloat#` y
419 (F# x) <= (F# y) = x `leFloat#` y
420 (F# x) >= (F# y) = x `geFloat#` y
421 (F# x) > (F# y) = x `gtFloat#` y
423 instance Num Float where
424 (+) x y = plusFloat x y
425 (-) x y = minusFloat x y
426 negate x = negateFloat x
427 (*) x y = timesFloat x y
429 | otherwise = negateFloat x
430 signum x | x == 0.0 = 0
432 | otherwise = negate 1
433 fromInteger n = encodeFloat n 0
434 fromInt i = int2Float i
436 instance Real Float where
437 toRational x = (m%1)*(b%1)^^n
438 where (m,n) = decodeFloat x
441 instance Fractional Float where
442 (/) x y = divideFloat x y
443 fromRational x = fromRat x
446 instance Floating Float where
447 pi = 3.141592653589793238
460 (**) x y = powerFloat x y
461 logBase x y = log y / log x
463 asinh x = log (x + sqrt (1.0+x*x))
464 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
465 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
467 instance RealFrac Float where
469 {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
470 {-# SPECIALIZE truncate :: Float -> Int #-}
471 {-# SPECIALIZE round :: Float -> Int #-}
472 {-# SPECIALIZE ceiling :: Float -> Int #-}
473 {-# SPECIALIZE floor :: Float -> Int #-}
475 {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
476 {-# SPECIALIZE truncate :: Float -> Integer #-}
477 {-# SPECIALIZE round :: Float -> Integer #-}
478 {-# SPECIALIZE ceiling :: Float -> Integer #-}
479 {-# SPECIALIZE floor :: Float -> Integer #-}
482 = case (decodeFloat x) of { (m,n) ->
483 let b = floatRadix x in
485 (fromInteger m * fromInteger b ^ n, 0.0)
487 case (quotRem m (b^(negate n))) of { (w,r) ->
488 (fromInteger w, encodeFloat r n)
492 truncate x = case properFraction x of
495 round x = case properFraction x of
497 m = if r < 0.0 then n - 1 else n + 1
498 half_down = abs r - 0.5
500 case (compare half_down 0.0) of
502 EQ -> if even n then n else m
505 ceiling x = case properFraction x of
506 (n,r) -> if r > 0.0 then n + 1 else n
508 floor x = case properFraction x of
509 (n,r) -> if r < 0.0 then n - 1 else n
511 instance RealFloat Float where
512 floatRadix _ = FLT_RADIX -- from float.h
513 floatDigits _ = FLT_MANT_DIG -- ditto
514 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
517 = case decodeFloat# f# of
518 ReturnIntAndGMP exp# a# s# d# ->
519 (J# a# s# d#, I# exp#)
521 encodeFloat (J# a# s# d#) (I# e#)
522 = case encodeFloat# a# s# d# e# of { flt# -> F# flt# }
524 exponent x = case decodeFloat x of
525 (m,n) -> if m == 0 then 0 else n + floatDigits x
527 significand x = case decodeFloat x of
528 (m,_) -> encodeFloat m (negate (floatDigits x))
530 scaleFloat k x = case decodeFloat x of
531 (m,n) -> encodeFloat m (n+k)
533 (0::Int) /= unsafePerformIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
535 (0::Int) /= unsafePerformIO (_ccall_ isFloatInfinite x) {- ditto! -}
537 (0::Int) /= unsafePerformIO (_ccall_ isFloatDenormalized x) -- ..
539 (0::Int) /= unsafePerformIO (_ccall_ isFloatNegativeZero x) -- ...
542 instance Show Float where
543 showsPrec x = showSigned showFloat x
544 showList = showList__ (showsPrec 0)
547 %*********************************************************
549 \subsection{Type @Double@}
551 %*********************************************************
554 instance Eq Double where
555 (D# x) == (D# y) = x ==## y
557 instance Ord Double where
558 (D# x) `compare` (D# y) | x <## y = LT
562 (D# x) < (D# y) = x <## y
563 (D# x) <= (D# y) = x <=## y
564 (D# x) >= (D# y) = x >=## y
565 (D# x) > (D# y) = x >## y
567 instance Num Double where
568 (+) x y = plusDouble x y
569 (-) x y = minusDouble x y
570 negate x = negateDouble x
571 (*) x y = timesDouble x y
573 | otherwise = negateDouble x
574 signum x | x == 0.0 = 0
576 | otherwise = negate 1
577 fromInteger n = encodeFloat n 0
578 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
580 instance Real Double where
581 toRational x = (m%1)*(b%1)^^n
582 where (m,n) = decodeFloat x
585 instance Fractional Double where
586 (/) x y = divideDouble x y
587 fromRational x = fromRat x
590 instance Floating Double where
591 pi = 3.141592653589793238
594 sqrt x = sqrtDouble x
598 asin x = asinDouble x
599 acos x = acosDouble x
600 atan x = atanDouble x
601 sinh x = sinhDouble x
602 cosh x = coshDouble x
603 tanh x = tanhDouble x
604 (**) x y = powerDouble x y
605 logBase x y = log y / log x
607 asinh x = log (x + sqrt (1.0+x*x))
608 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
609 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
611 instance RealFrac Double where
613 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
614 {-# SPECIALIZE truncate :: Double -> Int #-}
615 {-# SPECIALIZE round :: Double -> Int #-}
616 {-# SPECIALIZE ceiling :: Double -> Int #-}
617 {-# SPECIALIZE floor :: Double -> Int #-}
619 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
620 {-# SPECIALIZE truncate :: Double -> Integer #-}
621 {-# SPECIALIZE round :: Double -> Integer #-}
622 {-# SPECIALIZE ceiling :: Double -> Integer #-}
623 {-# SPECIALIZE floor :: Double -> Integer #-}
625 #if defined(__UNBOXED_INSTANCES__)
626 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
627 {-# SPECIALIZE truncate :: Double -> Int# #-}
628 {-# SPECIALIZE round :: Double -> Int# #-}
629 {-# SPECIALIZE ceiling :: Double -> Int# #-}
630 {-# SPECIALIZE floor :: Double -> Int# #-}
634 = case (decodeFloat x) of { (m,n) ->
635 let b = floatRadix x in
637 (fromInteger m * fromInteger b ^ n, 0.0)
639 case (quotRem m (b^(negate n))) of { (w,r) ->
640 (fromInteger w, encodeFloat r n)
644 truncate x = case properFraction x of
647 round x = case properFraction x of
649 m = if r < 0.0 then n - 1 else n + 1
650 half_down = abs r - 0.5
652 case (compare half_down 0.0) of
654 EQ -> if even n then n else m
657 ceiling x = case properFraction x of
658 (n,r) -> if r > 0.0 then n + 1 else n
660 floor x = case properFraction x of
661 (n,r) -> if r < 0.0 then n - 1 else n
663 instance RealFloat Double where
664 floatRadix _ = FLT_RADIX -- from float.h
665 floatDigits _ = DBL_MANT_DIG -- ditto
666 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
669 = case decodeDouble# d# of
670 ReturnIntAndGMP exp# a# s# d# ->
671 (J# a# s# d#, I# exp#)
673 encodeFloat (J# a# s# d#) (I# e#)
674 = case encodeDouble# a# s# d# e# of { dbl# -> D# dbl# }
676 exponent x = case decodeFloat x of
677 (m,n) -> if m == 0 then 0 else n + floatDigits x
679 significand x = case decodeFloat x of
680 (m,_) -> encodeFloat m (negate (floatDigits x))
682 scaleFloat k x = case decodeFloat x of
683 (m,n) -> encodeFloat m (n+k)
685 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
687 (0::Int) /= unsafePerformIO (_ccall_ isDoubleInfinite x) {- ditto -}
689 (0::Int) /= unsafePerformIO (_ccall_ isDoubleDenormalized x) -- ..
691 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNegativeZero x) -- ...
694 instance Show Double where
695 showsPrec x = showSigned showFloat x
696 showList = showList__ (showsPrec 0)
700 %*********************************************************
702 \subsection{Common code for @Float@ and @Double@}
704 %*********************************************************
706 The @Enum@ instances for Floats and Doubles are slightly unusual.
707 The @toEnum@ function truncates numbers to Int. The definitions
708 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
709 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
710 dubious. This example may have either 10 or 11 elements, depending on
711 how 0.1 is represented.
713 NOTE: The instances for Float and Double do not make use of the default
714 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
715 a `non-lossy' conversion to and from Ints. Instead we make use of the
716 1.2 default methods (back in the days when Enum had Ord as a superclass)
717 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
720 instance Enum Float where
721 toEnum = fromIntegral
722 fromEnum = fromInteger . truncate -- may overflow
723 enumFrom = numericEnumFrom
724 enumFromThen = numericEnumFromThen
725 enumFromThenTo = numericEnumFromThenTo
727 instance Enum Double where
728 toEnum = fromIntegral
729 fromEnum = fromInteger . truncate -- may overflow
730 enumFrom = numericEnumFrom
731 enumFromThen = numericEnumFromThen
732 enumFromThenTo = numericEnumFromThenTo
734 numericEnumFrom :: (Real a) => a -> [a]
735 numericEnumFromThen :: (Real a) => a -> a -> [a]
736 numericEnumFromThenTo :: (Real a) => a -> a -> a -> [a]
737 numericEnumFrom = iterate (+1)
738 numericEnumFromThen n m = iterate (+(m-n)) n
739 numericEnumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
740 (numericEnumFromThen n m)
744 %*********************************************************
746 \subsection{The @Ratio@ and @Rational@ types}
748 %*********************************************************
751 data (Eval a, Integral a) => Ratio a = !a :% !a deriving (Eq)
752 type Rational = Ratio Integer
756 {-# SPECIALISE (%) :: Integer -> Integer -> Rational #-}
758 (%) :: (Integral a) => a -> a -> Ratio a
759 numerator, denominator :: (Integral a) => Ratio a -> a
760 approxRational :: (RealFrac a) => a -> a -> Rational
764 \tr{reduce} is a subsidiary function used only in this module .
765 It normalises a ratio by dividing both numerator and denominator by
766 their greatest common divisor.
769 reduce x 0 = error "{Ratio.%}: zero denominator"
770 reduce x y = (x `quot` d) :% (y `quot` d)
775 x % y = reduce (x * signum y) (abs y)
779 denominator (x:%y) = y
783 @approxRational@, applied to two real fractional numbers x and epsilon,
784 returns the simplest rational number within epsilon of x. A rational
785 number n%d in reduced form is said to be simpler than another n'%d' if
786 abs n <= abs n' && d <= d'. Any real interval contains a unique
787 simplest rational; here, for simplicity, we assume a closed rational
788 interval. If such an interval includes at least one whole number, then
789 the simplest rational is the absolutely least whole number. Otherwise,
790 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
791 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
792 the simplest rational between d'%r' and d%r.
795 approxRational x eps = simplest (x-eps) (x+eps)
796 where simplest x y | y < x = simplest y x
798 | x > 0 = simplest' n d n' d'
799 | y < 0 = - simplest' (-n') d' (-n) d
801 where xr = toRational x
808 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
810 | q /= q' = (q+1) :% 1
811 | otherwise = (q*n''+d'') :% n''
812 where (q,r) = quotRem n d
813 (q',r') = quotRem n' d'
814 nd'' = simplest' d' r' d r
816 d'' = denominator nd''
821 instance (Integral a) => Ord (Ratio a) where
822 (x:%y) <= (x':%y') = x * y' <= x' * y
823 (x:%y) < (x':%y') = x * y' < x' * y
825 instance (Integral a) => Num (Ratio a) where
826 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
827 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
828 (x:%y) * (x':%y') = reduce (x * x') (y * y')
829 negate (x:%y) = (-x) :% y
830 abs (x:%y) = abs x :% y
831 signum (x:%y) = signum x :% 1
832 fromInteger x = fromInteger x :% 1
834 instance (Integral a) => Real (Ratio a) where
835 toRational (x:%y) = toInteger x :% toInteger y
837 instance (Integral a) => Fractional (Ratio a) where
838 (x:%y) / (x':%y') = (x*y') % (y*x')
839 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
840 fromRational (x:%y) = fromInteger x :% fromInteger y
842 instance (Integral a) => RealFrac (Ratio a) where
843 properFraction (x:%y) = (fromIntegral q, r:%y)
844 where (q,r) = quotRem x y
846 instance (Integral a) => Enum (Ratio a) where
847 enumFrom = iterate ((+)1)
848 enumFromThen n m = iterate ((+)(m-n)) n
849 toEnum n = fromIntegral n :% 1
850 fromEnum = fromInteger . truncate
855 instance (Integral a) => Show (Ratio a) where
856 showsPrec p (x:%y) = showParen (p > ratio_prec)
857 (shows x . showString " % " . shows y)
861 --Exported from std library Numeric, defined here to
862 --avoid mut. rec. between PrelNum and Numeric.
863 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
864 showSigned showPos p x = if x < 0 then showParen (p > 6)
865 (showChar '-' . showPos (-x))
868 showSignedInteger :: Int -> Integer -> ShowS
869 showSignedInteger p n r
870 = -- from HBC version; support code follows
871 if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
873 jtos :: Integer -> String
880 jtos' :: Integer -> String -> String
883 chr (fromInteger (n + ord_0)) : cs
885 jtos' q (chr (toInt r + (ord_0::Int)) : cs)
887 (q,r) = n `quotRem` 10
889 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
891 -- These are the format types. This type is not exported.
893 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
895 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
896 formatRealFloat fmt decs x = s
903 if x < 0 then "-Infinity" else "Infinity"
905 if x < 0 || isNegativeZero x then
906 '-':doFmt fmt (floatToDigits (toInteger base) (-x))
908 doFmt fmt (floatToDigits (toInteger base) x)
911 let ds = map intToDigit is in
914 doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
919 let e' = if e==0 then 0 else e-1 in
922 (d:ds) -> d : '.' : ds ++ "e") ++ show e'
924 let dec' = max dec 1 in
926 [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
929 (ei,is') = roundTo base (dec'+1) is
930 d:ds = map intToDigit (if ei > 0 then init is' else is')
932 d:'.':ds ++ 'e':show (e-1+ei)
935 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
940 f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
941 f n s "" = f (n-1) ('0':s) ""
942 f n s (d:ds) = f (n-1) (d:s) ds
946 let dec' = max dec 1 in
949 (ei,is') = roundTo base (dec' + e) is
950 (ls,rs) = splitAt (e+ei) (map intToDigit is')
952 mk0 ls ++ (if null rs then "" else '.':rs)
955 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
956 d:ds = map intToDigit (if ei > 0 then is' else 0:is')
961 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
972 f n [] = (0, replicate n 0)
973 f 0 (i:_) = (if i>=b2 then 1 else 0, [])
979 if i' == base then (1,0:ds) else (0,i':ds)
982 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
983 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
984 -- This version uses a much slower logarithm estimator. It should be improved.
986 -- This function returns a list of digits (Ints in [0..base-1]) and an
988 --floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
989 floatToDigits _ 0 = ([0], 0)
990 floatToDigits base x =
992 (f0, e0) = decodeFloat x
993 (minExp0, _) = floatRange x
996 minExp = minExp0 - p -- the real minimum exponent
997 -- Haskell requires that f be adjusted so denormalized numbers
998 -- will have an impossibly low exponent. Adjust for this.
1000 let n = minExp - e0 in
1001 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
1005 if f == b^(p-1) then
1006 (f*be*b*2, 2*b, be*b, b)
1010 if e > minExp && f == b^(p-1) then
1011 (f*b*2, b^(-e+1)*2, b, 1)
1013 (f*2, b^(-e)*2, 1, 1)
1017 if b == 2 && base == 10 then
1018 -- logBase 10 2 is slightly bigger than 3/10 so
1019 -- the following will err on the low side. Ignoring
1020 -- the fraction will make it err even more.
1021 -- Haskell promises that p-1 <= logBase b f < p.
1022 (p - 1 + e0) * 3 `div` 10
1024 ceiling ((log (fromInteger (f+1)) +
1025 fromInt e * log (fromInteger b)) /
1026 fromInt e * log (fromInteger b))
1030 if r + mUp <= expt base n * s then n else fixup (n+1)
1032 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
1036 gen ds rn sN mUpN mDnN =
1038 (dn, rn') = (rn * base) `divMod` sN
1042 case (rn' < mDnN', rn' + mUpN' > sN) of
1043 (True, False) -> dn : ds
1044 (False, True) -> dn+1 : ds
1045 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
1046 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
1050 gen [] r (s * expt base k) mUp mDn
1052 let bk = expt base (-k) in
1053 gen [] (r * bk) s (mUp * bk) (mDn * bk)
1055 (map toInt (reverse rds), k)
1059 @showRational@ converts a Rational to a string that looks like a
1060 floating point number, but without converting to any floating type
1061 (because of the possible overflow).
1063 From/by Lennart, 94/09/26
1066 showRational :: Int -> Rational -> String
1071 let (r', e) = normalize r
1074 startExpExp = 4 :: Int
1076 -- make sure 1 <= r < 10
1077 normalize :: Rational -> (Rational, Int)
1078 normalize r = if r < 1 then
1079 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1081 norm startExpExp r 0
1082 where norm :: Int -> Rational -> Int -> (Rational, Int)
1083 -- Invariant: r*10^e == original r
1088 in if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1091 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1093 prR :: Int -> Rational -> Int -> String
1094 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
1095 prR n r e | r >= 10 = prR n (r/10) (e+1)
1097 let s = show ((round (r * 10^n))::Integer)
1099 in if e > 0 && e < 8 then
1100 take e s ++ "." ++ drop0 (drop e s)
1101 else if e <= 0 && e > -3 then
1102 "0." ++ take (-e) (repeat '0') ++ drop0 s
1104 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1108 [In response to a request for documentation of how fromRational works,
1109 Joe Fasel writes:] A quite reasonable request! This code was added to
1110 the Prelude just before the 1.2 release, when Lennart, working with an
1111 early version of hbi, noticed that (read . show) was not the identity
1112 for floating-point numbers. (There was a one-bit error about half the
1113 time.) The original version of the conversion function was in fact
1114 simply a floating-point divide, as you suggest above. The new version
1115 is, I grant you, somewhat denser.
1117 Unfortunately, Joe's code doesn't work! Here's an example:
1119 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
1124 1.8217369128763981e-300
1126 Lennart's code follows, and it works...
1129 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
1130 fromRat :: (RealFloat a) => Rational -> a
1134 -- If the exponent of the nearest floating-point number to x
1135 -- is e, then the significand is the integer nearest xb^(-e),
1136 -- where b is the floating-point radix. We start with a good
1137 -- guess for e, and if it is correct, the exponent of the
1138 -- floating-point number we construct will again be e. If
1139 -- not, one more iteration is needed.
1141 f e = if e' == e then y else f e'
1142 where y = encodeFloat (round (x * (1 % b)^^e)) e
1143 (_,e') = decodeFloat y
1146 -- We obtain a trial exponent by doing a floating-point
1147 -- division of x's numerator by its denominator. The
1148 -- result of this division may not itself be the ultimate
1149 -- result, because of an accumulation of three rounding
1152 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
1153 / fromInteger (denominator x))
1156 Now, here's Lennart's code.
1159 {-# SPECIALISE fromRat ::
1161 Rational -> Float #-}
1163 --fromRat :: (RealFloat a) => Rational -> a
1165 if x == 0 then encodeFloat 0 0 -- Handle exceptional cases
1166 else if x < 0 then - fromRat' (-x) -- first.
1169 -- Conversion process:
1170 -- Scale the rational number by the RealFloat base until
1171 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
1172 -- Then round the rational to an Integer and encode it with the exponent
1173 -- that we got from the scaling.
1174 -- To speed up the scaling process we compute the log2 of the number to get
1175 -- a first guess of the exponent.
1177 fromRat' :: (RealFloat a) => Rational -> a
1179 where b = floatRadix r
1181 (minExp0, _) = floatRange r
1182 minExp = minExp0 - p -- the real minimum exponent
1183 xMin = toRational (expt b (p-1))
1184 xMax = toRational (expt b p)
1185 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
1186 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
1187 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
1188 r = encodeFloat (round x') p'
1190 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
1191 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
1192 scaleRat b minExp xMin xMax p x =
1195 else if x >= xMax then
1196 scaleRat b minExp xMin xMax (p+1) (x/b)
1197 else if x < xMin then
1198 scaleRat b minExp xMin xMax (p-1) (x*b)
1202 -- Exponentiation with a cache for the most common numbers.
1205 expt :: Integer -> Int -> Integer
1207 if base == 2 && n >= minExpt && n <= maxExpt then
1211 expts :: Array Int Integer
1212 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1214 -- Compute the (floor of the) log of i in base b.
1215 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1216 -- be very slow! We are just slightly more clever.
1217 integerLogBase :: Integer -> Integer -> Int
1218 integerLogBase b i =
1222 -- Try squaring the base first to cut down the number of divisions.
1223 let l = 2 * integerLogBase (b*b) i
1224 doDiv :: Integer -> Int -> Int
1225 doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
1226 in doDiv (i `div` (b^l)) l
1229 %*********************************************************
1231 \subsection{Numeric primops}
1233 %*********************************************************
1235 Definitions of the boxed PrimOps; these will be
1236 used in the case of partial applications, etc.
1239 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
1240 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
1241 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
1242 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1243 negateFloat (F# x) = F# (negateFloat# x)
1245 gtFloat (F# x) (F# y) = gtFloat# x y
1246 geFloat (F# x) (F# y) = geFloat# x y
1247 eqFloat (F# x) (F# y) = eqFloat# x y
1248 neFloat (F# x) (F# y) = neFloat# x y
1249 ltFloat (F# x) (F# y) = ltFloat# x y
1250 leFloat (F# x) (F# y) = leFloat# x y
1252 float2Int (F# x) = I# (float2Int# x)
1253 int2Float (I# x) = F# (int2Float# x)
1255 expFloat (F# x) = F# (expFloat# x)
1256 logFloat (F# x) = F# (logFloat# x)
1257 sqrtFloat (F# x) = F# (sqrtFloat# x)
1258 sinFloat (F# x) = F# (sinFloat# x)
1259 cosFloat (F# x) = F# (cosFloat# x)
1260 tanFloat (F# x) = F# (tanFloat# x)
1261 asinFloat (F# x) = F# (asinFloat# x)
1262 acosFloat (F# x) = F# (acosFloat# x)
1263 atanFloat (F# x) = F# (atanFloat# x)
1264 sinhFloat (F# x) = F# (sinhFloat# x)
1265 coshFloat (F# x) = F# (coshFloat# x)
1266 tanhFloat (F# x) = F# (tanhFloat# x)
1268 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
1270 -- definitions of the boxed PrimOps; these will be
1271 -- used in the case of partial applications, etc.
1273 plusDouble (D# x) (D# y) = D# (x +## y)
1274 minusDouble (D# x) (D# y) = D# (x -## y)
1275 timesDouble (D# x) (D# y) = D# (x *## y)
1276 divideDouble (D# x) (D# y) = D# (x /## y)
1277 negateDouble (D# x) = D# (negateDouble# x)
1279 gtDouble (D# x) (D# y) = x >## y
1280 geDouble (D# x) (D# y) = x >=## y
1281 eqDouble (D# x) (D# y) = x ==## y
1282 neDouble (D# x) (D# y) = x /=## y
1283 ltDouble (D# x) (D# y) = x <## y
1284 leDouble (D# x) (D# y) = x <=## y
1286 double2Int (D# x) = I# (double2Int# x)
1287 int2Double (I# x) = D# (int2Double# x)
1288 double2Float (D# x) = F# (double2Float# x)
1289 float2Double (F# x) = D# (float2Double# x)
1291 expDouble (D# x) = D# (expDouble# x)
1292 logDouble (D# x) = D# (logDouble# x)
1293 sqrtDouble (D# x) = D# (sqrtDouble# x)
1294 sinDouble (D# x) = D# (sinDouble# x)
1295 cosDouble (D# x) = D# (cosDouble# x)
1296 tanDouble (D# x) = D# (tanDouble# x)
1297 asinDouble (D# x) = D# (asinDouble# x)
1298 acosDouble (D# x) = D# (acosDouble# x)
1299 atanDouble (D# x) = D# (atanDouble# x)
1300 sinhDouble (D# x) = D# (sinhDouble# x)
1301 coshDouble (D# x) = D# (coshDouble# x)
1302 tanhDouble (D# x) = D# (tanhDouble# x)
1304 powerDouble (D# x) (D# y) = D# (x **## y)