2 % (c) The AQUA Project, Glasgow University, 1994-1996
5 \section[PrelNumExtra]{Module @PrelNumExtra@}
8 {-# OPTIONS -fno-implicit-prelude #-}
11 #include "../includes/ieee-flpt.h"
16 module PrelNumExtra where
21 import {-# SOURCE #-} PrelErr ( error )
24 import Maybe ( fromMaybe )
26 import PrelArr ( Array, array, (!) )
27 import PrelIOBase ( unsafePerformIO )
28 import PrelCCall () -- we need the definitions of CCallable and
29 -- CReturnable for the _ccall_s herein.
32 %*********************************************************
34 \subsection{Type @Float@}
36 %*********************************************************
39 instance Eq Float where
40 (F# x) == (F# y) = x `eqFloat#` y
42 instance Ord Float where
43 (F# x) `compare` (F# y) | x `ltFloat#` y = LT
47 (F# x) < (F# y) = x `ltFloat#` y
48 (F# x) <= (F# y) = x `leFloat#` y
49 (F# x) >= (F# y) = x `geFloat#` y
50 (F# x) > (F# y) = x `gtFloat#` y
52 instance Num Float where
53 (+) x y = plusFloat x y
54 (-) x y = minusFloat x y
55 negate x = negateFloat x
56 (*) x y = timesFloat x y
58 | otherwise = negateFloat x
59 signum x | x == 0.0 = 0
61 | otherwise = negate 1
62 fromInteger n = encodeFloat n 0
63 fromInt i = int2Float i
65 instance Real Float where
66 toRational x = (m%1)*(b%1)^^n
67 where (m,n) = decodeFloat x
70 instance Fractional Float where
71 (/) x y = divideFloat x y
72 fromRational x = fromRat x
75 instance Floating Float where
76 pi = 3.141592653589793238
89 (**) x y = powerFloat x y
90 logBase x y = log y / log x
92 asinh x = log (x + sqrt (1.0+x*x))
93 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
94 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
96 instance RealFrac Float where
98 {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
99 {-# SPECIALIZE truncate :: Float -> Int #-}
100 {-# SPECIALIZE round :: Float -> Int #-}
101 {-# SPECIALIZE ceiling :: Float -> Int #-}
102 {-# SPECIALIZE floor :: Float -> Int #-}
104 {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
105 {-# SPECIALIZE truncate :: Float -> Integer #-}
106 {-# SPECIALIZE round :: Float -> Integer #-}
107 {-# SPECIALIZE ceiling :: Float -> Integer #-}
108 {-# SPECIALIZE floor :: Float -> Integer #-}
111 = case (decodeFloat x) of { (m,n) ->
112 let b = floatRadix x in
114 (fromInteger m * fromInteger b ^ n, 0.0)
116 case (quotRem m (b^(negate n))) of { (w,r) ->
117 (fromInteger w, encodeFloat r n)
121 truncate x = case properFraction x of
124 round x = case properFraction x of
126 m = if r < 0.0 then n - 1 else n + 1
127 half_down = abs r - 0.5
129 case (compare half_down 0.0) of
131 EQ -> if even n then n else m
134 ceiling x = case properFraction x of
135 (n,r) -> if r > 0.0 then n + 1 else n
137 floor x = case properFraction x of
138 (n,r) -> if r < 0.0 then n - 1 else n
140 foreign import ccall "__encodeFloat" unsafe
141 encodeFloat# :: Int# -> ByteArray# -> Int -> Float
142 foreign import ccall "__int_encodeFloat" unsafe
143 int_encodeFloat# :: Int# -> Int -> Float
145 foreign import ccall "isFloatNaN" unsafe isFloatNaN :: Float -> Int
146 foreign import ccall "isFloatInfinite" unsafe isFloatInfinite :: Float -> Int
147 foreign import ccall "isFloatDenormalized" unsafe isFloatDenormalized :: Float -> Int
148 foreign import ccall "isFloatNegativeZero" unsafe isFloatNegativeZero :: Float -> Int
150 instance RealFloat Float where
151 floatRadix _ = FLT_RADIX -- from float.h
152 floatDigits _ = FLT_MANT_DIG -- ditto
153 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
156 = case decodeFloat# f# of
157 (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
159 encodeFloat (S# i) j = int_encodeFloat# i j
160 encodeFloat (J# s# d#) e = encodeFloat# s# d# e
162 exponent x = case decodeFloat x of
163 (m,n) -> if m == 0 then 0 else n + floatDigits x
165 significand x = case decodeFloat x of
166 (m,_) -> encodeFloat m (negate (floatDigits x))
168 scaleFloat k x = case decodeFloat x of
169 (m,n) -> encodeFloat m (n+k)
170 isNaN x = 0 /= isFloatNaN x
171 isInfinite x = 0 /= isFloatInfinite x
172 isDenormalized x = 0 /= isFloatDenormalized x
173 isNegativeZero x = 0 /= isFloatNegativeZero x
177 %*********************************************************
179 \subsection{Type @Double@}
181 %*********************************************************
184 instance Show Float where
185 showsPrec x = showSigned showFloat x
186 showList = showList__ (showsPrec 0)
188 instance Eq Double where
189 (D# x) == (D# y) = x ==## y
191 instance Ord Double where
192 (D# x) `compare` (D# y) | x <## y = LT
196 (D# x) < (D# y) = x <## y
197 (D# x) <= (D# y) = x <=## y
198 (D# x) >= (D# y) = x >=## y
199 (D# x) > (D# y) = x >## y
201 instance Num Double where
202 (+) x y = plusDouble x y
203 (-) x y = minusDouble x y
204 negate x = negateDouble x
205 (*) x y = timesDouble x y
207 | otherwise = negateDouble x
208 signum x | x == 0.0 = 0
210 | otherwise = negate 1
211 fromInteger n = encodeFloat n 0
212 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
214 instance Real Double where
215 toRational x = (m%1)*(b%1)^^n
216 where (m,n) = decodeFloat x
219 instance Fractional Double where
220 (/) x y = divideDouble x y
221 fromRational x = fromRat x
224 instance Floating Double where
225 pi = 3.141592653589793238
228 sqrt x = sqrtDouble x
232 asin x = asinDouble x
233 acos x = acosDouble x
234 atan x = atanDouble x
235 sinh x = sinhDouble x
236 cosh x = coshDouble x
237 tanh x = tanhDouble x
238 (**) x y = powerDouble x y
239 logBase x y = log y / log x
241 asinh x = log (x + sqrt (1.0+x*x))
242 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
243 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
245 instance RealFrac Double where
247 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
248 {-# SPECIALIZE truncate :: Double -> Int #-}
249 {-# SPECIALIZE round :: Double -> Int #-}
250 {-# SPECIALIZE ceiling :: Double -> Int #-}
251 {-# SPECIALIZE floor :: Double -> Int #-}
253 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
254 {-# SPECIALIZE truncate :: Double -> Integer #-}
255 {-# SPECIALIZE round :: Double -> Integer #-}
256 {-# SPECIALIZE ceiling :: Double -> Integer #-}
257 {-# SPECIALIZE floor :: Double -> Integer #-}
259 #if defined(__UNBOXED_INSTANCES__)
260 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
261 {-# SPECIALIZE truncate :: Double -> Int# #-}
262 {-# SPECIALIZE round :: Double -> Int# #-}
263 {-# SPECIALIZE ceiling :: Double -> Int# #-}
264 {-# SPECIALIZE floor :: Double -> Int# #-}
268 = case (decodeFloat x) of { (m,n) ->
269 let b = floatRadix x in
271 (fromInteger m * fromInteger b ^ n, 0.0)
273 case (quotRem m (b^(negate n))) of { (w,r) ->
274 (fromInteger w, encodeFloat r n)
278 truncate x = case properFraction x of
281 round x = case properFraction x of
283 m = if r < 0.0 then n - 1 else n + 1
284 half_down = abs r - 0.5
286 case (compare half_down 0.0) of
288 EQ -> if even n then n else m
291 ceiling x = case properFraction x of
292 (n,r) -> if r > 0.0 then n + 1 else n
294 floor x = case properFraction x of
295 (n,r) -> if r < 0.0 then n - 1 else n
297 foreign import ccall "__encodeDouble" unsafe
298 encodeDouble# :: Int# -> ByteArray# -> Int -> Double
299 foreign import ccall "__int_encodeDouble" unsafe
300 int_encodeDouble# :: Int# -> Int -> Double
302 foreign import ccall "isDoubleNaN" unsafe isDoubleNaN :: Double -> Int
303 foreign import ccall "isDoubleInfinite" unsafe isDoubleInfinite :: Double -> Int
304 foreign import ccall "isDoubleDenormalized" unsafe isDoubleDenormalized :: Double -> Int
305 foreign import ccall "isDoubleNegativeZero" unsafe isDoubleNegativeZero :: Double -> Int
307 instance RealFloat Double where
308 floatRadix _ = FLT_RADIX -- from float.h
309 floatDigits _ = DBL_MANT_DIG -- ditto
310 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
313 = case decodeDouble# x# of
314 (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
316 encodeFloat (S# i) j = int_encodeDouble# i j
317 encodeFloat (J# s# d#) e = encodeDouble# s# d# e
319 exponent x = case decodeFloat x of
320 (m,n) -> if m == 0 then 0 else n + floatDigits x
322 significand x = case decodeFloat x of
323 (m,_) -> encodeFloat m (negate (floatDigits x))
325 scaleFloat k x = case decodeFloat x of
326 (m,n) -> encodeFloat m (n+k)
328 isNaN x = 0 /= isDoubleNaN x
329 isInfinite x = 0 /= isDoubleInfinite x
330 isDenormalized x = 0 /= isDoubleDenormalized x
331 isNegativeZero x = 0 /= isDoubleNegativeZero x
334 instance Show Double where
335 showsPrec x = showSigned showFloat x
336 showList = showList__ (showsPrec 0)
339 %*********************************************************
341 \subsection{Coercions}
343 %*********************************************************
346 {-# SPECIALIZE fromIntegral ::
356 Integer -> Double #-}
357 fromIntegral :: (Integral a, Num b) => a -> b
358 fromIntegral = fromInteger . toInteger
360 {-# SPECIALIZE realToFrac ::
365 Rational -> Rational,
370 realToFrac :: (Real a, Fractional b) => a -> b
371 realToFrac = fromRational . toRational
374 %*********************************************************
376 \subsection{Common code for @Float@ and @Double@}
378 %*********************************************************
380 The @Enum@ instances for Floats and Doubles are slightly unusual.
381 The @toEnum@ function truncates numbers to Int. The definitions
382 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
383 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
384 dubious. This example may have either 10 or 11 elements, depending on
385 how 0.1 is represented.
387 NOTE: The instances for Float and Double do not make use of the default
388 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
389 a `non-lossy' conversion to and from Ints. Instead we make use of the
390 1.2 default methods (back in the days when Enum had Ord as a superclass)
391 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
394 instance Enum Float where
397 toEnum = fromIntegral
398 fromEnum = fromInteger . truncate -- may overflow
399 enumFrom = numericEnumFrom
400 enumFromTo = numericEnumFromTo
401 enumFromThen = numericEnumFromThen
402 enumFromThenTo = numericEnumFromThenTo
404 instance Enum Double where
407 toEnum = fromIntegral
408 fromEnum = fromInteger . truncate -- may overflow
409 enumFrom = numericEnumFrom
410 enumFromTo = numericEnumFromTo
411 enumFromThen = numericEnumFromThen
412 enumFromThenTo = numericEnumFromThenTo
414 numericEnumFrom :: (Fractional a) => a -> [a]
415 numericEnumFrom = iterate (+1)
417 numericEnumFromThen :: (Fractional a) => a -> a -> [a]
418 numericEnumFromThen n m = iterate (+(m-n)) n
420 numericEnumFromTo :: (Ord a, Fractional a) => a -> a -> [a]
421 numericEnumFromTo n m = takeWhile (<= m + 1/2) (numericEnumFrom n)
423 numericEnumFromThenTo :: (Ord a, Fractional a) => a -> a -> a -> [a]
424 numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
427 pred | e2 > e1 = (<= e3 + mid)
428 | otherwise = (>= e3 + mid)
432 @approxRational@, applied to two real fractional numbers x and epsilon,
433 returns the simplest rational number within epsilon of x. A rational
434 number n%d in reduced form is said to be simpler than another n'%d' if
435 abs n <= abs n' && d <= d'. Any real interval contains a unique
436 simplest rational; here, for simplicity, we assume a closed rational
437 interval. If such an interval includes at least one whole number, then
438 the simplest rational is the absolutely least whole number. Otherwise,
439 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
440 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
441 the simplest rational between d'%r' and d%r.
444 approxRational :: (RealFrac a) => a -> a -> Rational
445 approxRational rat eps = simplest (rat-eps) (rat+eps)
446 where simplest x y | y < x = simplest y x
448 | x > 0 = simplest' n d n' d'
449 | y < 0 = - simplest' (-n') d' (-n) d
451 where xr = toRational x
458 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
460 | q /= q' = (q+1) :% 1
461 | otherwise = (q*n''+d'') :% n''
462 where (q,r) = quotRem n d
463 (q',r') = quotRem n' d'
464 nd'' = simplest' d' r' d r
466 d'' = denominator nd''
471 instance (Integral a) => Ord (Ratio a) where
472 (x:%y) <= (x':%y') = x * y' <= x' * y
473 (x:%y) < (x':%y') = x * y' < x' * y
475 instance (Integral a) => Num (Ratio a) where
476 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
477 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
478 (x:%y) * (x':%y') = reduce (x * x') (y * y')
479 negate (x:%y) = (-x) :% y
480 abs (x:%y) = abs x :% y
481 signum (x:%_) = signum x :% 1
482 fromInteger x = fromInteger x :% 1
484 instance (Integral a) => Real (Ratio a) where
485 toRational (x:%y) = toInteger x :% toInteger y
487 instance (Integral a) => Fractional (Ratio a) where
488 (x:%y) / (x':%y') = (x*y') % (y*x')
489 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
490 fromRational (x:%y) = fromInteger x :% fromInteger y
492 instance (Integral a) => RealFrac (Ratio a) where
493 properFraction (x:%y) = (fromIntegral q, r:%y)
494 where (q,r) = quotRem x y
496 instance (Integral a) => Enum (Ratio a) where
500 toEnum n = fromIntegral n :% 1
501 fromEnum = fromInteger . truncate
503 enumFrom = bounded_iterator True (1)
504 enumFromThen n m = bounded_iterator (diff >= 0) diff n
508 bounded_iterator :: (Ord a, Num a) => Bool -> a -> a -> [a]
509 bounded_iterator inc step v
510 | inc && v > new_v = [v] -- oflow
511 | not inc && v < new_v = [v] -- uflow
512 | otherwise = v : bounded_iterator inc step new_v
519 instance (Integral a) => Show (Ratio a) where
520 showsPrec p (x:%y) = showParen (p > ratio_prec)
521 (shows x . showString " % " . shows y)
524 @showRational@ converts a Rational to a string that looks like a
525 floating point number, but without converting to any floating type
526 (because of the possible overflow).
528 From/by Lennart, 94/09/26
531 showRational :: Int -> Rational -> String
536 let (r', e) = normalize r
542 -- make sure 1 <= r < 10
543 normalize :: Rational -> (Rational, Int)
544 normalize r = if r < 1 then
545 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
548 where norm :: Int -> Rational -> Int -> (Rational, Int)
549 -- Invariant: x*10^e == original r
554 in if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
556 prR :: Int -> Rational -> Int -> String
557 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
558 prR n r e | r >= 10 = prR n (r/10) (e+1)
560 | e > 0 && e < 8 = takeN e s ('.' : drop0 (drop e s) [])
561 | e <= 0 && e > -3 = '0': '.' : takeN (-e) (repeat '0') (drop0 s [])
562 | otherwise = h : '.' : drop0 t ('e':show e0)
564 s@(h:t) = show ((round (r * 10^n))::Integer)
567 takeN (I# n#) ls rs = takeUInt_append n# ls rs
569 drop0 :: String -> String -> String
571 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
573 dropTrailing0s [] = Nothing
574 dropTrailing0s ('0':xs) =
575 case dropTrailing0s xs of
577 Just ls -> Just ('0':ls)
578 dropTrailing0s (x:xs) =
579 case dropTrailing0s xs of
581 Just ls -> Just (x:ls)
585 [In response to a request for documentation of how fromRational works,
586 Joe Fasel writes:] A quite reasonable request! This code was added to
587 the Prelude just before the 1.2 release, when Lennart, working with an
588 early version of hbi, noticed that (read . show) was not the identity
589 for floating-point numbers. (There was a one-bit error about half the
590 time.) The original version of the conversion function was in fact
591 simply a floating-point divide, as you suggest above. The new version
592 is, I grant you, somewhat denser.
594 Unfortunately, Joe's code doesn't work! Here's an example:
596 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
601 1.8217369128763981e-300
603 Lennart's code follows, and it works...
606 fromRat :: (RealFloat a) => Rational -> a
610 -- If the exponent of the nearest floating-point number to x
611 -- is e, then the significand is the integer nearest xb^(-e),
612 -- where b is the floating-point radix. We start with a good
613 -- guess for e, and if it is correct, the exponent of the
614 -- floating-point number we construct will again be e. If
615 -- not, one more iteration is needed.
617 f e = if e' == e then y else f e'
618 where y = encodeFloat (round (x * (1 % b)^^e)) e
619 (_,e') = decodeFloat y
622 -- We obtain a trial exponent by doing a floating-point
623 -- division of x's numerator by its denominator. The
624 -- result of this division may not itself be the ultimate
625 -- result, because of an accumulation of three rounding
628 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
629 / fromInteger (denominator x))
632 Now, here's Lennart's code.
635 {-# SPECIALISE fromRat ::
637 Rational -> Float #-}
638 fromRat :: (RealFloat a) => Rational -> a
640 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
641 | x < 0 = - fromRat' (-x) -- first.
642 | otherwise = fromRat' x
644 -- Conversion process:
645 -- Scale the rational number by the RealFloat base until
646 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
647 -- Then round the rational to an Integer and encode it with the exponent
648 -- that we got from the scaling.
649 -- To speed up the scaling process we compute the log2 of the number to get
650 -- a first guess of the exponent.
652 fromRat' :: (RealFloat a) => Rational -> a
654 where b = floatRadix r
656 (minExp0, _) = floatRange r
657 minExp = minExp0 - p -- the real minimum exponent
658 xMin = toRational (expt b (p-1))
659 xMax = toRational (expt b p)
660 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
661 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
662 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
663 r = encodeFloat (round x') p'
665 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
666 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
667 scaleRat b minExp xMin xMax p x
668 | p <= minExp = (x, p)
669 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
670 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
673 -- Exponentiation with a cache for the most common numbers.
674 minExpt, maxExpt :: Int
678 expt :: Integer -> Int -> Integer
680 if base == 2 && n >= minExpt && n <= maxExpt then
685 expts :: Array Int Integer
686 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
688 -- Compute the (floor of the) log of i in base b.
689 -- Simplest way would be just divide i by b until it's smaller then b, but that would
690 -- be very slow! We are just slightly more clever.
691 integerLogBase :: Integer -> Integer -> Int
694 | otherwise = doDiv (i `div` (b^l)) l
696 -- Try squaring the base first to cut down the number of divisions.
697 l = 2 * integerLogBase (b*b) i
699 doDiv :: Integer -> Int -> Int
702 | otherwise = doDiv (x `div` b) (y+1)
706 %*********************************************************
708 \subsection{Printing out numbers}
710 %*********************************************************
713 --Exported from std library Numeric, defined here to
714 --avoid mut. rec. between PrelNum and Numeric.
715 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
716 showSigned showPos p x
717 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
718 | otherwise = showPos x
720 showFloat :: (RealFloat a) => a -> ShowS
721 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
723 -- These are the format types. This type is not exported.
725 data FFFormat = FFExponent | FFFixed | FFGeneric
727 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
728 formatRealFloat fmt decs x
730 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
731 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
732 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
736 doFmt format (is, e) =
737 let ds = map intToDigit is in
740 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
745 let e' = if e==0 then 0 else e-1 in
748 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
750 let dec' = max dec 1 in
752 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
755 (ei,is') = roundTo base (dec'+1) is
756 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
758 d:'.':ds' ++ 'e':show (e-1+ei)
761 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
766 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
767 f n s "" = f (n-1) ('0':s) ""
768 f n s (r:rs) = f (n-1) (r:s) rs
772 let dec' = max dec 0 in
775 (ei,is') = roundTo base (dec' + e) is
776 (ls,rs) = splitAt (e+ei) (map intToDigit is')
778 mk0 ls ++ (if null rs then "" else '.':rs)
781 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
782 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
787 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
795 f n [] = (0, replicate n 0)
796 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
798 | i' == base = (1,0:ds)
799 | otherwise = (0,i':ds)
805 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
806 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
807 -- This version uses a much slower logarithm estimator. It should be improved.
809 -- This function returns a list of digits (Ints in [0..base-1]) and an
812 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
813 floatToDigits _ 0 = ([0], 0)
814 floatToDigits base x =
816 (f0, e0) = decodeFloat x
817 (minExp0, _) = floatRange x
820 minExp = minExp0 - p -- the real minimum exponent
821 -- Haskell requires that f be adjusted so denormalized numbers
822 -- will have an impossibly low exponent. Adjust for this.
824 let n = minExp - e0 in
825 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
830 (f*be*b*2, 2*b, be*b, b)
834 if e > minExp && f == b^(p-1) then
835 (f*b*2, b^(-e+1)*2, b, 1)
837 (f*2, b^(-e)*2, 1, 1)
841 if b == 2 && base == 10 then
842 -- logBase 10 2 is slightly bigger than 3/10 so
843 -- the following will err on the low side. Ignoring
844 -- the fraction will make it err even more.
845 -- Haskell promises that p-1 <= logBase b f < p.
846 (p - 1 + e0) * 3 `div` 10
848 ceiling ((log (fromInteger (f+1)) +
849 fromInt e * log (fromInteger b)) /
850 log (fromInteger base))
851 --WAS: fromInt e * log (fromInteger b))
855 if r + mUp <= expt base n * s then n else fixup (n+1)
857 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
861 gen ds rn sN mUpN mDnN =
863 (dn, rn') = (rn * base) `divMod` sN
867 case (rn' < mDnN', rn' + mUpN' > sN) of
868 (True, False) -> dn : ds
869 (False, True) -> dn+1 : ds
870 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
871 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
875 gen [] r (s * expt base k) mUp mDn
877 let bk = expt base (-k) in
878 gen [] (r * bk) s (mUp * bk) (mDn * bk)
880 (map toInt (reverse rds), k)
884 %*********************************************************
886 \subsection{Numeric primops}
888 %*********************************************************
890 Definitions of the boxed PrimOps; these will be
891 used in the case of partial applications, etc.
894 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
895 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
896 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
897 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
898 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
900 negateFloat :: Float -> Float
901 negateFloat (F# x) = F# (negateFloat# x)
903 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
904 gtFloat (F# x) (F# y) = gtFloat# x y
905 geFloat (F# x) (F# y) = geFloat# x y
906 eqFloat (F# x) (F# y) = eqFloat# x y
907 neFloat (F# x) (F# y) = neFloat# x y
908 ltFloat (F# x) (F# y) = ltFloat# x y
909 leFloat (F# x) (F# y) = leFloat# x y
911 float2Int :: Float -> Int
912 float2Int (F# x) = I# (float2Int# x)
914 int2Float :: Int -> Float
915 int2Float (I# x) = F# (int2Float# x)
917 expFloat, logFloat, sqrtFloat :: Float -> Float
918 sinFloat, cosFloat, tanFloat :: Float -> Float
919 asinFloat, acosFloat, atanFloat :: Float -> Float
920 sinhFloat, coshFloat, tanhFloat :: Float -> Float
921 expFloat (F# x) = F# (expFloat# x)
922 logFloat (F# x) = F# (logFloat# x)
923 sqrtFloat (F# x) = F# (sqrtFloat# x)
924 sinFloat (F# x) = F# (sinFloat# x)
925 cosFloat (F# x) = F# (cosFloat# x)
926 tanFloat (F# x) = F# (tanFloat# x)
927 asinFloat (F# x) = F# (asinFloat# x)
928 acosFloat (F# x) = F# (acosFloat# x)
929 atanFloat (F# x) = F# (atanFloat# x)
930 sinhFloat (F# x) = F# (sinhFloat# x)
931 coshFloat (F# x) = F# (coshFloat# x)
932 tanhFloat (F# x) = F# (tanhFloat# x)
934 powerFloat :: Float -> Float -> Float
935 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
937 -- definitions of the boxed PrimOps; these will be
938 -- used in the case of partial applications, etc.
940 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
941 plusDouble (D# x) (D# y) = D# (x +## y)
942 minusDouble (D# x) (D# y) = D# (x -## y)
943 timesDouble (D# x) (D# y) = D# (x *## y)
944 divideDouble (D# x) (D# y) = D# (x /## y)
946 negateDouble :: Double -> Double
947 negateDouble (D# x) = D# (negateDouble# x)
949 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
950 gtDouble (D# x) (D# y) = x >## y
951 geDouble (D# x) (D# y) = x >=## y
952 eqDouble (D# x) (D# y) = x ==## y
953 neDouble (D# x) (D# y) = x /=## y
954 ltDouble (D# x) (D# y) = x <## y
955 leDouble (D# x) (D# y) = x <=## y
957 double2Int :: Double -> Int
958 double2Int (D# x) = I# (double2Int# x)
960 int2Double :: Int -> Double
961 int2Double (I# x) = D# (int2Double# x)
963 double2Float :: Double -> Float
964 double2Float (D# x) = F# (double2Float# x)
965 float2Double :: Float -> Double
966 float2Double (F# x) = D# (float2Double# x)
968 expDouble, logDouble, sqrtDouble :: Double -> Double
969 sinDouble, cosDouble, tanDouble :: Double -> Double
970 asinDouble, acosDouble, atanDouble :: Double -> Double
971 sinhDouble, coshDouble, tanhDouble :: Double -> Double
972 expDouble (D# x) = D# (expDouble# x)
973 logDouble (D# x) = D# (logDouble# x)
974 sqrtDouble (D# x) = D# (sqrtDouble# x)
975 sinDouble (D# x) = D# (sinDouble# x)
976 cosDouble (D# x) = D# (cosDouble# x)
977 tanDouble (D# x) = D# (tanDouble# x)
978 asinDouble (D# x) = D# (asinDouble# x)
979 acosDouble (D# x) = D# (acosDouble# x)
980 atanDouble (D# x) = D# (atanDouble# x)
981 sinhDouble (D# x) = D# (sinhDouble# x)
982 coshDouble (D# x) = D# (coshDouble# x)
983 tanhDouble (D# x) = D# (tanhDouble# x)
985 powerDouble :: Double -> Double -> Double
986 powerDouble (D# x) (D# y) = D# (x **## y)