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 #-} IOBase ( error )
27 import ArrBase ( Array, array, (!) )
28 import STBase ( unsafePerformPrimIO )
30 import Foreign () -- This import tells the dependency analyser to compile Foreign first.
31 -- There's an implicit dependency on Foreign because the ccalls in
32 -- PrelNum implicitly mention CCallable.
35 infixl 7 /, %, `quot`, `rem`, `div`, `mod`
39 %*********************************************************
41 \subsection{Standard numeric classes}
43 %*********************************************************
46 class (Num a, Ord a) => Real a where
47 toRational :: a -> Rational
49 class (Real a, Enum a) => Integral a where
50 quot, rem, div, mod :: a -> a -> a
51 quotRem, divMod :: a -> a -> (a,a)
52 toInteger :: a -> Integer
53 toInt :: a -> Int -- partain: Glasgow extension
55 n `quot` d = q where (q,r) = quotRem n d
56 n `rem` d = r where (q,r) = quotRem n d
57 n `div` d = q where (q,r) = divMod n d
58 n `mod` d = r where (q,r) = divMod n d
59 divMod n d = if signum r == negate (signum d) then (q-1, r+d) else qr
60 where qr@(q,r) = quotRem n d
62 class (Num a) => Fractional a where
65 fromRational :: Rational -> a
69 class (Fractional a) => Floating a where
71 exp, log, sqrt :: a -> a
72 (**), logBase :: a -> a -> a
73 sin, cos, tan :: a -> a
74 asin, acos, atan :: a -> a
75 sinh, cosh, tanh :: a -> a
76 asinh, acosh, atanh :: a -> a
78 x ** y = exp (log x * y)
79 logBase x y = log y / log x
82 tanh x = sinh x / cosh x
84 class (Real a, Fractional a) => RealFrac a where
85 properFraction :: (Integral b) => a -> (b,a)
86 truncate, round :: (Integral b) => a -> b
87 ceiling, floor :: (Integral b) => a -> b
89 truncate x = m where (m,_) = properFraction x
91 round x = let (n,r) = properFraction x
92 m = if r < 0 then n - 1 else n + 1
93 in case signum (abs r - 0.5) of
95 0 -> if even n then n else m
98 ceiling x = if r > 0 then n + 1 else n
99 where (n,r) = properFraction x
101 floor x = if r < 0 then n - 1 else n
102 where (n,r) = properFraction x
104 class (RealFrac a, Floating a) => RealFloat a where
105 floatRadix :: a -> Integer
106 floatDigits :: a -> Int
107 floatRange :: a -> (Int,Int)
108 decodeFloat :: a -> (Integer,Int)
109 encodeFloat :: Integer -> Int -> a
111 significand :: a -> a
112 scaleFloat :: Int -> a -> a
113 isNaN, isInfinite, isDenormalized, isNegativeZero, isIEEE
116 exponent x = if m == 0 then 0 else n + floatDigits x
117 where (m,n) = decodeFloat x
119 significand x = encodeFloat m (negate (floatDigits x))
120 where (m,_) = decodeFloat x
122 scaleFloat k x = encodeFloat m (n+k)
123 where (m,n) = decodeFloat x
126 %*********************************************************
128 \subsection{Overloaded numeric functions}
130 %*********************************************************
133 even, odd :: (Integral a) => a -> Bool
134 even n = n `rem` 2 == 0
137 {-# GENERATE_SPECS gcd a{Int#,Int,Integer} #-}
138 gcd :: (Integral a) => a -> a -> a
139 gcd 0 0 = error "Prelude.gcd: gcd 0 0 is undefined"
140 gcd x y = gcd' (abs x) (abs y)
142 gcd' x y = gcd' y (x `rem` y)
144 {-# GENERATE_SPECS lcm a{Int#,Int,Integer} #-}
145 lcm :: (Integral a) => a -> a -> a
148 lcm x y = abs ((x `quot` (gcd x y)) * y)
150 (^) :: (Num a, Integral b) => a -> b -> a
152 x ^ n | n > 0 = f x (n-1) x
154 f x n y = g x n where
155 g x n | even n = g (x*x) (n `quot` 2)
156 | otherwise = f x (n-1) (x*y)
157 _ ^ _ = error "Prelude.^: negative exponent"
159 (^^) :: (Fractional a, Integral b) => a -> b -> a
160 x ^^ n = if n >= 0 then x^n else recip (x^(negate n))
162 fromIntegral :: (Integral a, Num b) => a -> b
163 fromIntegral = fromInteger . toInteger
165 fromRealFrac :: (RealFrac a, Fractional b) => a -> b
166 fromRealFrac = fromRational . toRational
168 atan2 :: (RealFloat a) => a -> a -> a
169 atan2 y x = case (signum y, signum x) of
173 (-1, 0) -> (negate pi)/2
174 ( _, 1) -> atan (y/x)
175 ( _,-1) -> atan (y/x) + pi
176 ( 0, 0) -> error "Prelude.atan2: atan2 of origin"
180 %*********************************************************
182 \subsection{Instances for @Int@}
184 %*********************************************************
187 instance Real Int where
188 toRational x = toInteger x % 1
190 instance Integral Int where
191 a@(I# _) `quotRem` b@(I# _) = (a `quotInt` b, a `remInt` b)
192 -- OK, so I made it a little stricter. Shoot me. (WDP 94/10)
194 -- following chks for zero divisor are non-standard (WDP)
195 a `quot` b = if b /= 0
197 else error "Integral.Int.quot{PreludeCore}: divide by 0\n"
198 a `rem` b = if b /= 0
200 else error "Integral.Int.rem{PreludeCore}: divide by 0\n"
202 x `div` y = if x > 0 && y < 0 then quotInt (x-y-1) y
203 else if x < 0 && y > 0 then quotInt (x-y+1) y
205 x `mod` y = if x > 0 && y < 0 || x < 0 && y > 0 then
206 if r/=0 then r+y else 0
211 divMod x@(I# _) y@(I# _) = (x `div` y, x `mod` y)
212 -- Stricter. Sorry if you don't like it. (WDP 94/10)
214 --OLD: even x = eqInt (x `mod` 2) 0
215 --OLD: odd x = neInt (x `mod` 2) 0
217 toInteger (I# n#) = int2Integer# n# -- give back a full-blown Integer
223 %*********************************************************
225 \subsection{Type @Integer@}
227 %*********************************************************
229 These types are used to return from integer primops
232 data Return2GMPs = Return2GMPs Int# Int# ByteArray# Int# Int# ByteArray#
233 data ReturnIntAndGMP = ReturnIntAndGMP Int# Int# Int# ByteArray#
239 instance Eq Integer where
240 (J# a1 s1 d1) == (J# a2 s2 d2)
241 = (cmpInteger# a1 s1 d1 a2 s2 d2) ==# 0#
243 (J# a1 s1 d1) /= (J# a2 s2 d2)
244 = (cmpInteger# a1 s1 d1 a2 s2 d2) /=# 0#
246 instance Ord Integer where
247 (J# a1 s1 d1) <= (J# a2 s2 d2)
248 = (cmpInteger# a1 s1 d1 a2 s2 d2) <=# 0#
250 (J# a1 s1 d1) < (J# a2 s2 d2)
251 = (cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#
253 (J# a1 s1 d1) >= (J# a2 s2 d2)
254 = (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
256 (J# a1 s1 d1) > (J# a2 s2 d2)
257 = (cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#
259 x@(J# a1 s1 d1) `max` y@(J# a2 s2 d2)
260 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#) then x else y
262 x@(J# a1 s1 d1) `min` y@(J# a2 s2 d2)
263 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#) then x else y
265 compare (J# a1 s1 d1) (J# a2 s2 d2)
266 = case cmpInteger# a1 s1 d1 a2 s2 d2 of { res# ->
267 if res# <# 0# then LT else
268 if res# ># 0# then GT else EQ
271 instance Num Integer where
272 (+) (J# a1 s1 d1) (J# a2 s2 d2)
273 = plusInteger# a1 s1 d1 a2 s2 d2
275 (-) (J# a1 s1 d1) (J# a2 s2 d2)
276 = minusInteger# a1 s1 d1 a2 s2 d2
278 negate (J# a s d) = negateInteger# a s d
280 (*) (J# a1 s1 d1) (J# a2 s2 d2)
281 = timesInteger# a1 s1 d1 a2 s2 d2
283 -- ORIG: abs n = if n >= 0 then n else -n
286 = case 0 of { J# a2 s2 d2 ->
287 if (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
289 else negateInteger# a1 s1 d1
292 signum n@(J# a1 s1 d1)
293 = case 0 of { J# a2 s2 d2 ->
295 cmp = cmpInteger# a1 s1 d1 a2 s2 d2
298 else if cmp ==# 0# then 0
304 fromInt (I# n#) = int2Integer# n# -- gives back a full-blown Integer
306 instance Real Integer where
309 instance Integral Integer where
310 quotRem (J# a1 s1 d1) (J# a2 s2 d2)
311 = case (quotRemInteger# a1 s1 d1 a2 s2 d2) of
312 Return2GMPs a3 s3 d3 a4 s4 d4
313 -> (J# a3 s3 d3, J# a4 s4 d4)
315 {- USING THE UNDERLYING "GMP" CODE IS DUBIOUS FOR NOW:
317 divMod (J# a1 s1 d1) (J# a2 s2 d2)
318 = case (divModInteger# a1 s1 d1 a2 s2 d2) of
319 Return2GMPs a3 s3 d3 a4 s4 d4
320 -> (J# a3 s3 d3, J# a4 s4 d4)
323 toInt (J# a s d) = case (integer2Int# a s d) of { n# -> I# n# }
325 -- the rest are identical to the report default methods;
326 -- you get slightly better code if you let the compiler
327 -- see them right here:
328 n `quot` d = q where (q,r) = quotRem n d
329 n `rem` d = r where (q,r) = quotRem n d
330 n `div` d = q where (q,r) = divMod n d
331 n `mod` d = r where (q,r) = divMod n d
333 divMod n d = case (quotRem n d) of { qr@(q,r) ->
334 if signum r == negate (signum d) then (q - 1, r+d) else qr }
335 -- Case-ified by WDP 94/10
337 instance Enum Integer where
338 toEnum n = toInteger n
340 enumFrom n = n : enumFrom (n + 1)
341 enumFromThen m n = en' m (n - m)
342 where en' m n = m : en' (m + n) n
343 enumFromTo n m = takeWhile (<= m) (enumFrom n)
344 enumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
347 instance Show Integer where
348 showsPrec x = showSignedInteger x
349 showList = showList__ (showsPrec 0)
351 instance Ix Integer where
354 | inRange b i = fromInteger (i - m)
355 | otherwise = error "Integer.index: Index out of range."
356 inRange (m,n) i = m <= i && i <= n
358 integer_0, integer_1, integer_2, integer_m1 :: Integer
359 integer_0 = int2Integer# 0#
360 integer_1 = int2Integer# 1#
361 integer_2 = int2Integer# 2#
362 integer_m1 = int2Integer# (negateInt# 1#)
366 %*********************************************************
368 \subsection{Type @Float@}
370 %*********************************************************
373 instance Eq Float where
374 (F# x) == (F# y) = x `eqFloat#` y
376 instance Ord Float where
377 (F# x) `compare` (F# y) | x `ltFloat#` y = LT
378 | x `eqFloat#` y = EQ
381 (F# x) < (F# y) = x `ltFloat#` y
382 (F# x) <= (F# y) = x `leFloat#` y
383 (F# x) >= (F# y) = x `geFloat#` y
384 (F# x) > (F# y) = x `gtFloat#` y
386 instance Num Float where
387 (+) x y = plusFloat x y
388 (-) x y = minusFloat x y
389 negate x = negateFloat x
390 (*) x y = timesFloat x y
392 | otherwise = negateFloat x
393 signum x | x == 0.0 = 0
395 | otherwise = negate 1
396 fromInteger n = encodeFloat n 0
397 fromInt i = int2Float i
399 instance Real Float where
400 toRational x = (m%1)*(b%1)^^n
401 where (m,n) = decodeFloat x
404 instance Fractional Float where
405 (/) x y = divideFloat x y
406 fromRational x = fromRat x
409 instance Floating Float where
410 pi = 3.141592653589793238
423 (**) x y = powerFloat x y
424 logBase x y = log y / log x
426 asinh x = log (x + sqrt (1.0+x*x))
427 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
428 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
430 instance RealFrac Float where
432 {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
433 {-# SPECIALIZE truncate :: Float -> Int #-}
434 {-# SPECIALIZE round :: Float -> Int #-}
435 {-# SPECIALIZE ceiling :: Float -> Int #-}
436 {-# SPECIALIZE floor :: Float -> Int #-}
438 {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
439 {-# SPECIALIZE truncate :: Float -> Integer #-}
440 {-# SPECIALIZE round :: Float -> Integer #-}
441 {-# SPECIALIZE ceiling :: Float -> Integer #-}
442 {-# SPECIALIZE floor :: Float -> Integer #-}
445 = case (decodeFloat x) of { (m,n) ->
446 let b = floatRadix x in
448 (fromInteger m * fromInteger b ^ n, 0.0)
450 case (quotRem m (b^(negate n))) of { (w,r) ->
451 (fromInteger w, encodeFloat r n)
455 truncate x = case properFraction x of
458 round x = case properFraction x of
460 m = if r < 0.0 then n - 1 else n + 1
461 half_down = abs r - 0.5
463 case (compare half_down 0.0) of
465 EQ -> if even n then n else m
468 ceiling x = case properFraction x of
469 (n,r) -> if r > 0.0 then n + 1 else n
471 floor x = case properFraction x of
472 (n,r) -> if r < 0.0 then n - 1 else n
474 instance RealFloat Float where
475 floatRadix _ = FLT_RADIX -- from float.h
476 floatDigits _ = FLT_MANT_DIG -- ditto
477 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
480 = case decodeFloat# f# of
481 ReturnIntAndGMP exp# a# s# d# ->
482 (J# a# s# d#, I# exp#)
484 encodeFloat (J# a# s# d#) (I# e#)
485 = case encodeFloat# a# s# d# e# of { flt# -> F# flt# }
487 exponent x = case decodeFloat x of
488 (m,n) -> if m == 0 then 0 else n + floatDigits x
490 significand x = case decodeFloat x of
491 (m,_) -> encodeFloat m (negate (floatDigits x))
493 scaleFloat k x = case decodeFloat x of
494 (m,n) -> encodeFloat m (n+k)
496 (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
498 (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatInfinite x) {- ditto! -}
500 (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatDenormalized x) -- ..
502 (0::Int) /= unsafePerformPrimIO (_ccall_ isFloatNegativeZero x) -- ...
505 instance Show Float where
506 showsPrec x = showSigned showFloat x
507 showList = showList__ (showsPrec 0)
510 %*********************************************************
512 \subsection{Type @Double@}
514 %*********************************************************
517 instance Eq Double where
518 (D# x) == (D# y) = x ==## y
520 instance Ord Double where
521 (D# x) `compare` (D# y) | x <## y = LT
525 (D# x) < (D# y) = x <## y
526 (D# x) <= (D# y) = x <=## y
527 (D# x) >= (D# y) = x >=## y
528 (D# x) > (D# y) = x >## y
530 instance Num Double where
531 (+) x y = plusDouble x y
532 (-) x y = minusDouble x y
533 negate x = negateDouble x
534 (*) x y = timesDouble x y
536 | otherwise = negateDouble x
537 signum x | x == 0.0 = 0
539 | otherwise = negate 1
540 fromInteger n = encodeFloat n 0
541 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
543 instance Real Double where
544 toRational x = (m%1)*(b%1)^^n
545 where (m,n) = decodeFloat x
548 instance Fractional Double where
549 (/) x y = divideDouble x y
550 fromRational x = fromRat x
553 instance Floating Double where
554 pi = 3.141592653589793238
557 sqrt x = sqrtDouble x
561 asin x = asinDouble x
562 acos x = acosDouble x
563 atan x = atanDouble x
564 sinh x = sinhDouble x
565 cosh x = coshDouble x
566 tanh x = tanhDouble x
567 (**) x y = powerDouble x y
568 logBase x y = log y / log x
570 asinh x = log (x + sqrt (1.0+x*x))
571 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
572 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
574 instance RealFrac Double where
576 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
577 {-# SPECIALIZE truncate :: Double -> Int #-}
578 {-# SPECIALIZE round :: Double -> Int #-}
579 {-# SPECIALIZE ceiling :: Double -> Int #-}
580 {-# SPECIALIZE floor :: Double -> Int #-}
582 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
583 {-# SPECIALIZE truncate :: Double -> Integer #-}
584 {-# SPECIALIZE round :: Double -> Integer #-}
585 {-# SPECIALIZE ceiling :: Double -> Integer #-}
586 {-# SPECIALIZE floor :: Double -> Integer #-}
588 #if defined(__UNBOXED_INSTANCES__)
589 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
590 {-# SPECIALIZE truncate :: Double -> Int# #-}
591 {-# SPECIALIZE round :: Double -> Int# #-}
592 {-# SPECIALIZE ceiling :: Double -> Int# #-}
593 {-# SPECIALIZE floor :: Double -> Int# #-}
597 = case (decodeFloat x) of { (m,n) ->
598 let b = floatRadix x in
600 (fromInteger m * fromInteger b ^ n, 0.0)
602 case (quotRem m (b^(negate n))) of { (w,r) ->
603 (fromInteger w, encodeFloat r n)
607 truncate x = case properFraction x of
610 round x = case properFraction x of
612 m = if r < 0.0 then n - 1 else n + 1
613 half_down = abs r - 0.5
615 case (compare half_down 0.0) of
617 EQ -> if even n then n else m
620 ceiling x = case properFraction x of
621 (n,r) -> if r > 0.0 then n + 1 else n
623 floor x = case properFraction x of
624 (n,r) -> if r < 0.0 then n - 1 else n
626 instance RealFloat Double where
627 floatRadix _ = FLT_RADIX -- from float.h
628 floatDigits _ = DBL_MANT_DIG -- ditto
629 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
632 = case decodeDouble# d# of
633 ReturnIntAndGMP exp# a# s# d# ->
634 (J# a# s# d#, I# exp#)
636 encodeFloat (J# a# s# d#) (I# e#)
637 = case encodeDouble# a# s# d# e# of { dbl# -> D# dbl# }
639 exponent x = case decodeFloat x of
640 (m,n) -> if m == 0 then 0 else n + floatDigits x
642 significand x = case decodeFloat x of
643 (m,_) -> encodeFloat m (negate (floatDigits x))
645 scaleFloat k x = case decodeFloat x of
646 (m,n) -> encodeFloat m (n+k)
648 (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
650 (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleInfinite x) {- ditto -}
652 (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleDenormalized x) -- ..
654 (0::Int) /= unsafePerformPrimIO (_ccall_ isDoubleNegativeZero x) -- ...
657 instance Show Double where
658 showsPrec x = showSigned showFloat x
659 showList = showList__ (showsPrec 0)
663 %*********************************************************
665 \subsection{Common code for @Float@ and @Double@}
667 %*********************************************************
669 The Enum instances for Floats and Doubles are slightly unusual.
670 The `toEnum' function truncates numbers to Int. The definitions
671 of enumFrom and enumFromThen allow floats to be used in arithmetic
672 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
673 dubious. This example may have either 10 or 11 elements, depending on
674 how 0.1 is represented.
677 instance Enum Float where
678 toEnum = fromIntegral
679 fromEnum = fromInteger . truncate -- may overflow
680 enumFrom = numericEnumFrom
681 enumFromThen = numericEnumFromThen
683 instance Enum Double where
684 toEnum = fromIntegral
685 fromEnum = fromInteger . truncate -- may overflow
686 enumFrom = numericEnumFrom
687 enumFromThen = numericEnumFromThen
689 numericEnumFrom :: (Real a) => a -> [a]
690 numericEnumFromThen :: (Real a) => a -> a -> [a]
691 numericEnumFrom = iterate (+1)
692 numericEnumFromThen n m = iterate (+(m-n)) n
696 %*********************************************************
698 \subsection{The @Ratio@ and @Rational@ types}
700 %*********************************************************
703 data (Integral a) => Ratio a = !a :% !a deriving (Eq)
704 type Rational = Ratio Integer
708 (%) :: (Integral a) => a -> a -> Ratio a
709 numerator, denominator :: (Integral a) => Ratio a -> a
710 approxRational :: (RealFrac a) => a -> a -> Rational
714 \tr{reduce} is a subsidiary function used only in this module .
715 It normalises a ratio by dividing both numerator and denominator by
716 their greatest common divisor.
719 reduce x 0 = error "{Ratio.%}: zero denominator"
720 reduce x y = (x `quot` d) :% (y `quot` d)
725 x % y = reduce (x * signum y) (abs y)
729 denominator (x:%y) = y
733 @approxRational@, applied to two real fractional numbers x and epsilon,
734 returns the simplest rational number within epsilon of x. A rational
735 number n%d in reduced form is said to be simpler than another n'%d' if
736 abs n <= abs n' && d <= d'. Any real interval contains a unique
737 simplest rational; here, for simplicity, we assume a closed rational
738 interval. If such an interval includes at least one whole number, then
739 the simplest rational is the absolutely least whole number. Otherwise,
740 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
741 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
742 the simplest rational between d'%r' and d%r.
745 approxRational x eps = simplest (x-eps) (x+eps)
746 where simplest x y | y < x = simplest y x
748 | x > 0 = simplest' n d n' d'
749 | y < 0 = - simplest' (-n') d' (-n) d
751 where xr = toRational x
758 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
760 | q /= q' = (q+1) :% 1
761 | otherwise = (q*n''+d'') :% n''
762 where (q,r) = quotRem n d
763 (q',r') = quotRem n' d'
764 nd'' = simplest' d' r' d r
766 d'' = denominator nd''
771 instance (Integral a) => Ord (Ratio a) where
772 (x:%y) <= (x':%y') = x * y' <= x' * y
773 (x:%y) < (x':%y') = x * y' < x' * y
775 instance (Integral a) => Num (Ratio a) where
776 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
777 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
778 (x:%y) * (x':%y') = reduce (x * x') (y * y')
779 negate (x:%y) = (-x) :% y
780 abs (x:%y) = abs x :% y
781 signum (x:%y) = signum x :% 1
782 fromInteger x = fromInteger x :% 1
784 instance (Integral a) => Real (Ratio a) where
785 toRational (x:%y) = toInteger x :% toInteger y
787 instance (Integral a) => Fractional (Ratio a) where
788 (x:%y) / (x':%y') = (x*y') % (y*x')
789 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
790 fromRational (x:%y) = fromInteger x :% fromInteger y
792 instance (Integral a) => RealFrac (Ratio a) where
793 properFraction (x:%y) = (fromIntegral q, r:%y)
794 where (q,r) = quotRem x y
796 instance (Integral a) => Enum (Ratio a) where
797 enumFrom = iterate ((+)1)
798 enumFromThen n m = iterate ((+)(m-n)) n
799 toEnum n = fromIntegral n :% 1
800 fromEnum = fromInteger . truncate
805 instance (Integral a) => Show (Ratio a) where
806 showsPrec p (x:%y) = showParen (p > ratio_prec)
807 (shows x . showString " % " . shows y)
811 --Exported from std library Numeric, defined here to
812 --avoid mut. rec. between PrelNum and Numeric.
813 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
814 showSigned showPos p x = if x < 0 then showParen (p > 6)
815 (showChar '-' . showPos (-x))
818 showSignedInteger :: Int -> Integer -> ShowS
819 showSignedInteger p n r
820 = -- from HBC version; support code follows
821 if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
823 jtos :: Integer -> String
830 jtos' :: Integer -> String -> String
833 chr (fromInteger (n + ord_0)) : cs
835 jtos' (n `quot` 10) (chr (fromInteger (n `rem` 10 + ord_0)) : cs)
837 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
839 -- These are the format types. This type is not exported.
841 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
843 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
844 formatRealFloat fmt decs x = s
851 if x < 0 then "-Infinity" else "Infinity"
853 if x < 0 || isNegativeZero x then
854 '-':doFmt fmt (floatToDigits (toInteger base) (-x))
856 doFmt fmt (floatToDigits (toInteger base) x)
859 let ds = map intToDigit is in
862 doFmt (if e <0 || e > 7 then FFExponent else FFFixed)
867 let e' = if e==0 then 0 else e-1 in
870 (d:ds) -> d : '.' : ds ++ "e") ++ show e'
872 let dec' = max dec 1 in
874 [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
877 (ei,is') = roundTo base (dec'+1) is
878 d:ds = map intToDigit (if ei > 0 then init is' else is')
880 d:'.':ds ++ 'e':show (e-1+ei)
883 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
888 f 0 s ds = mk0 (reverse s) ++ '.':mk0 ds
889 f n s "" = f (n-1) ('0':s) ""
890 f n s (d:ds) = f (n-1) (d:s) ds
894 let dec' = max dec 1 in
897 (ei,is') = roundTo base (dec' + e) is
898 (ls,rs) = splitAt (e+ei) (map intToDigit is')
900 mk0 ls ++ (if null rs then "" else '.':rs)
903 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
904 d:ds = map intToDigit (if ei > 0 then is' else 0:is')
909 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
920 f n [] = (0, replicate n 0)
921 f 0 (i:_) = (if i>=b2 then 1 else 0, [])
927 if i' == base then (1,0:ds) else (0,i':ds)
930 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
931 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
932 -- This version uses a much slower logarithm estimator. It should be improved.
934 -- This function returns a list of digits (Ints in [0..base-1]) and an
936 --floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
937 floatToDigits _ 0 = ([0], 0)
938 floatToDigits base x =
940 (f0, e0) = decodeFloat x
941 (minExp0, _) = floatRange x
944 minExp = minExp0 - p -- the real minimum exponent
945 -- Haskell requires that f be adjusted so denormalized numbers
946 -- will have an impossibly low exponent. Adjust for this.
948 let n = minExp - e0 in
949 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
954 (f*be*b*2, 2*b, be*b, b)
958 if e > minExp && f == b^(p-1) then
959 (f*b*2, b^(-e+1)*2, b, 1)
961 (f*2, b^(-e)*2, 1, 1)
965 if b == 2 && base == 10 then
966 -- logBase 10 2 is slightly bigger than 3/10 so
967 -- the following will err on the low side. Ignoring
968 -- the fraction will make it err even more.
969 -- Haskell promises that p-1 <= logBase b f < p.
970 (p - 1 + e0) * 3 `div` 10
972 ceiling ((log (fromInteger (f+1)) +
973 fromInt e * log (fromInteger b)) /
974 fromInt e * log (fromInteger b))
978 if r + mUp <= expt base n * s then n else fixup (n+1)
980 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
984 gen ds rn sN mUpN mDnN =
986 (dn, rn') = (rn * base) `divMod` sN
990 case (rn' < mDnN', rn' + mUpN' > sN) of
991 (True, False) -> dn : ds
992 (False, True) -> dn+1 : ds
993 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
994 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
998 gen [] r (s * expt base k) mUp mDn
1000 let bk = expt base (-k) in
1001 gen [] (r * bk) s (mUp * bk) (mDn * bk)
1003 (map toInt (reverse rds), k)
1007 @showRational@ converts a Rational to a string that looks like a
1008 floating point number, but without converting to any floating type
1009 (because of the possible overflow).
1011 From/by Lennart, 94/09/26
1014 showRational :: Int -> Rational -> String
1019 let (r', e) = normalize r
1022 startExpExp = 4 :: Int
1024 -- make sure 1 <= r < 10
1025 normalize :: Rational -> (Rational, Int)
1026 normalize r = if r < 1 then
1027 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1029 norm startExpExp r 0
1030 where norm :: Int -> Rational -> Int -> (Rational, Int)
1031 -- Invariant: r*10^e == original r
1036 in if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1039 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1041 prR :: Int -> Rational -> Int -> String
1042 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
1043 prR n r e | r >= 10 = prR n (r/10) (e+1)
1045 let s = show ((round (r * 10^n))::Integer)
1047 in if e > 0 && e < 8 then
1048 take e s ++ "." ++ drop0 (drop e s)
1049 else if e <= 0 && e > -3 then
1050 "0." ++ take (-e) (repeat '0') ++ drop0 s
1052 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1056 [In response to a request for documentation of how fromRational works,
1057 Joe Fasel writes:] A quite reasonable request! This code was added to
1058 the Prelude just before the 1.2 release, when Lennart, working with an
1059 early version of hbi, noticed that (read . show) was not the identity
1060 for floating-point numbers. (There was a one-bit error about half the
1061 time.) The original version of the conversion function was in fact
1062 simply a floating-point divide, as you suggest above. The new version
1063 is, I grant you, somewhat denser.
1065 Unfortunately, Joe's code doesn't work! Here's an example:
1067 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
1072 1.8217369128763981e-300
1074 Lennart's code follows, and it works...
1077 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
1078 fromRat :: (RealFloat a) => Rational -> a
1082 -- If the exponent of the nearest floating-point number to x
1083 -- is e, then the significand is the integer nearest xb^(-e),
1084 -- where b is the floating-point radix. We start with a good
1085 -- guess for e, and if it is correct, the exponent of the
1086 -- floating-point number we construct will again be e. If
1087 -- not, one more iteration is needed.
1089 f e = if e' == e then y else f e'
1090 where y = encodeFloat (round (x * (1 % b)^^e)) e
1091 (_,e') = decodeFloat y
1094 -- We obtain a trial exponent by doing a floating-point
1095 -- division of x's numerator by its denominator. The
1096 -- result of this division may not itself be the ultimate
1097 -- result, because of an accumulation of three rounding
1100 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
1101 / fromInteger (denominator x))
1104 Now, here's Lennart's code.
1107 --fromRat :: (RealFloat a) => Rational -> a
1109 if x == 0 then encodeFloat 0 0 -- Handle exceptional cases
1110 else if x < 0 then - fromRat' (-x) -- first.
1113 -- Conversion process:
1114 -- Scale the rational number by the RealFloat base until
1115 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
1116 -- Then round the rational to an Integer and encode it with the exponent
1117 -- that we got from the scaling.
1118 -- To speed up the scaling process we compute the log2 of the number to get
1119 -- a first guess of the exponent.
1121 fromRat' :: (RealFloat a) => Rational -> a
1123 where b = floatRadix r
1125 (minExp0, _) = floatRange r
1126 minExp = minExp0 - p -- the real minimum exponent
1127 xMin = toRational (expt b (p-1))
1128 xMax = toRational (expt b p)
1129 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
1130 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
1131 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
1132 r = encodeFloat (round x') p'
1134 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
1135 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
1136 scaleRat b minExp xMin xMax p x =
1139 else if x >= xMax then
1140 scaleRat b minExp xMin xMax (p+1) (x/b)
1141 else if x < xMin then
1142 scaleRat b minExp xMin xMax (p-1) (x*b)
1146 -- Exponentiation with a cache for the most common numbers.
1149 expt :: Integer -> Int -> Integer
1151 if base == 2 && n >= minExpt && n <= maxExpt then
1155 expts :: Array Int Integer
1156 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1158 -- Compute the (floor of the) log of i in base b.
1159 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1160 -- be very slow! We are just slightly more clever.
1161 integerLogBase :: Integer -> Integer -> Int
1162 integerLogBase b i =
1166 -- Try squaring the base first to cut down the number of divisions.
1167 let l = 2 * integerLogBase (b*b) i
1168 doDiv :: Integer -> Int -> Int
1169 doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
1170 in doDiv (i `div` (b^l)) l
1174 %*********************************************************
1176 \subsection{Numeric primops}
1178 %*********************************************************
1180 Definitions of the boxed PrimOps; these will be
1181 used in the case of partial applications, etc.
1184 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
1185 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
1186 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
1187 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1188 negateFloat (F# x) = F# (negateFloat# x)
1190 gtFloat (F# x) (F# y) = gtFloat# x y
1191 geFloat (F# x) (F# y) = geFloat# x y
1192 eqFloat (F# x) (F# y) = eqFloat# x y
1193 neFloat (F# x) (F# y) = neFloat# x y
1194 ltFloat (F# x) (F# y) = ltFloat# x y
1195 leFloat (F# x) (F# y) = leFloat# x y
1197 float2Int (F# x) = I# (float2Int# x)
1198 int2Float (I# x) = F# (int2Float# x)
1200 expFloat (F# x) = F# (expFloat# x)
1201 logFloat (F# x) = F# (logFloat# x)
1202 sqrtFloat (F# x) = F# (sqrtFloat# x)
1203 sinFloat (F# x) = F# (sinFloat# x)
1204 cosFloat (F# x) = F# (cosFloat# x)
1205 tanFloat (F# x) = F# (tanFloat# x)
1206 asinFloat (F# x) = F# (asinFloat# x)
1207 acosFloat (F# x) = F# (acosFloat# x)
1208 atanFloat (F# x) = F# (atanFloat# x)
1209 sinhFloat (F# x) = F# (sinhFloat# x)
1210 coshFloat (F# x) = F# (coshFloat# x)
1211 tanhFloat (F# x) = F# (tanhFloat# x)
1213 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
1215 -- definitions of the boxed PrimOps; these will be
1216 -- used in the case of partial applications, etc.
1218 plusDouble (D# x) (D# y) = D# (x +## y)
1219 minusDouble (D# x) (D# y) = D# (x -## y)
1220 timesDouble (D# x) (D# y) = D# (x *## y)
1221 divideDouble (D# x) (D# y) = D# (x /## y)
1222 negateDouble (D# x) = D# (negateDouble# x)
1224 gtDouble (D# x) (D# y) = x >## y
1225 geDouble (D# x) (D# y) = x >=## y
1226 eqDouble (D# x) (D# y) = x ==## y
1227 neDouble (D# x) (D# y) = x /=## y
1228 ltDouble (D# x) (D# y) = x <## y
1229 leDouble (D# x) (D# y) = x <=## y
1231 double2Int (D# x) = I# (double2Int# x)
1232 int2Double (I# x) = D# (int2Double# x)
1233 double2Float (D# x) = F# (double2Float# x)
1234 float2Double (F# x) = D# (float2Double# x)
1236 expDouble (D# x) = D# (expDouble# x)
1237 logDouble (D# x) = D# (logDouble# x)
1238 sqrtDouble (D# x) = D# (sqrtDouble# x)
1239 sinDouble (D# x) = D# (sinDouble# x)
1240 cosDouble (D# x) = D# (cosDouble# x)
1241 tanDouble (D# x) = D# (tanDouble# x)
1242 asinDouble (D# x) = D# (asinDouble# x)
1243 acosDouble (D# x) = D# (acosDouble# x)
1244 atanDouble (D# x) = D# (atanDouble# x)
1245 sinhDouble (D# x) = D# (sinhDouble# x)
1246 coshDouble (D# x) = D# (coshDouble# x)
1247 tanhDouble (D# x) = D# (tanhDouble# x)
1249 powerDouble (D# x) (D# y) = D# (x **## y)