2 % (c) The AQUA Project, Glasgow University, 1994-1996
5 \section[PrelNum]{Module @PrelNum@}
7 Numeric part of the prelude.
13 #include "../includes/ieee-flpt.h"
17 {-# OPTIONS -fno-implicit-prelude #-}
21 import {-# SOURCE #-} IOBase ( error )
24 import ArrBase ( Array, array, (!) )
29 infixl 7 /, %, `quot`, `rem`, `div`, `mod`
33 %*********************************************************
35 \subsection{Standard numeric classes}
37 %*********************************************************
40 class (Num a, Ord a) => Real a where
41 toRational :: a -> Rational
43 class (Real a, Enum a) => Integral a where
44 quot, rem, div, mod :: a -> a -> a
45 quotRem, divMod :: a -> a -> (a,a)
46 toInteger :: a -> Integer
47 toInt :: a -> Int -- partain: Glasgow extension
49 n `quot` d = q where (q,r) = quotRem n d
50 n `rem` d = r where (q,r) = quotRem n d
51 n `div` d = q where (q,r) = divMod n d
52 n `mod` d = r where (q,r) = divMod n d
53 divMod n d = if signum r == negate (signum d) then (q-1, r+d) else qr
54 where qr@(q,r) = quotRem n d
56 class (Num a) => Fractional a where
59 fromRational :: Rational -> a
63 class (Fractional a) => Floating a where
65 exp, log, sqrt :: a -> a
66 (**), logBase :: a -> a -> a
67 sin, cos, tan :: a -> a
68 asin, acos, atan :: a -> a
69 sinh, cosh, tanh :: a -> a
70 asinh, acosh, atanh :: a -> a
72 x ** y = exp (log x * y)
73 logBase x y = log y / log x
76 tanh x = sinh x / cosh x
78 class (Real a, Fractional a) => RealFrac a where
79 properFraction :: (Integral b) => a -> (b,a)
80 truncate, round :: (Integral b) => a -> b
81 ceiling, floor :: (Integral b) => a -> b
83 truncate x = m where (m,_) = properFraction x
85 round x = let (n,r) = properFraction x
86 m = if r < 0 then n - 1 else n + 1
87 in case signum (abs r - 0.5) of
89 0 -> if even n then n else m
92 ceiling x = if r > 0 then n + 1 else n
93 where (n,r) = properFraction x
95 floor x = if r < 0 then n - 1 else n
96 where (n,r) = properFraction x
98 class (RealFrac a, Floating a) => RealFloat a where
99 floatRadix :: a -> Integer
100 floatDigits :: a -> Int
101 floatRange :: a -> (Int,Int)
102 decodeFloat :: a -> (Integer,Int)
103 encodeFloat :: Integer -> Int -> a
105 significand :: a -> a
106 scaleFloat :: Int -> a -> a
107 isNaN, isInfinite, isDenormalized, isNegativeZero, isIEEE
110 exponent x = if m == 0 then 0 else n + floatDigits x
111 where (m,n) = decodeFloat x
113 significand x = encodeFloat m (negate (floatDigits x))
114 where (m,_) = decodeFloat x
116 scaleFloat k x = encodeFloat m (n+k)
117 where (m,n) = decodeFloat x
120 %*********************************************************
122 \subsection{Overloaded numeric functions}
124 %*********************************************************
127 even, odd :: (Integral a) => a -> Bool
128 even n = n `rem` 2 == 0
131 {-# GENERATE_SPECS gcd a{Int#,Int,Integer} #-}
132 gcd :: (Integral a) => a -> a -> a
133 gcd 0 0 = error "Prelude.gcd: gcd 0 0 is undefined"
134 gcd x y = gcd' (abs x) (abs y)
136 gcd' x y = gcd' y (x `rem` y)
138 {-# GENERATE_SPECS lcm a{Int#,Int,Integer} #-}
139 lcm :: (Integral a) => a -> a -> a
142 lcm x y = abs ((x `quot` (gcd x y)) * y)
144 (^) :: (Num a, Integral b) => a -> b -> a
146 x ^ n | n > 0 = f x (n-1) x
148 f x n y = g x n where
149 g x n | even n = g (x*x) (n `quot` 2)
150 | otherwise = f x (n-1) (x*y)
151 _ ^ _ = error "Prelude.^: negative exponent"
153 (^^) :: (Fractional a, Integral b) => a -> b -> a
154 x ^^ n = if n >= 0 then x^n else recip (x^(negate n))
156 fromIntegral :: (Integral a, Num b) => a -> b
157 fromIntegral = fromInteger . toInteger
159 fromRealFrac :: (RealFrac a, Fractional b) => a -> b
160 fromRealFrac = fromRational . toRational
162 atan2 :: (RealFloat a) => a -> a -> a
163 atan2 y x = case (signum y, signum x) of
167 (-1, 0) -> (negate pi)/2
168 ( _, 1) -> atan (y/x)
169 ( _,-1) -> atan (y/x) + pi
170 ( 0, 0) -> error "Prelude.atan2: atan2 of origin"
174 %*********************************************************
176 \subsection{Instances for @Int@}
178 %*********************************************************
181 instance Real Int where
182 toRational x = toInteger x % 1
184 instance Integral Int where
185 a@(I# _) `quotRem` b@(I# _) = (a `quotInt` b, a `remInt` b)
186 -- OK, so I made it a little stricter. Shoot me. (WDP 94/10)
188 -- following chks for zero divisor are non-standard (WDP)
189 a `quot` b = if b /= 0
191 else error "Integral.Int.quot{PreludeCore}: divide by 0\n"
192 a `rem` b = if b /= 0
194 else error "Integral.Int.rem{PreludeCore}: divide by 0\n"
196 x `div` y = if x > 0 && y < 0 then quotInt (x-y-1) y
197 else if x < 0 && y > 0 then quotInt (x-y+1) y
199 x `mod` y = if x > 0 && y < 0 || x < 0 && y > 0 then
200 if r/=0 then r+y else 0
205 divMod x@(I# _) y@(I# _) = (x `div` y, x `mod` y)
206 -- Stricter. Sorry if you don't like it. (WDP 94/10)
208 --OLD: even x = eqInt (x `mod` 2) 0
209 --OLD: odd x = neInt (x `mod` 2) 0
211 toInteger (I# n#) = int2Integer# n# -- give back a full-blown Integer
217 %*********************************************************
219 \subsection{Type @Integer@}
221 %*********************************************************
223 These types are used to return from integer primops
226 data Return2GMPs = Return2GMPs Int# Int# ByteArray# Int# Int# ByteArray#
227 data ReturnIntAndGMP = ReturnIntAndGMP Int# Int# Int# ByteArray#
233 instance Eq Integer where
234 (J# a1 s1 d1) == (J# a2 s2 d2)
235 = (cmpInteger# a1 s1 d1 a2 s2 d2) ==# 0#
237 (J# a1 s1 d1) /= (J# a2 s2 d2)
238 = (cmpInteger# a1 s1 d1 a2 s2 d2) /=# 0#
240 instance Ord Integer where
241 (J# a1 s1 d1) <= (J# a2 s2 d2)
242 = (cmpInteger# a1 s1 d1 a2 s2 d2) <=# 0#
244 (J# a1 s1 d1) < (J# a2 s2 d2)
245 = (cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#
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 x@(J# a1 s1 d1) `max` y@(J# a2 s2 d2)
254 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) ># 0#) then x else y
256 x@(J# a1 s1 d1) `min` y@(J# a2 s2 d2)
257 = if ((cmpInteger# a1 s1 d1 a2 s2 d2) <# 0#) then x else y
259 compare (J# a1 s1 d1) (J# a2 s2 d2)
260 = case cmpInteger# a1 s1 d1 a2 s2 d2 of { res# ->
261 if res# <# 0# then LT else
262 if res# ># 0# then GT else EQ
265 instance Num Integer where
266 (+) (J# a1 s1 d1) (J# a2 s2 d2)
267 = plusInteger# a1 s1 d1 a2 s2 d2
269 (-) (J# a1 s1 d1) (J# a2 s2 d2)
270 = minusInteger# a1 s1 d1 a2 s2 d2
272 negate (J# a s d) = negateInteger# a s d
274 (*) (J# a1 s1 d1) (J# a2 s2 d2)
275 = timesInteger# a1 s1 d1 a2 s2 d2
277 -- ORIG: abs n = if n >= 0 then n else -n
280 = case 0 of { J# a2 s2 d2 ->
281 if (cmpInteger# a1 s1 d1 a2 s2 d2) >=# 0#
283 else negateInteger# a1 s1 d1
286 signum n@(J# a1 s1 d1)
287 = case 0 of { J# a2 s2 d2 ->
289 cmp = cmpInteger# a1 s1 d1 a2 s2 d2
292 else if cmp ==# 0# then 0
298 fromInt (I# n#) = int2Integer# n# -- gives back a full-blown Integer
300 instance Real Integer where
303 instance Integral Integer where
304 quotRem (J# a1 s1 d1) (J# a2 s2 d2)
305 = case (quotRemInteger# a1 s1 d1 a2 s2 d2) of
306 Return2GMPs a3 s3 d3 a4 s4 d4
307 -> (J# a3 s3 d3, J# a4 s4 d4)
309 {- USING THE UNDERLYING "GMP" CODE IS DUBIOUS FOR NOW:
311 divMod (J# a1 s1 d1) (J# a2 s2 d2)
312 = case (divModInteger# a1 s1 d1 a2 s2 d2) of
313 Return2GMPs a3 s3 d3 a4 s4 d4
314 -> (J# a3 s3 d3, J# a4 s4 d4)
317 toInt (J# a s d) = case (integer2Int# a s d) of { n# -> I# n# }
319 -- the rest are identical to the report default methods;
320 -- you get slightly better code if you let the compiler
321 -- see them right here:
322 n `quot` d = q where (q,r) = quotRem n d
323 n `rem` d = r where (q,r) = quotRem n d
324 n `div` d = q where (q,r) = divMod n d
325 n `mod` d = r where (q,r) = divMod n d
327 divMod n d = case (quotRem n d) of { qr@(q,r) ->
328 if signum r == negate (signum d) then (q - 1, r+d) else qr }
329 -- Case-ified by WDP 94/10
331 instance Enum Integer where
332 enumFrom n = n : enumFrom (n + 1)
333 enumFromThen m n = en' m (n - m)
334 where en' m n = m : en' (m + n) n
335 enumFromTo n m = takeWhile (<= m) (enumFrom n)
336 enumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
339 instance Show Integer where
340 showsPrec x = showSignedInteger x
341 showList = showList__ (showsPrec 0)
343 instance Ix Integer where
346 | inRange b i = fromInteger (i - m)
347 | otherwise = error "Integer.index: Index out of range."
348 inRange (m,n) i = m <= i && i <= n
350 integer_0, integer_1, integer_2, integer_m1 :: Integer
351 integer_0 = int2Integer# 0#
352 integer_1 = int2Integer# 1#
353 integer_2 = int2Integer# 2#
354 integer_m1 = int2Integer# (negateInt# 1#)
358 %*********************************************************
360 \subsection{Type @Float@}
362 %*********************************************************
365 instance Eq Float where
366 (F# x) == (F# y) = x `eqFloat#` y
368 instance Ord Float where
369 (F# x) `compare` (F# y) | x `ltFloat#` y = LT
370 | x `eqFloat#` y = EQ
373 (F# x) < (F# y) = x `ltFloat#` y
374 (F# x) <= (F# y) = x `leFloat#` y
375 (F# x) >= (F# y) = x `geFloat#` y
376 (F# x) > (F# y) = x `gtFloat#` y
378 instance Num Float where
379 (+) x y = plusFloat x y
380 (-) x y = minusFloat x y
381 negate x = negateFloat x
382 (*) x y = timesFloat x y
384 | otherwise = negateFloat x
385 signum x | x == 0.0 = 0
387 | otherwise = negate 1
388 fromInteger n = encodeFloat n 0
389 fromInt i = int2Float i
391 instance Real Float where
392 toRational x = (m%1)*(b%1)^^n
393 where (m,n) = decodeFloat x
396 instance Fractional Float where
397 (/) x y = divideFloat x y
398 fromRational x = fromRational__ x
401 instance Floating Float where
402 pi = 3.141592653589793238
415 (**) x y = powerFloat x y
416 logBase x y = log y / log x
418 asinh x = log (x + sqrt (1.0+x*x))
419 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
420 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
422 instance RealFrac Float where
424 {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
425 {-# SPECIALIZE truncate :: Float -> Int #-}
426 {-# SPECIALIZE round :: Float -> Int #-}
427 {-# SPECIALIZE ceiling :: Float -> Int #-}
428 {-# SPECIALIZE floor :: Float -> Int #-}
430 {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
431 {-# SPECIALIZE truncate :: Float -> Integer #-}
432 {-# SPECIALIZE round :: Float -> Integer #-}
433 {-# SPECIALIZE ceiling :: Float -> Integer #-}
434 {-# SPECIALIZE floor :: Float -> Integer #-}
437 = case (decodeFloat x) of { (m,n) ->
438 let b = floatRadix x in
440 (fromInteger m * fromInteger b ^ n, 0.0)
442 case (quotRem m (b^(negate n))) of { (w,r) ->
443 (fromInteger w, encodeFloat r n)
447 truncate x = case properFraction x of
450 round x = case properFraction x of
452 m = if r < 0.0 then n - 1 else n + 1
453 half_down = abs r - 0.5
455 case (compare half_down 0.0) of
457 EQ -> if even n then n else m
460 ceiling x = case properFraction x of
461 (n,r) -> if r > 0.0 then n + 1 else n
463 floor x = case properFraction x of
464 (n,r) -> if r < 0.0 then n - 1 else n
466 instance RealFloat Float where
467 floatRadix _ = FLT_RADIX -- from float.h
468 floatDigits _ = FLT_MANT_DIG -- ditto
469 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
472 = case decodeFloat# f# of
473 ReturnIntAndGMP exp# a# s# d# ->
474 (J# a# s# d#, I# exp#)
476 encodeFloat (J# a# s# d#) (I# e#)
477 = case encodeFloat# a# s# d# e# of { flt# -> F# flt# }
479 exponent x = case decodeFloat x of
480 (m,n) -> if m == 0 then 0 else n + floatDigits x
482 significand x = case decodeFloat x of
483 (m,_) -> encodeFloat m (negate (floatDigits x))
485 scaleFloat k x = case decodeFloat x of
486 (m,n) -> encodeFloat m (n+k)
488 instance Show Float where
489 showsPrec x = showSigned showFloat x
490 showList = showList__ (showsPrec 0)
493 %*********************************************************
495 \subsection{Type @Double@}
497 %*********************************************************
500 instance Eq Double where
501 (D# x) == (D# y) = x ==## y
503 instance Ord Double where
504 (D# x) `compare` (D# y) | x <## y = LT
508 (D# x) < (D# y) = x <## y
509 (D# x) <= (D# y) = x <=## y
510 (D# x) >= (D# y) = x >=## y
511 (D# x) > (D# y) = x >## y
513 instance Num Double where
514 (+) x y = plusDouble x y
515 (-) x y = minusDouble x y
516 negate x = negateDouble x
517 (*) x y = timesDouble x y
519 | otherwise = negateDouble x
520 signum x | x == 0.0 = 0
522 | otherwise = negate 1
523 fromInteger n = encodeFloat n 0
524 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
526 instance Real Double where
527 toRational x = (m%1)*(b%1)^^n
528 where (m,n) = decodeFloat x
531 instance Fractional Double where
532 (/) x y = divideDouble x y
533 fromRational x = fromRational__ x
536 instance Floating Double where
537 pi = 3.141592653589793238
540 sqrt x = sqrtDouble x
544 asin x = asinDouble x
545 acos x = acosDouble x
546 atan x = atanDouble x
547 sinh x = sinhDouble x
548 cosh x = coshDouble x
549 tanh x = tanhDouble x
550 (**) x y = powerDouble x y
551 logBase x y = log y / log x
553 asinh x = log (x + sqrt (1.0+x*x))
554 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
555 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
557 instance RealFrac Double where
559 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
560 {-# SPECIALIZE truncate :: Double -> Int #-}
561 {-# SPECIALIZE round :: Double -> Int #-}
562 {-# SPECIALIZE ceiling :: Double -> Int #-}
563 {-# SPECIALIZE floor :: Double -> Int #-}
565 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
566 {-# SPECIALIZE truncate :: Double -> Integer #-}
567 {-# SPECIALIZE round :: Double -> Integer #-}
568 {-# SPECIALIZE ceiling :: Double -> Integer #-}
569 {-# SPECIALIZE floor :: Double -> Integer #-}
571 #if defined(__UNBOXED_INSTANCES__)
572 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
573 {-# SPECIALIZE truncate :: Double -> Int# #-}
574 {-# SPECIALIZE round :: Double -> Int# #-}
575 {-# SPECIALIZE ceiling :: Double -> Int# #-}
576 {-# SPECIALIZE floor :: Double -> Int# #-}
580 = case (decodeFloat x) of { (m,n) ->
581 let b = floatRadix x in
583 (fromInteger m * fromInteger b ^ n, 0.0)
585 case (quotRem m (b^(negate n))) of { (w,r) ->
586 (fromInteger w, encodeFloat r n)
590 truncate x = case properFraction x of
593 round x = case properFraction x of
595 m = if r < 0.0 then n - 1 else n + 1
596 half_down = abs r - 0.5
598 case (compare half_down 0.0) of
600 EQ -> if even n then n else m
603 ceiling x = case properFraction x of
604 (n,r) -> if r > 0.0 then n + 1 else n
606 floor x = case properFraction x of
607 (n,r) -> if r < 0.0 then n - 1 else n
609 instance RealFloat Double where
610 floatRadix _ = FLT_RADIX -- from float.h
611 floatDigits _ = DBL_MANT_DIG -- ditto
612 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
615 = case decodeDouble# d# of
616 ReturnIntAndGMP exp# a# s# d# ->
617 (J# a# s# d#, I# exp#)
619 encodeFloat (J# a# s# d#) (I# e#)
620 = case encodeDouble# a# s# d# e# of { dbl# -> D# dbl# }
622 exponent x = case decodeFloat x of
623 (m,n) -> if m == 0 then 0 else n + floatDigits x
625 significand x = case decodeFloat x of
626 (m,_) -> encodeFloat m (negate (floatDigits x))
628 scaleFloat k x = case decodeFloat x of
629 (m,n) -> encodeFloat m (n+k)
631 instance Show Double where
632 showsPrec x = showSigned showFloat x
633 showList = showList__ (showsPrec 0)
637 %*********************************************************
639 \subsection{Common code for @Float@ and @Double@}
641 %*********************************************************
643 The Enum instances for Floats and Doubles are slightly unusual.
644 The `toEnum' function truncates numbers to Int. The definitions
645 of enumFrom and enumFromThen allow floats to be used in arithmetic
646 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
647 dubious. This example may have either 10 or 11 elements, depending on
648 how 0.1 is represented.
651 instance Enum Float where
652 toEnum = fromIntegral
653 fromEnum = fromInteger . truncate -- may overflow
654 enumFrom = numericEnumFrom
655 enumFromThen = numericEnumFromThen
657 instance Enum Double where
658 toEnum = fromIntegral
659 fromEnum = fromInteger . truncate -- may overflow
660 enumFrom = numericEnumFrom
661 enumFromThen = numericEnumFromThen
663 numericEnumFrom :: (Real a) => a -> [a]
664 numericEnumFromThen :: (Real a) => a -> a -> [a]
665 numericEnumFrom = iterate (+1)
666 numericEnumFromThen n m = iterate (+(m-n)) n
670 %*********************************************************
672 \subsection{The @Ratio@ and @Rational@ types}
674 %*********************************************************
677 data (Integral a) => Ratio a = !a :% !a deriving (Eq)
678 type Rational = Ratio Integer
682 (%) :: (Integral a) => a -> a -> Ratio a
683 numerator, denominator :: (Integral a) => Ratio a -> a
684 approxRational :: (RealFrac a) => a -> a -> Rational
688 \tr{reduce} is a subsidiary function used only in this module .
689 It normalises a ratio by dividing both numerator and denominator by
690 their greatest common divisor.
693 reduce _ 0 = error "{Ratio.%}: zero denominator"
694 reduce x y = (x `quot` d) :% (y `quot` d)
699 x % y = reduce (x * signum y) (abs y)
703 denominator (x:%y) = y
707 @approxRational@, applied to two real fractional numbers x and epsilon,
708 returns the simplest rational number within epsilon of x. A rational
709 number n%d in reduced form is said to be simpler than another n'%d' if
710 abs n <= abs n' && d <= d'. Any real interval contains a unique
711 simplest rational; here, for simplicity, we assume a closed rational
712 interval. If such an interval includes at least one whole number, then
713 the simplest rational is the absolutely least whole number. Otherwise,
714 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
715 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
716 the simplest rational between d'%r' and d%r.
719 approxRational x eps = simplest (x-eps) (x+eps)
720 where simplest x y | y < x = simplest y x
722 | x > 0 = simplest' n d n' d'
723 | y < 0 = - simplest' (-n') d' (-n) d
725 where xr@(n:%d) = toRational x
726 (n':%d') = toRational y
728 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
730 | q /= q' = (q+1) :% 1
731 | otherwise = (q*n''+d'') :% n''
732 where (q,r) = quotRem n d
733 (q',r') = quotRem n' d'
734 (n'':%d'') = simplest' d' r' d r
739 instance (Integral a) => Ord (Ratio a) where
740 (x:%y) <= (x':%y') = x * y' <= x' * y
741 (x:%y) < (x':%y') = x * y' < x' * y
743 instance (Integral a) => Num (Ratio a) where
744 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
745 (x:%y) * (x':%y') = reduce (x * x') (y * y')
746 negate (x:%y) = (-x) :% y
747 abs (x:%y) = abs x :% y
748 signum (x:%y) = signum x :% 1
749 fromInteger x = fromInteger x :% 1
751 instance (Integral a) => Real (Ratio a) where
752 toRational (x:%y) = toInteger x :% toInteger y
754 instance (Integral a) => Fractional (Ratio a) where
755 (x:%y) / (x':%y') = (x*y') % (y*x')
756 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
757 fromRational (x:%y) = fromInteger x :% fromInteger y
759 instance (Integral a) => RealFrac (Ratio a) where
760 properFraction (x:%y) = (fromIntegral q, r:%y)
761 where (q,r) = quotRem x y
763 instance (Integral a) => Enum (Ratio a) where
764 enumFrom = iterate ((+)1)
765 enumFromThen n m = iterate ((+)(m-n)) n
766 toEnum n = fromIntegral n :% 1
767 fromEnum = fromInteger . truncate
772 instance (Integral a) => Show (Ratio a) where
773 showsPrec p (x:%y) = showParen (p > ratio_prec)
774 (shows x . showString " % " . shows y)
777 [In response to a request for documentation of how fromRational works,
778 Joe Fasel writes:] A quite reasonable request! This code was added to
779 the Prelude just before the 1.2 release, when Lennart, working with an
780 early version of hbi, noticed that (read . show) was not the identity
781 for floating-point numbers. (There was a one-bit error about half the
782 time.) The original version of the conversion function was in fact
783 simply a floating-point divide, as you suggest above. The new version
784 is, I grant you, somewhat denser.
786 Unfortunately, Joe's code doesn't work! Here's an example:
788 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
793 1.8217369128763981e-300
795 Lennart's code follows, and it works...
798 {-# GENERATE_SPECS fromRational__ a{Double#,Double} #-}
799 fromRational__ :: (RealFloat a) => Rational -> a
800 fromRational__ x = x'
803 -- If the exponent of the nearest floating-point number to x
804 -- is e, then the significand is the integer nearest xb^(-e),
805 -- where b is the floating-point radix. We start with a good
806 -- guess for e, and if it is correct, the exponent of the
807 -- floating-point number we construct will again be e. If
808 -- not, one more iteration is needed.
810 f e = if e' == e then y else f e'
811 where y = encodeFloat (round (x * (1 % b)^^e)) e
812 (_,e') = decodeFloat y
815 -- We obtain a trial exponent by doing a floating-point
816 -- division of x's numerator by its denominator. The
817 -- result of this division may not itself be the ultimate
818 -- result, because of an accumulation of three rounding
821 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
822 / fromInteger (denominator x))
825 Now, here's Lennart's code.
828 fromRational__ :: (RealFloat a) => Rational -> a
830 if x == 0 then encodeFloat 0 0 -- Handle exceptional cases
831 else if x < 0 then - fromRat' (-x) -- first.
834 -- Conversion process:
835 -- Scale the rational number by the RealFloat base until
836 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
837 -- Then round the rational to an Integer and encode it with the exponent
838 -- that we got from the scaling.
839 -- To speed up the scaling process we compute the log2 of the number to get
840 -- a first guess of the exponent.
842 fromRat' :: (RealFloat a) => Rational -> a
844 where b = floatRadix r
846 (minExp0, _) = floatRange r
847 minExp = minExp0 - p -- the real minimum exponent
848 xMin = toRational (expt b (p-1))
849 xMax = toRational (expt b p)
850 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
851 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
852 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
853 r = encodeFloat (round x') p'
855 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
856 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
857 scaleRat b minExp xMin xMax p x =
860 else if x >= xMax then
861 scaleRat b minExp xMin xMax (p+1) (x/b)
862 else if x < xMin then
863 scaleRat b minExp xMin xMax (p-1) (x*b)
867 -- Exponentiation with a cache for the most common numbers.
870 expt :: Integer -> Int -> Integer
872 if base == 2 && n >= minExpt && n <= maxExpt then
876 expts :: Array Int Integer
877 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
879 -- Compute the (floor of the) log of i in base b.
880 -- Simplest way would be just divide i by b until it's smaller then b, but that would
881 -- be very slow! We are just slightly more clever.
882 integerLogBase :: Integer -> Integer -> Int
887 -- Try squaring the base first to cut down the number of divisions.
888 let l = 2 * integerLogBase (b*b) i
889 doDiv :: Integer -> Int -> Int
890 doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
891 in doDiv (i `div` (b^l)) l
894 %*********************************************************
896 \subsection{Showing numbers}
898 %*********************************************************
902 = case quotRem n 10 of { (n', d) ->
903 case (chr (ord_0 + fromIntegral d)) of { C# c# -> -- stricter than necessary
907 if n' == 0 then r' else showInteger n' r'
910 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
911 showSigned showPos p x = if x < 0 then showParen (p > 6)
912 (showChar '-' . showPos (-x))
915 showSignedInteger :: Int -> Integer -> ShowS
916 showSignedInteger p n r
917 = -- from HBC version; support code follows
918 if n < 0 && p > 6 then '(':jtos n++(')':r) else jtos n ++ r
920 jtos :: Integer -> String
927 jtos' :: Integer -> String -> String
930 chr (fromInteger (n + ord_0)) : cs
932 jtos' (n `quot` 10) (chr (fromInteger (n `rem` 10 + ord_0)) : cs)
935 The functions showFloat below uses rational arithmetic
936 to insure correct conversion between the floating-point radix and
937 decimal. It is often possible to use a higher-precision floating-
938 point type to obtain the same results.
941 {-# GENERATE_SPECS showFloat a{Double#,Double} #-}
942 showFloat:: (RealFloat a) => a -> ShowS
944 if isNaN x then showString "NaN"
945 else if isInfinite x then
946 ( if x < 0 then showString "-" else id) . showString "Infinite"
947 else if x < 0 || isNegativeZero x then showChar '-' . showFloat' (-x) else showFloat' x
949 showFloat' :: (RealFloat a) => a -> ShowS
951 if x == 0 then showString ("0." ++ take (m-1) zeros)
952 else if e >= m-1 || e < 0 then showSci else showFix
954 showFix = showString whole . showChar '.' . showString frac
955 where (whole,frac) = splitAt (e+1) (show sig)
956 showSci = showChar d . showChar '.' . showString frac
957 . showChar 'e' . shows e
958 where (d:frac) = show sig
959 (m, sig, e) = if b == 10 then (w, s, n+w-1)
962 ((fromInt w * log (fromInteger b)) / log 10 :: Double)
964 (sig', e') = if sig1 >= 10^m' then (round (t/10), e1+1)
965 else if sig1 < 10^(m'-1) then (round (t*10), e1-1)
968 t = s%1 * (b%1)^^n * 10^^(m'-e1-1)
969 e1 = floor (logBase 10 x)
970 (s, n) = decodeFloat x
977 @showRational@ converts a Rational to a string that looks like a
978 floating point number, but without converting to any floating type
979 (because of the possible overflow).
981 From/by Lennart, 94/09/26
984 showRational :: Int -> Rational -> String
989 let (r', e) = normalize r
992 startExpExp = 4 :: Int
994 -- make sure 1 <= r < 10
995 normalize :: Rational -> (Rational, Int)
996 normalize r = if r < 1 then
997 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
1000 where norm :: Int -> Rational -> Int -> (Rational, Int)
1001 -- Invariant: r*10^e == original r
1006 in if r >= tn then norm ee (r/tn) (e+n) else norm (ee-1) r e
1009 drop0 (c:cs) = c : reverse (dropWhile (=='0') (reverse cs))
1011 prR :: Int -> Rational -> Int -> String
1012 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
1013 prR n r e | r >= 10 = prR n (r/10) (e+1)
1015 let s = show ((round (r * 10^n))::Integer)
1017 in if e > 0 && e < 8 then
1018 take e s ++ "." ++ drop0 (drop e s)
1019 else if e <= 0 && e > -3 then
1020 "0." ++ take (-e) (repeat '0') ++ drop0 s
1022 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
1025 %*********************************************************
1027 \subsection{Numeric primops}
1029 %*********************************************************
1031 Definitions of the boxed PrimOps; these will be
1032 used in the case of partial applications, etc.
1035 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
1036 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
1037 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
1038 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1039 negateFloat (F# x) = F# (negateFloat# x)
1041 gtFloat (F# x) (F# y) = gtFloat# x y
1042 geFloat (F# x) (F# y) = geFloat# x y
1043 eqFloat (F# x) (F# y) = eqFloat# x y
1044 neFloat (F# x) (F# y) = neFloat# x y
1045 ltFloat (F# x) (F# y) = ltFloat# x y
1046 leFloat (F# x) (F# y) = leFloat# x y
1048 float2Int (F# x) = I# (float2Int# x)
1049 int2Float (I# x) = F# (int2Float# x)
1051 expFloat (F# x) = F# (expFloat# x)
1052 logFloat (F# x) = F# (logFloat# x)
1053 sqrtFloat (F# x) = F# (sqrtFloat# x)
1054 sinFloat (F# x) = F# (sinFloat# x)
1055 cosFloat (F# x) = F# (cosFloat# x)
1056 tanFloat (F# x) = F# (tanFloat# x)
1057 asinFloat (F# x) = F# (asinFloat# x)
1058 acosFloat (F# x) = F# (acosFloat# x)
1059 atanFloat (F# x) = F# (atanFloat# x)
1060 sinhFloat (F# x) = F# (sinhFloat# x)
1061 coshFloat (F# x) = F# (coshFloat# x)
1062 tanhFloat (F# x) = F# (tanhFloat# x)
1064 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
1066 -- definitions of the boxed PrimOps; these will be
1067 -- used in the case of partial applications, etc.
1069 plusDouble (D# x) (D# y) = D# (x +## y)
1070 minusDouble (D# x) (D# y) = D# (x -## y)
1071 timesDouble (D# x) (D# y) = D# (x *## y)
1072 divideDouble (D# x) (D# y) = D# (x /## y)
1073 negateDouble (D# x) = D# (negateDouble# x)
1075 gtDouble (D# x) (D# y) = x >## y
1076 geDouble (D# x) (D# y) = x >=## y
1077 eqDouble (D# x) (D# y) = x ==## y
1078 neDouble (D# x) (D# y) = x /=## y
1079 ltDouble (D# x) (D# y) = x <## y
1080 leDouble (D# x) (D# y) = x <=## y
1082 double2Int (D# x) = I# (double2Int# x)
1083 int2Double (I# x) = D# (int2Double# x)
1084 double2Float (D# x) = F# (double2Float# x)
1085 float2Double (F# x) = D# (float2Double# x)
1087 expDouble (D# x) = D# (expDouble# x)
1088 logDouble (D# x) = D# (logDouble# x)
1089 sqrtDouble (D# x) = D# (sqrtDouble# x)
1090 sinDouble (D# x) = D# (sinDouble# x)
1091 cosDouble (D# x) = D# (cosDouble# x)
1092 tanDouble (D# x) = D# (tanDouble# x)
1093 asinDouble (D# x) = D# (asinDouble# x)
1094 acosDouble (D# x) = D# (acosDouble# x)
1095 atanDouble (D# x) = D# (atanDouble# x)
1096 sinhDouble (D# x) = D# (sinhDouble# x)
1097 coshDouble (D# x) = D# (coshDouble# x)
1098 tanhDouble (D# x) = D# (tanhDouble# x)
1100 powerDouble (D# x) (D# y) = D# (x **## y)