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 instance RealFloat Float where
141 floatRadix _ = FLT_RADIX -- from float.h
142 floatDigits _ = FLT_MANT_DIG -- ditto
143 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
146 = case decodeFloat# f# of
147 (# exp#, a#, s#, d# #) -> (J# a# s# d#, I# exp#)
149 encodeFloat (J# a# s# d#) (I# e#)
150 = case encodeFloat# a# s# d# e# of { flt# -> F# flt# }
152 exponent x = case decodeFloat x of
153 (m,n) -> if m == 0 then 0 else n + floatDigits x
155 significand x = case decodeFloat x of
156 (m,_) -> encodeFloat m (negate (floatDigits x))
158 scaleFloat k x = case decodeFloat x of
159 (m,n) -> encodeFloat m (n+k)
161 (0::Int) /= unsafePerformIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
163 (0::Int) /= unsafePerformIO (_ccall_ isFloatInfinite x) {- ditto! -}
165 (0::Int) /= unsafePerformIO (_ccall_ isFloatDenormalized x) -- ..
167 (0::Int) /= unsafePerformIO (_ccall_ isFloatNegativeZero x) -- ...
171 %*********************************************************
173 \subsection{Type @Double@}
175 %*********************************************************
178 instance Show Float where
179 showsPrec x = showSigned showFloat x
180 showList = showList__ (showsPrec 0)
182 instance Eq Double where
183 (D# x) == (D# y) = x ==## y
185 instance Ord Double where
186 (D# x) `compare` (D# y) | x <## y = LT
190 (D# x) < (D# y) = x <## y
191 (D# x) <= (D# y) = x <=## y
192 (D# x) >= (D# y) = x >=## y
193 (D# x) > (D# y) = x >## y
195 instance Num Double where
196 (+) x y = plusDouble x y
197 (-) x y = minusDouble x y
198 negate x = negateDouble x
199 (*) x y = timesDouble x y
201 | otherwise = negateDouble x
202 signum x | x == 0.0 = 0
204 | otherwise = negate 1
205 fromInteger n = encodeFloat n 0
206 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
208 instance Real Double where
209 toRational x = (m%1)*(b%1)^^n
210 where (m,n) = decodeFloat x
213 instance Fractional Double where
214 (/) x y = divideDouble x y
215 fromRational x = fromRat x
218 instance Floating Double where
219 pi = 3.141592653589793238
222 sqrt x = sqrtDouble x
226 asin x = asinDouble x
227 acos x = acosDouble x
228 atan x = atanDouble x
229 sinh x = sinhDouble x
230 cosh x = coshDouble x
231 tanh x = tanhDouble x
232 (**) x y = powerDouble x y
233 logBase x y = log y / log x
235 asinh x = log (x + sqrt (1.0+x*x))
236 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
237 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
239 instance RealFrac Double where
241 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
242 {-# SPECIALIZE truncate :: Double -> Int #-}
243 {-# SPECIALIZE round :: Double -> Int #-}
244 {-# SPECIALIZE ceiling :: Double -> Int #-}
245 {-# SPECIALIZE floor :: Double -> Int #-}
247 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
248 {-# SPECIALIZE truncate :: Double -> Integer #-}
249 {-# SPECIALIZE round :: Double -> Integer #-}
250 {-# SPECIALIZE ceiling :: Double -> Integer #-}
251 {-# SPECIALIZE floor :: Double -> Integer #-}
253 #if defined(__UNBOXED_INSTANCES__)
254 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
255 {-# SPECIALIZE truncate :: Double -> Int# #-}
256 {-# SPECIALIZE round :: Double -> Int# #-}
257 {-# SPECIALIZE ceiling :: Double -> Int# #-}
258 {-# SPECIALIZE floor :: Double -> Int# #-}
262 = case (decodeFloat x) of { (m,n) ->
263 let b = floatRadix x in
265 (fromInteger m * fromInteger b ^ n, 0.0)
267 case (quotRem m (b^(negate n))) of { (w,r) ->
268 (fromInteger w, encodeFloat r n)
272 truncate x = case properFraction x of
275 round x = case properFraction x of
277 m = if r < 0.0 then n - 1 else n + 1
278 half_down = abs r - 0.5
280 case (compare half_down 0.0) of
282 EQ -> if even n then n else m
285 ceiling x = case properFraction x of
286 (n,r) -> if r > 0.0 then n + 1 else n
288 floor x = case properFraction x of
289 (n,r) -> if r < 0.0 then n - 1 else n
291 instance RealFloat Double where
292 floatRadix _ = FLT_RADIX -- from float.h
293 floatDigits _ = DBL_MANT_DIG -- ditto
294 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
297 = case decodeDouble# x# of
298 (# exp#, a#, s#, d# #) -> (J# a# s# d#, I# exp#)
300 encodeFloat (J# a# s# d#) (I# e#)
301 = case encodeDouble# a# s# d# e# of { dbl# -> D# dbl# }
303 exponent x = case decodeFloat x of
304 (m,n) -> if m == 0 then 0 else n + floatDigits x
306 significand x = case decodeFloat x of
307 (m,_) -> encodeFloat m (negate (floatDigits x))
309 scaleFloat k x = case decodeFloat x of
310 (m,n) -> encodeFloat m (n+k)
312 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
314 (0::Int) /= unsafePerformIO (_ccall_ isDoubleInfinite x) {- ditto -}
316 (0::Int) /= unsafePerformIO (_ccall_ isDoubleDenormalized x) -- ..
318 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNegativeZero x) -- ...
321 instance Show Double where
322 showsPrec x = showSigned showFloat x
323 showList = showList__ (showsPrec 0)
326 %*********************************************************
328 \subsection{Coercions}
330 %*********************************************************
333 {- SPECIALIZE fromIntegral ::
343 Integer -> Double #-}
344 fromIntegral :: (Integral a, Num b) => a -> b
345 fromIntegral = fromInteger . toInteger
347 {- SPECIALIZE realToFrac ::
352 Rational -> Rational,
357 realToFrac :: (Real a, Fractional b) => a -> b
358 realToFrac = fromRational . toRational
361 %*********************************************************
363 \subsection{Common code for @Float@ and @Double@}
365 %*********************************************************
367 The @Enum@ instances for Floats and Doubles are slightly unusual.
368 The @toEnum@ function truncates numbers to Int. The definitions
369 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
370 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
371 dubious. This example may have either 10 or 11 elements, depending on
372 how 0.1 is represented.
374 NOTE: The instances for Float and Double do not make use of the default
375 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
376 a `non-lossy' conversion to and from Ints. Instead we make use of the
377 1.2 default methods (back in the days when Enum had Ord as a superclass)
378 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
381 instance Enum Float where
384 toEnum = fromIntegral
385 fromEnum = fromInteger . truncate -- may overflow
386 enumFrom = numericEnumFrom
387 enumFromTo = numericEnumFromTo
388 enumFromThen = numericEnumFromThen
389 enumFromThenTo = numericEnumFromThenTo
391 instance Enum Double where
394 toEnum = fromIntegral
395 fromEnum = fromInteger . truncate -- may overflow
396 enumFrom = numericEnumFrom
397 enumFromTo = numericEnumFromTo
398 enumFromThen = numericEnumFromThen
399 enumFromThenTo = numericEnumFromThenTo
401 numericEnumFrom :: (Fractional a) => a -> [a]
402 numericEnumFrom = iterate (+1)
404 numericEnumFromThen :: (Fractional a) => a -> a -> [a]
405 numericEnumFromThen n m = iterate (+(m-n)) n
407 numericEnumFromTo :: (Ord a, Fractional a) => a -> a -> [a]
408 numericEnumFromTo n m = takeWhile (<= m + 1/2) (numericEnumFrom n)
410 numericEnumFromThenTo :: (Ord a, Fractional a) => a -> a -> a -> [a]
411 numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
414 pred | e2 > e1 = (<= e3 + mid)
415 | otherwise = (>= e3 + mid)
419 @approxRational@, applied to two real fractional numbers x and epsilon,
420 returns the simplest rational number within epsilon of x. A rational
421 number n%d in reduced form is said to be simpler than another n'%d' if
422 abs n <= abs n' && d <= d'. Any real interval contains a unique
423 simplest rational; here, for simplicity, we assume a closed rational
424 interval. If such an interval includes at least one whole number, then
425 the simplest rational is the absolutely least whole number. Otherwise,
426 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
427 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
428 the simplest rational between d'%r' and d%r.
431 approxRational :: (RealFrac a) => a -> a -> Rational
432 approxRational rat eps = simplest (rat-eps) (rat+eps)
433 where simplest x y | y < x = simplest y x
435 | x > 0 = simplest' n d n' d'
436 | y < 0 = - simplest' (-n') d' (-n) d
438 where xr = toRational x
445 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
447 | q /= q' = (q+1) :% 1
448 | otherwise = (q*n''+d'') :% n''
449 where (q,r) = quotRem n d
450 (q',r') = quotRem n' d'
451 nd'' = simplest' d' r' d r
453 d'' = denominator nd''
458 instance (Integral a) => Ord (Ratio a) where
459 (x:%y) <= (x':%y') = x * y' <= x' * y
460 (x:%y) < (x':%y') = x * y' < x' * y
462 instance (Integral a) => Num (Ratio a) where
463 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
464 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
465 (x:%y) * (x':%y') = reduce (x * x') (y * y')
466 negate (x:%y) = (-x) :% y
467 abs (x:%y) = abs x :% y
468 signum (x:%_) = signum x :% 1
469 fromInteger x = fromInteger x :% 1
471 instance (Integral a) => Real (Ratio a) where
472 toRational (x:%y) = toInteger x :% toInteger y
474 instance (Integral a) => Fractional (Ratio a) where
475 (x:%y) / (x':%y') = (x*y') % (y*x')
476 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
477 fromRational (x:%y) = fromInteger x :% fromInteger y
479 instance (Integral a) => RealFrac (Ratio a) where
480 properFraction (x:%y) = (fromIntegral q, r:%y)
481 where (q,r) = quotRem x y
483 instance (Integral a) => Enum (Ratio a) where
487 toEnum n = fromIntegral n :% 1
488 fromEnum = fromInteger . truncate
490 enumFrom = bounded_iterator True (1)
491 enumFromThen n m = bounded_iterator (diff >= 0) diff n
495 bounded_iterator :: (Ord a, Num a) => Bool -> a -> a -> [a]
496 bounded_iterator inc step v
497 | inc && v > new_v = [v] -- oflow
498 | not inc && v < new_v = [v] -- uflow
499 | otherwise = v : bounded_iterator inc step new_v
506 instance (Integral a) => Show (Ratio a) where
507 showsPrec p (x:%y) = showParen (p > ratio_prec)
508 (shows x . showString " % " . shows y)
511 @showRational@ converts a Rational to a string that looks like a
512 floating point number, but without converting to any floating type
513 (because of the possible overflow).
515 From/by Lennart, 94/09/26
518 showRational :: Int -> Rational -> String
523 let (r', e) = normalize r
529 -- make sure 1 <= r < 10
530 normalize :: Rational -> (Rational, Int)
531 normalize r = if r < 1 then
532 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
535 where norm :: Int -> Rational -> Int -> (Rational, Int)
536 -- Invariant: x*10^e == original r
541 in if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
543 prR :: Int -> Rational -> Int -> String
544 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
545 prR n r e | r >= 10 = prR n (r/10) (e+1)
547 | e > 0 && e < 8 = takeN e s ('.' : drop0 (drop e s) [])
548 | e <= 0 && e > -3 = '0': '.' : takeN (-e) (repeat '0') (drop0 s [])
549 | otherwise = h : '.' : drop0 t ('e':show e0)
551 s@(h:t) = show ((round (r * 10^n))::Integer)
554 takeN (I# n#) ls rs = takeUInt n# ls rs
556 drop0 :: String -> String -> String
558 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
560 dropTrailing0s [] = Nothing
561 dropTrailing0s ('0':xs) =
562 case dropTrailing0s xs of
564 Just ls -> Just ('0':ls)
565 dropTrailing0s (x:xs) =
566 case dropTrailing0s xs of
568 Just ls -> Just (x:ls)
572 [In response to a request for documentation of how fromRational works,
573 Joe Fasel writes:] A quite reasonable request! This code was added to
574 the Prelude just before the 1.2 release, when Lennart, working with an
575 early version of hbi, noticed that (read . show) was not the identity
576 for floating-point numbers. (There was a one-bit error about half the
577 time.) The original version of the conversion function was in fact
578 simply a floating-point divide, as you suggest above. The new version
579 is, I grant you, somewhat denser.
581 Unfortunately, Joe's code doesn't work! Here's an example:
583 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
588 1.8217369128763981e-300
590 Lennart's code follows, and it works...
593 {-# SPECIALISE fromRat ::
595 Rational -> Float #-}
596 fromRat :: (RealFloat a) => Rational -> a
600 -- If the exponent of the nearest floating-point number to x
601 -- is e, then the significand is the integer nearest xb^(-e),
602 -- where b is the floating-point radix. We start with a good
603 -- guess for e, and if it is correct, the exponent of the
604 -- floating-point number we construct will again be e. If
605 -- not, one more iteration is needed.
607 f e = if e' == e then y else f e'
608 where y = encodeFloat (round (x * (1 % b)^^e)) e
609 (_,e') = decodeFloat y
612 -- We obtain a trial exponent by doing a floating-point
613 -- division of x's numerator by its denominator. The
614 -- result of this division may not itself be the ultimate
615 -- result, because of an accumulation of three rounding
618 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
619 / fromInteger (denominator x))
622 Now, here's Lennart's code.
625 fromRat :: (RealFloat a) => Rational -> a
627 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
628 | x < 0 = - fromRat' (-x) -- first.
629 | otherwise = fromRat' x
631 -- Conversion process:
632 -- Scale the rational number by the RealFloat base until
633 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
634 -- Then round the rational to an Integer and encode it with the exponent
635 -- that we got from the scaling.
636 -- To speed up the scaling process we compute the log2 of the number to get
637 -- a first guess of the exponent.
639 fromRat' :: (RealFloat a) => Rational -> a
641 where b = floatRadix r
643 (minExp0, _) = floatRange r
644 minExp = minExp0 - p -- the real minimum exponent
645 xMin = toRational (expt b (p-1))
646 xMax = toRational (expt b p)
647 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
648 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
649 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
650 r = encodeFloat (round x') p'
652 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
653 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
654 scaleRat b minExp xMin xMax p x
655 | p <= minExp = (x, p)
656 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
657 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
660 -- Exponentiation with a cache for the most common numbers.
661 minExpt, maxExpt :: Int
665 expt :: Integer -> Int -> Integer
667 if base == 2 && n >= minExpt && n <= maxExpt then
672 expts :: Array Int Integer
673 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
675 -- Compute the (floor of the) log of i in base b.
676 -- Simplest way would be just divide i by b until it's smaller then b, but that would
677 -- be very slow! We are just slightly more clever.
678 integerLogBase :: Integer -> Integer -> Int
681 | otherwise = doDiv (i `div` (b^l)) l
683 -- Try squaring the base first to cut down the number of divisions.
684 l = 2 * integerLogBase (b*b) i
686 doDiv :: Integer -> Int -> Int
689 | otherwise = doDiv (x `div` b) (y+1)
693 %*********************************************************
695 \subsection{Printing out numbers}
697 %*********************************************************
700 --Exported from std library Numeric, defined here to
701 --avoid mut. rec. between PrelNum and Numeric.
702 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
703 showSigned showPos p x
704 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
705 | otherwise = showPos x
707 showFloat :: (RealFloat a) => a -> ShowS
708 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
710 -- These are the format types. This type is not exported.
712 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
714 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
715 formatRealFloat fmt decs x
717 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
718 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
719 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
723 doFmt format (is, e) =
724 let ds = map intToDigit is in
727 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
732 let e' = if e==0 then 0 else e-1 in
735 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
737 let dec' = max dec 1 in
739 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
742 (ei,is') = roundTo base (dec'+1) is
743 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
745 d:'.':ds' ++ 'e':show (e-1+ei)
748 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
753 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
754 f n s "" = f (n-1) ('0':s) ""
755 f n s (r:rs) = f (n-1) (r:s) rs
759 let dec' = max dec 0 in
762 (ei,is') = roundTo base (dec' + e) is
763 (ls,rs) = splitAt (e+ei) (map intToDigit is')
765 mk0 ls ++ (if null rs then "" else '.':rs)
768 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
769 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
774 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
782 f n [] = (0, replicate n 0)
783 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
785 | i' == base = (1,0:ds)
786 | otherwise = (0,i':ds)
792 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
793 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
794 -- This version uses a much slower logarithm estimator. It should be improved.
796 -- This function returns a list of digits (Ints in [0..base-1]) and an
799 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
800 floatToDigits _ 0 = ([0], 0)
801 floatToDigits base x =
803 (f0, e0) = decodeFloat x
804 (minExp0, _) = floatRange x
807 minExp = minExp0 - p -- the real minimum exponent
808 -- Haskell requires that f be adjusted so denormalized numbers
809 -- will have an impossibly low exponent. Adjust for this.
811 let n = minExp - e0 in
812 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
817 (f*be*b*2, 2*b, be*b, b)
821 if e > minExp && f == b^(p-1) then
822 (f*b*2, b^(-e+1)*2, b, 1)
824 (f*2, b^(-e)*2, 1, 1)
828 if b == 2 && base == 10 then
829 -- logBase 10 2 is slightly bigger than 3/10 so
830 -- the following will err on the low side. Ignoring
831 -- the fraction will make it err even more.
832 -- Haskell promises that p-1 <= logBase b f < p.
833 (p - 1 + e0) * 3 `div` 10
835 ceiling ((log (fromInteger (f+1)) +
836 fromInt e * log (fromInteger b)) /
837 log (fromInteger base))
838 --WAS: fromInt e * log (fromInteger b))
842 if r + mUp <= expt base n * s then n else fixup (n+1)
844 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
848 gen ds rn sN mUpN mDnN =
850 (dn, rn') = (rn * base) `divMod` sN
854 case (rn' < mDnN', rn' + mUpN' > sN) of
855 (True, False) -> dn : ds
856 (False, True) -> dn+1 : ds
857 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
858 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
862 gen [] r (s * expt base k) mUp mDn
864 let bk = expt base (-k) in
865 gen [] (r * bk) s (mUp * bk) (mDn * bk)
867 (map toInt (reverse rds), k)
871 %*********************************************************
873 \subsection{Numeric primops}
875 %*********************************************************
877 Definitions of the boxed PrimOps; these will be
878 used in the case of partial applications, etc.
881 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
882 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
883 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
884 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
885 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
887 negateFloat :: Float -> Float
888 negateFloat (F# x) = F# (negateFloat# x)
890 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
891 gtFloat (F# x) (F# y) = gtFloat# x y
892 geFloat (F# x) (F# y) = geFloat# x y
893 eqFloat (F# x) (F# y) = eqFloat# x y
894 neFloat (F# x) (F# y) = neFloat# x y
895 ltFloat (F# x) (F# y) = ltFloat# x y
896 leFloat (F# x) (F# y) = leFloat# x y
898 float2Int :: Float -> Int
899 float2Int (F# x) = I# (float2Int# x)
901 int2Float :: Int -> Float
902 int2Float (I# x) = F# (int2Float# x)
904 expFloat, logFloat, sqrtFloat :: Float -> Float
905 sinFloat, cosFloat, tanFloat :: Float -> Float
906 asinFloat, acosFloat, atanFloat :: Float -> Float
907 sinhFloat, coshFloat, tanhFloat :: Float -> Float
908 expFloat (F# x) = F# (expFloat# x)
909 logFloat (F# x) = F# (logFloat# x)
910 sqrtFloat (F# x) = F# (sqrtFloat# x)
911 sinFloat (F# x) = F# (sinFloat# x)
912 cosFloat (F# x) = F# (cosFloat# x)
913 tanFloat (F# x) = F# (tanFloat# x)
914 asinFloat (F# x) = F# (asinFloat# x)
915 acosFloat (F# x) = F# (acosFloat# x)
916 atanFloat (F# x) = F# (atanFloat# x)
917 sinhFloat (F# x) = F# (sinhFloat# x)
918 coshFloat (F# x) = F# (coshFloat# x)
919 tanhFloat (F# x) = F# (tanhFloat# x)
921 powerFloat :: Float -> Float -> Float
922 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
924 -- definitions of the boxed PrimOps; these will be
925 -- used in the case of partial applications, etc.
927 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
928 plusDouble (D# x) (D# y) = D# (x +## y)
929 minusDouble (D# x) (D# y) = D# (x -## y)
930 timesDouble (D# x) (D# y) = D# (x *## y)
931 divideDouble (D# x) (D# y) = D# (x /## y)
933 negateDouble :: Double -> Double
934 negateDouble (D# x) = D# (negateDouble# x)
936 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
937 gtDouble (D# x) (D# y) = x >## y
938 geDouble (D# x) (D# y) = x >=## y
939 eqDouble (D# x) (D# y) = x ==## y
940 neDouble (D# x) (D# y) = x /=## y
941 ltDouble (D# x) (D# y) = x <## y
942 leDouble (D# x) (D# y) = x <=## y
944 double2Int :: Double -> Int
945 double2Int (D# x) = I# (double2Int# x)
947 int2Double :: Int -> Double
948 int2Double (I# x) = D# (int2Double# x)
950 double2Float :: Double -> Float
951 double2Float (D# x) = F# (double2Float# x)
952 float2Double :: Float -> Double
953 float2Double (F# x) = D# (float2Double# x)
955 expDouble, logDouble, sqrtDouble :: Double -> Double
956 sinDouble, cosDouble, tanDouble :: Double -> Double
957 asinDouble, acosDouble, atanDouble :: Double -> Double
958 sinhDouble, coshDouble, tanhDouble :: Double -> Double
959 expDouble (D# x) = D# (expDouble# x)
960 logDouble (D# x) = D# (logDouble# x)
961 sqrtDouble (D# x) = D# (sqrtDouble# x)
962 sinDouble (D# x) = D# (sinDouble# x)
963 cosDouble (D# x) = D# (cosDouble# x)
964 tanDouble (D# x) = D# (tanDouble# x)
965 asinDouble (D# x) = D# (asinDouble# x)
966 acosDouble (D# x) = D# (acosDouble# x)
967 atanDouble (D# x) = D# (atanDouble# x)
968 sinhDouble (D# x) = D# (sinhDouble# x)
969 coshDouble (D# x) = D# (coshDouble# x)
970 tanhDouble (D# x) = D# (tanhDouble# x)
972 powerDouble :: Double -> Double -> Double
973 powerDouble (D# x) (D# y) = D# (x **## y)