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 enumFromThen = numericEnumFromThen
388 enumFromThenTo = numericEnumFromThenTo
390 instance Enum Double where
393 toEnum = fromIntegral
394 fromEnum = fromInteger . truncate -- may overflow
395 enumFrom = numericEnumFrom
396 enumFromThen = numericEnumFromThen
397 enumFromThenTo = numericEnumFromThenTo
399 numericEnumFrom :: (Real a) => a -> [a]
400 numericEnumFromThen :: (Real a) => a -> a -> [a]
401 numericEnumFromThenTo :: (Real a) => a -> a -> a -> [a]
402 numericEnumFrom = iterate (+1)
403 numericEnumFromThen n m = iterate (+(m-n)) n
404 numericEnumFromThenTo n m p = takeWhile (if m >= n then (<= p) else (>= p))
405 (numericEnumFromThen n m)
408 @approxRational@, applied to two real fractional numbers x and epsilon,
409 returns the simplest rational number within epsilon of x. A rational
410 number n%d in reduced form is said to be simpler than another n'%d' if
411 abs n <= abs n' && d <= d'. Any real interval contains a unique
412 simplest rational; here, for simplicity, we assume a closed rational
413 interval. If such an interval includes at least one whole number, then
414 the simplest rational is the absolutely least whole number. Otherwise,
415 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
416 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
417 the simplest rational between d'%r' and d%r.
420 approxRational :: (RealFrac a) => a -> a -> Rational
421 approxRational rat eps = simplest (rat-eps) (rat+eps)
422 where simplest x y | y < x = simplest y x
424 | x > 0 = simplest' n d n' d'
425 | y < 0 = - simplest' (-n') d' (-n) d
427 where xr = toRational x
434 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
436 | q /= q' = (q+1) :% 1
437 | otherwise = (q*n''+d'') :% n''
438 where (q,r) = quotRem n d
439 (q',r') = quotRem n' d'
440 nd'' = simplest' d' r' d r
442 d'' = denominator nd''
447 instance (Integral a) => Ord (Ratio a) where
448 (x:%y) <= (x':%y') = x * y' <= x' * y
449 (x:%y) < (x':%y') = x * y' < x' * y
451 instance (Integral a) => Num (Ratio a) where
452 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
453 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
454 (x:%y) * (x':%y') = reduce (x * x') (y * y')
455 negate (x:%y) = (-x) :% y
456 abs (x:%y) = abs x :% y
457 signum (x:%_) = signum x :% 1
458 fromInteger x = fromInteger x :% 1
460 instance (Integral a) => Real (Ratio a) where
461 toRational (x:%y) = toInteger x :% toInteger y
463 instance (Integral a) => Fractional (Ratio a) where
464 (x:%y) / (x':%y') = (x*y') % (y*x')
465 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
466 fromRational (x:%y) = fromInteger x :% fromInteger y
468 instance (Integral a) => RealFrac (Ratio a) where
469 properFraction (x:%y) = (fromIntegral q, r:%y)
470 where (q,r) = quotRem x y
472 instance (Integral a) => Enum (Ratio a) where
475 enumFrom = iterate ((+)1)
476 enumFromThen n m = iterate ((+)(m-n)) n
477 toEnum n = fromIntegral n :% 1
478 fromEnum = fromInteger . truncate
483 instance (Integral a) => Show (Ratio a) where
484 showsPrec p (x:%y) = showParen (p > ratio_prec)
485 (shows x . showString " % " . shows y)
488 @showRational@ converts a Rational to a string that looks like a
489 floating point number, but without converting to any floating type
490 (because of the possible overflow).
492 From/by Lennart, 94/09/26
495 showRational :: Int -> Rational -> String
500 let (r', e) = normalize r
506 -- make sure 1 <= r < 10
507 normalize :: Rational -> (Rational, Int)
508 normalize r = if r < 1 then
509 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
512 where norm :: Int -> Rational -> Int -> (Rational, Int)
513 -- Invariant: x*10^e == original r
518 in if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
520 drop0 :: String -> String
522 drop0 (c:cs) = c : fromMaybe [] (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
524 dropTrailing0s [] = Nothing
525 dropTrailing0s ('0':xs) =
526 case dropTrailing0s xs of
528 Just ls -> Just ('0':ls)
529 dropTrailing0s (x:xs) =
530 case dropTrailing0s xs of
532 Just ls -> Just (x:ls)
534 prR :: Int -> Rational -> Int -> String
535 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
536 prR n r e | r >= 10 = prR n (r/10) (e+1)
538 let s = show ((round (r * 10^n))::Integer)
540 in if e > 0 && e < 8 then
541 take e s ++ "." ++ drop0 (drop e s)
542 else if e <= 0 && e > -3 then
543 "0." ++ take (-e) (repeat '0') ++ drop0 s
545 head s : "."++ drop0 (tail s) ++ "e" ++ show e0
548 [In response to a request for documentation of how fromRational works,
549 Joe Fasel writes:] A quite reasonable request! This code was added to
550 the Prelude just before the 1.2 release, when Lennart, working with an
551 early version of hbi, noticed that (read . show) was not the identity
552 for floating-point numbers. (There was a one-bit error about half the
553 time.) The original version of the conversion function was in fact
554 simply a floating-point divide, as you suggest above. The new version
555 is, I grant you, somewhat denser.
557 Unfortunately, Joe's code doesn't work! Here's an example:
559 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
564 1.8217369128763981e-300
566 Lennart's code follows, and it works...
569 {-# SPECIALISE fromRat ::
571 Rational -> Float #-}
572 fromRat :: (RealFloat a) => Rational -> a
576 -- If the exponent of the nearest floating-point number to x
577 -- is e, then the significand is the integer nearest xb^(-e),
578 -- where b is the floating-point radix. We start with a good
579 -- guess for e, and if it is correct, the exponent of the
580 -- floating-point number we construct will again be e. If
581 -- not, one more iteration is needed.
583 f e = if e' == e then y else f e'
584 where y = encodeFloat (round (x * (1 % b)^^e)) e
585 (_,e') = decodeFloat y
588 -- We obtain a trial exponent by doing a floating-point
589 -- division of x's numerator by its denominator. The
590 -- result of this division may not itself be the ultimate
591 -- result, because of an accumulation of three rounding
594 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
595 / fromInteger (denominator x))
598 Now, here's Lennart's code.
601 fromRat :: (RealFloat a) => Rational -> a
603 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
604 | x < 0 = - fromRat' (-x) -- first.
605 | otherwise = fromRat' x
607 -- Conversion process:
608 -- Scale the rational number by the RealFloat base until
609 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
610 -- Then round the rational to an Integer and encode it with the exponent
611 -- that we got from the scaling.
612 -- To speed up the scaling process we compute the log2 of the number to get
613 -- a first guess of the exponent.
615 fromRat' :: (RealFloat a) => Rational -> a
617 where b = floatRadix r
619 (minExp0, _) = floatRange r
620 minExp = minExp0 - p -- the real minimum exponent
621 xMin = toRational (expt b (p-1))
622 xMax = toRational (expt b p)
623 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
624 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
625 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
626 r = encodeFloat (round x') p'
628 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
629 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
630 scaleRat b minExp xMin xMax p x
631 | p <= minExp = (x, p)
632 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
633 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
636 -- Exponentiation with a cache for the most common numbers.
637 minExpt, maxExpt :: Int
641 expt :: Integer -> Int -> Integer
643 if base == 2 && n >= minExpt && n <= maxExpt then
648 expts :: Array Int Integer
649 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
651 -- Compute the (floor of the) log of i in base b.
652 -- Simplest way would be just divide i by b until it's smaller then b, but that would
653 -- be very slow! We are just slightly more clever.
654 integerLogBase :: Integer -> Integer -> Int
657 | otherwise = doDiv (i `div` (b^l)) l
659 -- Try squaring the base first to cut down the number of divisions.
660 l = 2 * integerLogBase (b*b) i
662 doDiv :: Integer -> Int -> Int
665 | otherwise = doDiv (x `div` b) (y+1)
669 %*********************************************************
671 \subsection{Printing out numbers}
673 %*********************************************************
676 --Exported from std library Numeric, defined here to
677 --avoid mut. rec. between PrelNum and Numeric.
678 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
679 showSigned showPos p x
680 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
681 | otherwise = showPos x
683 showFloat :: (RealFloat a) => a -> ShowS
684 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
686 -- These are the format types. This type is not exported.
688 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
690 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
691 formatRealFloat fmt decs x
693 | isInfinite x && x < 0 = if x < 0 then "-Infinity" else "Infinity"
694 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
695 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
699 doFmt format (is, e) =
700 let ds = map intToDigit is in
703 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
708 let e' = if e==0 then 0 else e-1 in
710 [d] -> d : ".0e" ++ show e'
711 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
713 let dec' = max dec 1 in
715 [0] -> '0':'.':take dec' (repeat '0') ++ "e0"
718 (ei,is') = roundTo base (dec'+1) is
719 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
721 d:'.':ds' ++ 'e':show (e-1+ei)
724 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
729 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
730 f n s "" = f (n-1) ('0':s) ""
731 f n s (r:rs) = f (n-1) (r:s) rs
735 let dec' = max dec 0 in
738 (ei,is') = roundTo base (dec' + e) is
739 (ls,rs) = splitAt (e+ei) (map intToDigit is')
741 mk0 ls ++ (if null rs then "" else '.':rs)
744 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
745 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
750 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
758 f n [] = (0, replicate n 0)
759 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
761 | i' == base = (1,0:ds)
762 | otherwise = (0,i':ds)
768 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
769 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
770 -- This version uses a much slower logarithm estimator. It should be improved.
772 -- This function returns a list of digits (Ints in [0..base-1]) and an
775 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
776 floatToDigits _ 0 = ([0], 0)
777 floatToDigits base x =
779 (f0, e0) = decodeFloat x
780 (minExp0, _) = floatRange x
783 minExp = minExp0 - p -- the real minimum exponent
784 -- Haskell requires that f be adjusted so denormalized numbers
785 -- will have an impossibly low exponent. Adjust for this.
787 let n = minExp - e0 in
788 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
793 (f*be*b*2, 2*b, be*b, b)
797 if e > minExp && f == b^(p-1) then
798 (f*b*2, b^(-e+1)*2, b, 1)
800 (f*2, b^(-e)*2, 1, 1)
804 if b == 2 && base == 10 then
805 -- logBase 10 2 is slightly bigger than 3/10 so
806 -- the following will err on the low side. Ignoring
807 -- the fraction will make it err even more.
808 -- Haskell promises that p-1 <= logBase b f < p.
809 (p - 1 + e0) * 3 `div` 10
811 ceiling ((log (fromInteger (f+1)) +
812 fromInt e * log (fromInteger b)) /
813 log (fromInteger base))
814 --WAS: fromInt e * log (fromInteger b))
818 if r + mUp <= expt base n * s then n else fixup (n+1)
820 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
824 gen ds rn sN mUpN mDnN =
826 (dn, rn') = (rn * base) `divMod` sN
830 case (rn' < mDnN', rn' + mUpN' > sN) of
831 (True, False) -> dn : ds
832 (False, True) -> dn+1 : ds
833 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
834 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
838 gen [] r (s * expt base k) mUp mDn
840 let bk = expt base (-k) in
841 gen [] (r * bk) s (mUp * bk) (mDn * bk)
843 (map toInt (reverse rds), k)
847 %*********************************************************
849 \subsection{Numeric primops}
851 %*********************************************************
853 Definitions of the boxed PrimOps; these will be
854 used in the case of partial applications, etc.
857 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
858 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
859 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
860 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
861 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
863 negateFloat :: Float -> Float
864 negateFloat (F# x) = F# (negateFloat# x)
866 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
867 gtFloat (F# x) (F# y) = gtFloat# x y
868 geFloat (F# x) (F# y) = geFloat# x y
869 eqFloat (F# x) (F# y) = eqFloat# x y
870 neFloat (F# x) (F# y) = neFloat# x y
871 ltFloat (F# x) (F# y) = ltFloat# x y
872 leFloat (F# x) (F# y) = leFloat# x y
874 float2Int :: Float -> Int
875 float2Int (F# x) = I# (float2Int# x)
877 int2Float :: Int -> Float
878 int2Float (I# x) = F# (int2Float# x)
880 expFloat, logFloat, sqrtFloat :: Float -> Float
881 sinFloat, cosFloat, tanFloat :: Float -> Float
882 asinFloat, acosFloat, atanFloat :: Float -> Float
883 sinhFloat, coshFloat, tanhFloat :: Float -> Float
884 expFloat (F# x) = F# (expFloat# x)
885 logFloat (F# x) = F# (logFloat# x)
886 sqrtFloat (F# x) = F# (sqrtFloat# x)
887 sinFloat (F# x) = F# (sinFloat# x)
888 cosFloat (F# x) = F# (cosFloat# x)
889 tanFloat (F# x) = F# (tanFloat# x)
890 asinFloat (F# x) = F# (asinFloat# x)
891 acosFloat (F# x) = F# (acosFloat# x)
892 atanFloat (F# x) = F# (atanFloat# x)
893 sinhFloat (F# x) = F# (sinhFloat# x)
894 coshFloat (F# x) = F# (coshFloat# x)
895 tanhFloat (F# x) = F# (tanhFloat# x)
897 powerFloat :: Float -> Float -> Float
898 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
900 -- definitions of the boxed PrimOps; these will be
901 -- used in the case of partial applications, etc.
903 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
904 plusDouble (D# x) (D# y) = D# (x +## y)
905 minusDouble (D# x) (D# y) = D# (x -## y)
906 timesDouble (D# x) (D# y) = D# (x *## y)
907 divideDouble (D# x) (D# y) = D# (x /## y)
909 negateDouble :: Double -> Double
910 negateDouble (D# x) = D# (negateDouble# x)
912 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
913 gtDouble (D# x) (D# y) = x >## y
914 geDouble (D# x) (D# y) = x >=## y
915 eqDouble (D# x) (D# y) = x ==## y
916 neDouble (D# x) (D# y) = x /=## y
917 ltDouble (D# x) (D# y) = x <## y
918 leDouble (D# x) (D# y) = x <=## y
920 double2Int :: Double -> Int
921 double2Int (D# x) = I# (double2Int# x)
923 int2Double :: Int -> Double
924 int2Double (I# x) = D# (int2Double# x)
926 double2Float :: Double -> Float
927 double2Float (D# x) = F# (double2Float# x)
928 float2Double :: Float -> Double
929 float2Double (F# x) = D# (float2Double# x)
931 expDouble, logDouble, sqrtDouble :: Double -> Double
932 sinDouble, cosDouble, tanDouble :: Double -> Double
933 asinDouble, acosDouble, atanDouble :: Double -> Double
934 sinhDouble, coshDouble, tanhDouble :: Double -> Double
935 expDouble (D# x) = D# (expDouble# x)
936 logDouble (D# x) = D# (logDouble# x)
937 sqrtDouble (D# x) = D# (sqrtDouble# x)
938 sinDouble (D# x) = D# (sinDouble# x)
939 cosDouble (D# x) = D# (cosDouble# x)
940 tanDouble (D# x) = D# (tanDouble# x)
941 asinDouble (D# x) = D# (asinDouble# x)
942 acosDouble (D# x) = D# (acosDouble# x)
943 atanDouble (D# x) = D# (atanDouble# x)
944 sinhDouble (D# x) = D# (sinhDouble# x)
945 coshDouble (D# x) = D# (coshDouble# x)
946 tanhDouble (D# x) = D# (tanhDouble# x)
948 powerDouble :: Double -> Double -> Double
949 powerDouble (D# x) (D# y) = D# (x **## y)