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 enumFromThen = numericEnumFromThen
398 enumFromThenTo = numericEnumFromThenTo
400 numericEnumFrom :: (Fractional a) => a -> [a]
401 numericEnumFrom = iterate (+1)
403 numericEnumFromThen :: (Fractional a) => a -> a -> [a]
404 numericEnumFromThen n m = iterate (+(m-n)) n
406 numericEnumFromTo :: (Ord a, Fractional a) => a -> a -> [a]
407 numericEnumFromTo n m = takeWhile (<= m + 1/2) (numericEnumFrom n)
409 numericEnumFromThenTo :: (Ord a, Fractional a) => a -> a -> a -> [a]
410 numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
413 pred | e2 > e1 = (<= e3 + mid)
414 | otherwise = (>= e3 + mid)
418 @approxRational@, applied to two real fractional numbers x and epsilon,
419 returns the simplest rational number within epsilon of x. A rational
420 number n%d in reduced form is said to be simpler than another n'%d' if
421 abs n <= abs n' && d <= d'. Any real interval contains a unique
422 simplest rational; here, for simplicity, we assume a closed rational
423 interval. If such an interval includes at least one whole number, then
424 the simplest rational is the absolutely least whole number. Otherwise,
425 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
426 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
427 the simplest rational between d'%r' and d%r.
430 approxRational :: (RealFrac a) => a -> a -> Rational
431 approxRational rat eps = simplest (rat-eps) (rat+eps)
432 where simplest x y | y < x = simplest y x
434 | x > 0 = simplest' n d n' d'
435 | y < 0 = - simplest' (-n') d' (-n) d
437 where xr = toRational x
444 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
446 | q /= q' = (q+1) :% 1
447 | otherwise = (q*n''+d'') :% n''
448 where (q,r) = quotRem n d
449 (q',r') = quotRem n' d'
450 nd'' = simplest' d' r' d r
452 d'' = denominator nd''
457 instance (Integral a) => Ord (Ratio a) where
458 (x:%y) <= (x':%y') = x * y' <= x' * y
459 (x:%y) < (x':%y') = x * y' < x' * y
461 instance (Integral a) => Num (Ratio a) where
462 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
463 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
464 (x:%y) * (x':%y') = reduce (x * x') (y * y')
465 negate (x:%y) = (-x) :% y
466 abs (x:%y) = abs x :% y
467 signum (x:%_) = signum x :% 1
468 fromInteger x = fromInteger x :% 1
470 instance (Integral a) => Real (Ratio a) where
471 toRational (x:%y) = toInteger x :% toInteger y
473 instance (Integral a) => Fractional (Ratio a) where
474 (x:%y) / (x':%y') = (x*y') % (y*x')
475 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
476 fromRational (x:%y) = fromInteger x :% fromInteger y
478 instance (Integral a) => RealFrac (Ratio a) where
479 properFraction (x:%y) = (fromIntegral q, r:%y)
480 where (q,r) = quotRem x y
482 instance (Integral a) => Enum (Ratio a) where
486 toEnum n = fromIntegral n :% 1
487 fromEnum = fromInteger . truncate
489 enumFrom = bounded_iterator True (1)
490 enumFromThen n m = bounded_iterator (diff >= 0) diff n
494 bounded_iterator :: (Ord a, Num a) => Bool -> a -> a -> [a]
495 bounded_iterator inc step v
496 | inc && v > new_v = [v] -- oflow
497 | not inc && v < new_v = [v] -- uflow
498 | otherwise = v : bounded_iterator inc step new_v
505 instance (Integral a) => Show (Ratio a) where
506 showsPrec p (x:%y) = showParen (p > ratio_prec)
507 (shows x . showString " % " . shows y)
510 @showRational@ converts a Rational to a string that looks like a
511 floating point number, but without converting to any floating type
512 (because of the possible overflow).
514 From/by Lennart, 94/09/26
517 showRational :: Int -> Rational -> String
522 let (r', e) = normalize r
528 -- make sure 1 <= r < 10
529 normalize :: Rational -> (Rational, Int)
530 normalize r = if r < 1 then
531 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
534 where norm :: Int -> Rational -> Int -> (Rational, Int)
535 -- Invariant: x*10^e == original r
540 in if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
542 prR :: Int -> Rational -> Int -> String
543 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
544 prR n r e | r >= 10 = prR n (r/10) (e+1)
546 | e > 0 && e < 8 = takeN e s ('.' : drop0 (drop e s) [])
547 | e <= 0 && e > -3 = '0': '.' : takeN (-e) (repeat '0') (drop0 s [])
548 | otherwise = h : '.' : drop0 t ('e':show e0)
550 s@(h:t) = show ((round (r * 10^n))::Integer)
553 takeN (I# n#) ls rs = takeUInt n# ls rs
555 drop0 :: String -> String -> String
557 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
559 dropTrailing0s [] = Nothing
560 dropTrailing0s ('0':xs) =
561 case dropTrailing0s xs of
563 Just ls -> Just ('0':ls)
564 dropTrailing0s (x:xs) =
565 case dropTrailing0s xs of
567 Just ls -> Just (x:ls)
571 [In response to a request for documentation of how fromRational works,
572 Joe Fasel writes:] A quite reasonable request! This code was added to
573 the Prelude just before the 1.2 release, when Lennart, working with an
574 early version of hbi, noticed that (read . show) was not the identity
575 for floating-point numbers. (There was a one-bit error about half the
576 time.) The original version of the conversion function was in fact
577 simply a floating-point divide, as you suggest above. The new version
578 is, I grant you, somewhat denser.
580 Unfortunately, Joe's code doesn't work! Here's an example:
582 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
587 1.8217369128763981e-300
589 Lennart's code follows, and it works...
592 {-# SPECIALISE fromRat ::
594 Rational -> Float #-}
595 fromRat :: (RealFloat a) => Rational -> a
599 -- If the exponent of the nearest floating-point number to x
600 -- is e, then the significand is the integer nearest xb^(-e),
601 -- where b is the floating-point radix. We start with a good
602 -- guess for e, and if it is correct, the exponent of the
603 -- floating-point number we construct will again be e. If
604 -- not, one more iteration is needed.
606 f e = if e' == e then y else f e'
607 where y = encodeFloat (round (x * (1 % b)^^e)) e
608 (_,e') = decodeFloat y
611 -- We obtain a trial exponent by doing a floating-point
612 -- division of x's numerator by its denominator. The
613 -- result of this division may not itself be the ultimate
614 -- result, because of an accumulation of three rounding
617 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
618 / fromInteger (denominator x))
621 Now, here's Lennart's code.
624 fromRat :: (RealFloat a) => Rational -> a
626 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
627 | x < 0 = - fromRat' (-x) -- first.
628 | otherwise = fromRat' x
630 -- Conversion process:
631 -- Scale the rational number by the RealFloat base until
632 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
633 -- Then round the rational to an Integer and encode it with the exponent
634 -- that we got from the scaling.
635 -- To speed up the scaling process we compute the log2 of the number to get
636 -- a first guess of the exponent.
638 fromRat' :: (RealFloat a) => Rational -> a
640 where b = floatRadix r
642 (minExp0, _) = floatRange r
643 minExp = minExp0 - p -- the real minimum exponent
644 xMin = toRational (expt b (p-1))
645 xMax = toRational (expt b p)
646 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
647 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
648 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
649 r = encodeFloat (round x') p'
651 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
652 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
653 scaleRat b minExp xMin xMax p x
654 | p <= minExp = (x, p)
655 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
656 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
659 -- Exponentiation with a cache for the most common numbers.
660 minExpt, maxExpt :: Int
664 expt :: Integer -> Int -> Integer
666 if base == 2 && n >= minExpt && n <= maxExpt then
671 expts :: Array Int Integer
672 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
674 -- Compute the (floor of the) log of i in base b.
675 -- Simplest way would be just divide i by b until it's smaller then b, but that would
676 -- be very slow! We are just slightly more clever.
677 integerLogBase :: Integer -> Integer -> Int
680 | otherwise = doDiv (i `div` (b^l)) l
682 -- Try squaring the base first to cut down the number of divisions.
683 l = 2 * integerLogBase (b*b) i
685 doDiv :: Integer -> Int -> Int
688 | otherwise = doDiv (x `div` b) (y+1)
692 %*********************************************************
694 \subsection{Printing out numbers}
696 %*********************************************************
699 --Exported from std library Numeric, defined here to
700 --avoid mut. rec. between PrelNum and Numeric.
701 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
702 showSigned showPos p x
703 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
704 | otherwise = showPos x
706 showFloat :: (RealFloat a) => a -> ShowS
707 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
709 -- These are the format types. This type is not exported.
711 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
713 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
714 formatRealFloat fmt decs x
716 | isInfinite x && x < 0 = if x < 0 then "-Infinity" else "Infinity"
717 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
718 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
722 doFmt format (is, e) =
723 let ds = map intToDigit is in
726 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
731 let e' = if e==0 then 0 else e-1 in
733 [d] -> d : ".0e" ++ show e'
734 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
736 let dec' = max dec 1 in
738 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
741 (ei,is') = roundTo base (dec'+1) is
742 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
744 d:'.':ds' ++ 'e':show (e-1+ei)
747 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
752 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
753 f n s "" = f (n-1) ('0':s) ""
754 f n s (r:rs) = f (n-1) (r:s) rs
758 let dec' = max dec 0 in
761 (ei,is') = roundTo base (dec' + e) is
762 (ls,rs) = splitAt (e+ei) (map intToDigit is')
764 mk0 ls ++ (if null rs then "" else '.':rs)
767 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
768 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
773 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
781 f n [] = (0, replicate n 0)
782 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
784 | i' == base = (1,0:ds)
785 | otherwise = (0,i':ds)
791 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
792 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
793 -- This version uses a much slower logarithm estimator. It should be improved.
795 -- This function returns a list of digits (Ints in [0..base-1]) and an
798 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
799 floatToDigits _ 0 = ([0], 0)
800 floatToDigits base x =
802 (f0, e0) = decodeFloat x
803 (minExp0, _) = floatRange x
806 minExp = minExp0 - p -- the real minimum exponent
807 -- Haskell requires that f be adjusted so denormalized numbers
808 -- will have an impossibly low exponent. Adjust for this.
810 let n = minExp - e0 in
811 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
816 (f*be*b*2, 2*b, be*b, b)
820 if e > minExp && f == b^(p-1) then
821 (f*b*2, b^(-e+1)*2, b, 1)
823 (f*2, b^(-e)*2, 1, 1)
827 if b == 2 && base == 10 then
828 -- logBase 10 2 is slightly bigger than 3/10 so
829 -- the following will err on the low side. Ignoring
830 -- the fraction will make it err even more.
831 -- Haskell promises that p-1 <= logBase b f < p.
832 (p - 1 + e0) * 3 `div` 10
834 ceiling ((log (fromInteger (f+1)) +
835 fromInt e * log (fromInteger b)) /
836 log (fromInteger base))
837 --WAS: fromInt e * log (fromInteger b))
841 if r + mUp <= expt base n * s then n else fixup (n+1)
843 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
847 gen ds rn sN mUpN mDnN =
849 (dn, rn') = (rn * base) `divMod` sN
853 case (rn' < mDnN', rn' + mUpN' > sN) of
854 (True, False) -> dn : ds
855 (False, True) -> dn+1 : ds
856 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
857 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
861 gen [] r (s * expt base k) mUp mDn
863 let bk = expt base (-k) in
864 gen [] (r * bk) s (mUp * bk) (mDn * bk)
866 (map toInt (reverse rds), k)
870 %*********************************************************
872 \subsection{Numeric primops}
874 %*********************************************************
876 Definitions of the boxed PrimOps; these will be
877 used in the case of partial applications, etc.
880 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
881 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
882 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
883 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
884 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
886 negateFloat :: Float -> Float
887 negateFloat (F# x) = F# (negateFloat# x)
889 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
890 gtFloat (F# x) (F# y) = gtFloat# x y
891 geFloat (F# x) (F# y) = geFloat# x y
892 eqFloat (F# x) (F# y) = eqFloat# x y
893 neFloat (F# x) (F# y) = neFloat# x y
894 ltFloat (F# x) (F# y) = ltFloat# x y
895 leFloat (F# x) (F# y) = leFloat# x y
897 float2Int :: Float -> Int
898 float2Int (F# x) = I# (float2Int# x)
900 int2Float :: Int -> Float
901 int2Float (I# x) = F# (int2Float# x)
903 expFloat, logFloat, sqrtFloat :: Float -> Float
904 sinFloat, cosFloat, tanFloat :: Float -> Float
905 asinFloat, acosFloat, atanFloat :: Float -> Float
906 sinhFloat, coshFloat, tanhFloat :: Float -> Float
907 expFloat (F# x) = F# (expFloat# x)
908 logFloat (F# x) = F# (logFloat# x)
909 sqrtFloat (F# x) = F# (sqrtFloat# x)
910 sinFloat (F# x) = F# (sinFloat# x)
911 cosFloat (F# x) = F# (cosFloat# x)
912 tanFloat (F# x) = F# (tanFloat# x)
913 asinFloat (F# x) = F# (asinFloat# x)
914 acosFloat (F# x) = F# (acosFloat# x)
915 atanFloat (F# x) = F# (atanFloat# x)
916 sinhFloat (F# x) = F# (sinhFloat# x)
917 coshFloat (F# x) = F# (coshFloat# x)
918 tanhFloat (F# x) = F# (tanhFloat# x)
920 powerFloat :: Float -> Float -> Float
921 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
923 -- definitions of the boxed PrimOps; these will be
924 -- used in the case of partial applications, etc.
926 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
927 plusDouble (D# x) (D# y) = D# (x +## y)
928 minusDouble (D# x) (D# y) = D# (x -## y)
929 timesDouble (D# x) (D# y) = D# (x *## y)
930 divideDouble (D# x) (D# y) = D# (x /## y)
932 negateDouble :: Double -> Double
933 negateDouble (D# x) = D# (negateDouble# x)
935 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
936 gtDouble (D# x) (D# y) = x >## y
937 geDouble (D# x) (D# y) = x >=## y
938 eqDouble (D# x) (D# y) = x ==## y
939 neDouble (D# x) (D# y) = x /=## y
940 ltDouble (D# x) (D# y) = x <## y
941 leDouble (D# x) (D# y) = x <=## y
943 double2Int :: Double -> Int
944 double2Int (D# x) = I# (double2Int# x)
946 int2Double :: Int -> Double
947 int2Double (I# x) = D# (int2Double# x)
949 double2Float :: Double -> Float
950 double2Float (D# x) = F# (double2Float# x)
951 float2Double :: Float -> Double
952 float2Double (F# x) = D# (float2Double# x)
954 expDouble, logDouble, sqrtDouble :: Double -> Double
955 sinDouble, cosDouble, tanDouble :: Double -> Double
956 asinDouble, acosDouble, atanDouble :: Double -> Double
957 sinhDouble, coshDouble, tanhDouble :: Double -> Double
958 expDouble (D# x) = D# (expDouble# x)
959 logDouble (D# x) = D# (logDouble# x)
960 sqrtDouble (D# x) = D# (sqrtDouble# x)
961 sinDouble (D# x) = D# (sinDouble# x)
962 cosDouble (D# x) = D# (cosDouble# x)
963 tanDouble (D# x) = D# (tanDouble# x)
964 asinDouble (D# x) = D# (asinDouble# x)
965 acosDouble (D# x) = D# (acosDouble# x)
966 atanDouble (D# x) = D# (atanDouble# x)
967 sinhDouble (D# x) = D# (sinhDouble# x)
968 coshDouble (D# x) = D# (coshDouble# x)
969 tanhDouble (D# x) = D# (tanhDouble# x)
971 powerDouble :: Double -> Double -> Double
972 powerDouble (D# x) (D# y) = D# (x **## y)