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#, s#, d# #) -> (J# s# d#, I# exp#)
149 encodeFloat i@(S# _) j = encodeFloat (toBig i) j
150 encodeFloat (J# s# d#) (I# e#)
151 = case encodeFloat# s# d# e# of { flt# -> F# flt# }
153 exponent x = case decodeFloat x of
154 (m,n) -> if m == 0 then 0 else n + floatDigits x
156 significand x = case decodeFloat x of
157 (m,_) -> encodeFloat m (negate (floatDigits x))
159 scaleFloat k x = case decodeFloat x of
160 (m,n) -> encodeFloat m (n+k)
162 (0::Int) /= unsafePerformIO (_ccall_ isFloatNaN x) {- a _pure_function! -}
164 (0::Int) /= unsafePerformIO (_ccall_ isFloatInfinite x) {- ditto! -}
166 (0::Int) /= unsafePerformIO (_ccall_ isFloatDenormalized x) -- ..
168 (0::Int) /= unsafePerformIO (_ccall_ isFloatNegativeZero x) -- ...
172 %*********************************************************
174 \subsection{Type @Double@}
176 %*********************************************************
179 instance Show Float where
180 showsPrec x = showSigned showFloat x
181 showList = showList__ (showsPrec 0)
183 instance Eq Double where
184 (D# x) == (D# y) = x ==## y
186 instance Ord Double where
187 (D# x) `compare` (D# y) | x <## y = LT
191 (D# x) < (D# y) = x <## y
192 (D# x) <= (D# y) = x <=## y
193 (D# x) >= (D# y) = x >=## y
194 (D# x) > (D# y) = x >## y
196 instance Num Double where
197 (+) x y = plusDouble x y
198 (-) x y = minusDouble x y
199 negate x = negateDouble x
200 (*) x y = timesDouble x y
202 | otherwise = negateDouble x
203 signum x | x == 0.0 = 0
205 | otherwise = negate 1
206 fromInteger n = encodeFloat n 0
207 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
209 instance Real Double where
210 toRational x = (m%1)*(b%1)^^n
211 where (m,n) = decodeFloat x
214 instance Fractional Double where
215 (/) x y = divideDouble x y
216 fromRational x = fromRat x
219 instance Floating Double where
220 pi = 3.141592653589793238
223 sqrt x = sqrtDouble x
227 asin x = asinDouble x
228 acos x = acosDouble x
229 atan x = atanDouble x
230 sinh x = sinhDouble x
231 cosh x = coshDouble x
232 tanh x = tanhDouble x
233 (**) x y = powerDouble x y
234 logBase x y = log y / log x
236 asinh x = log (x + sqrt (1.0+x*x))
237 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
238 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
240 instance RealFrac Double where
242 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
243 {-# SPECIALIZE truncate :: Double -> Int #-}
244 {-# SPECIALIZE round :: Double -> Int #-}
245 {-# SPECIALIZE ceiling :: Double -> Int #-}
246 {-# SPECIALIZE floor :: Double -> Int #-}
248 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
249 {-# SPECIALIZE truncate :: Double -> Integer #-}
250 {-# SPECIALIZE round :: Double -> Integer #-}
251 {-# SPECIALIZE ceiling :: Double -> Integer #-}
252 {-# SPECIALIZE floor :: Double -> Integer #-}
254 #if defined(__UNBOXED_INSTANCES__)
255 {-# SPECIALIZE properFraction :: Double -> (Int#, Double) #-}
256 {-# SPECIALIZE truncate :: Double -> Int# #-}
257 {-# SPECIALIZE round :: Double -> Int# #-}
258 {-# SPECIALIZE ceiling :: Double -> Int# #-}
259 {-# SPECIALIZE floor :: Double -> Int# #-}
263 = case (decodeFloat x) of { (m,n) ->
264 let b = floatRadix x in
266 (fromInteger m * fromInteger b ^ n, 0.0)
268 case (quotRem m (b^(negate n))) of { (w,r) ->
269 (fromInteger w, encodeFloat r n)
273 truncate x = case properFraction x of
276 round x = case properFraction x of
278 m = if r < 0.0 then n - 1 else n + 1
279 half_down = abs r - 0.5
281 case (compare half_down 0.0) of
283 EQ -> if even n then n else m
286 ceiling x = case properFraction x of
287 (n,r) -> if r > 0.0 then n + 1 else n
289 floor x = case properFraction x of
290 (n,r) -> if r < 0.0 then n - 1 else n
292 instance RealFloat Double where
293 floatRadix _ = FLT_RADIX -- from float.h
294 floatDigits _ = DBL_MANT_DIG -- ditto
295 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
298 = case decodeDouble# x# of
299 (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
301 encodeFloat i@(S# _) j = encodeFloat (toBig i) j
302 encodeFloat (J# s# d#) (I# e#)
303 = case encodeDouble# s# d# e# of { dbl# -> D# dbl# }
305 exponent x = case decodeFloat x of
306 (m,n) -> if m == 0 then 0 else n + floatDigits x
308 significand x = case decodeFloat x of
309 (m,_) -> encodeFloat m (negate (floatDigits x))
311 scaleFloat k x = case decodeFloat x of
312 (m,n) -> encodeFloat m (n+k)
314 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNaN x) {- a _pure_function! -}
316 (0::Int) /= unsafePerformIO (_ccall_ isDoubleInfinite x) {- ditto -}
318 (0::Int) /= unsafePerformIO (_ccall_ isDoubleDenormalized x) -- ..
320 (0::Int) /= unsafePerformIO (_ccall_ isDoubleNegativeZero x) -- ...
323 instance Show Double where
324 showsPrec x = showSigned showFloat x
325 showList = showList__ (showsPrec 0)
328 %*********************************************************
330 \subsection{Coercions}
332 %*********************************************************
335 {-# SPECIALIZE fromIntegral ::
345 Integer -> Double #-}
346 fromIntegral :: (Integral a, Num b) => a -> b
347 fromIntegral = fromInteger . toInteger
349 {-# SPECIALIZE realToFrac ::
354 Rational -> Rational,
359 realToFrac :: (Real a, Fractional b) => a -> b
360 realToFrac = fromRational . toRational
363 %*********************************************************
365 \subsection{Common code for @Float@ and @Double@}
367 %*********************************************************
369 The @Enum@ instances for Floats and Doubles are slightly unusual.
370 The @toEnum@ function truncates numbers to Int. The definitions
371 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
372 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
373 dubious. This example may have either 10 or 11 elements, depending on
374 how 0.1 is represented.
376 NOTE: The instances for Float and Double do not make use of the default
377 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
378 a `non-lossy' conversion to and from Ints. Instead we make use of the
379 1.2 default methods (back in the days when Enum had Ord as a superclass)
380 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
383 instance Enum Float where
386 toEnum = fromIntegral
387 fromEnum = fromInteger . truncate -- may overflow
388 enumFrom = numericEnumFrom
389 enumFromTo = numericEnumFromTo
390 enumFromThen = numericEnumFromThen
391 enumFromThenTo = numericEnumFromThenTo
393 instance Enum Double where
396 toEnum = fromIntegral
397 fromEnum = fromInteger . truncate -- may overflow
398 enumFrom = numericEnumFrom
399 enumFromTo = numericEnumFromTo
400 enumFromThen = numericEnumFromThen
401 enumFromThenTo = numericEnumFromThenTo
403 numericEnumFrom :: (Fractional a) => a -> [a]
404 numericEnumFrom = iterate (+1)
406 numericEnumFromThen :: (Fractional a) => a -> a -> [a]
407 numericEnumFromThen n m = iterate (+(m-n)) n
409 numericEnumFromTo :: (Ord a, Fractional a) => a -> a -> [a]
410 numericEnumFromTo n m = takeWhile (<= m + 1/2) (numericEnumFrom n)
412 numericEnumFromThenTo :: (Ord a, Fractional a) => a -> a -> a -> [a]
413 numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
416 pred | e2 > e1 = (<= e3 + mid)
417 | otherwise = (>= e3 + mid)
421 @approxRational@, applied to two real fractional numbers x and epsilon,
422 returns the simplest rational number within epsilon of x. A rational
423 number n%d in reduced form is said to be simpler than another n'%d' if
424 abs n <= abs n' && d <= d'. Any real interval contains a unique
425 simplest rational; here, for simplicity, we assume a closed rational
426 interval. If such an interval includes at least one whole number, then
427 the simplest rational is the absolutely least whole number. Otherwise,
428 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
429 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
430 the simplest rational between d'%r' and d%r.
433 approxRational :: (RealFrac a) => a -> a -> Rational
434 approxRational rat eps = simplest (rat-eps) (rat+eps)
435 where simplest x y | y < x = simplest y x
437 | x > 0 = simplest' n d n' d'
438 | y < 0 = - simplest' (-n') d' (-n) d
440 where xr = toRational x
447 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
449 | q /= q' = (q+1) :% 1
450 | otherwise = (q*n''+d'') :% n''
451 where (q,r) = quotRem n d
452 (q',r') = quotRem n' d'
453 nd'' = simplest' d' r' d r
455 d'' = denominator nd''
460 instance (Integral a) => Ord (Ratio a) where
461 (x:%y) <= (x':%y') = x * y' <= x' * y
462 (x:%y) < (x':%y') = x * y' < x' * y
464 instance (Integral a) => Num (Ratio a) where
465 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
466 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
467 (x:%y) * (x':%y') = reduce (x * x') (y * y')
468 negate (x:%y) = (-x) :% y
469 abs (x:%y) = abs x :% y
470 signum (x:%_) = signum x :% 1
471 fromInteger x = fromInteger x :% 1
473 instance (Integral a) => Real (Ratio a) where
474 toRational (x:%y) = toInteger x :% toInteger y
476 instance (Integral a) => Fractional (Ratio a) where
477 (x:%y) / (x':%y') = (x*y') % (y*x')
478 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
479 fromRational (x:%y) = fromInteger x :% fromInteger y
481 instance (Integral a) => RealFrac (Ratio a) where
482 properFraction (x:%y) = (fromIntegral q, r:%y)
483 where (q,r) = quotRem x y
485 instance (Integral a) => Enum (Ratio a) where
489 toEnum n = fromIntegral n :% 1
490 fromEnum = fromInteger . truncate
492 enumFrom = bounded_iterator True (1)
493 enumFromThen n m = bounded_iterator (diff >= 0) diff n
497 bounded_iterator :: (Ord a, Num a) => Bool -> a -> a -> [a]
498 bounded_iterator inc step v
499 | inc && v > new_v = [v] -- oflow
500 | not inc && v < new_v = [v] -- uflow
501 | otherwise = v : bounded_iterator inc step new_v
508 instance (Integral a) => Show (Ratio a) where
509 showsPrec p (x:%y) = showParen (p > ratio_prec)
510 (shows x . showString " % " . shows y)
513 @showRational@ converts a Rational to a string that looks like a
514 floating point number, but without converting to any floating type
515 (because of the possible overflow).
517 From/by Lennart, 94/09/26
520 showRational :: Int -> Rational -> String
525 let (r', e) = normalize r
531 -- make sure 1 <= r < 10
532 normalize :: Rational -> (Rational, Int)
533 normalize r = if r < 1 then
534 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
537 where norm :: Int -> Rational -> Int -> (Rational, Int)
538 -- Invariant: x*10^e == original r
543 in if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
545 prR :: Int -> Rational -> Int -> String
546 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
547 prR n r e | r >= 10 = prR n (r/10) (e+1)
549 | e > 0 && e < 8 = takeN e s ('.' : drop0 (drop e s) [])
550 | e <= 0 && e > -3 = '0': '.' : takeN (-e) (repeat '0') (drop0 s [])
551 | otherwise = h : '.' : drop0 t ('e':show e0)
553 s@(h:t) = show ((round (r * 10^n))::Integer)
556 takeN (I# n#) ls rs = takeUInt_append n# ls rs
558 drop0 :: String -> String -> String
560 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
562 dropTrailing0s [] = Nothing
563 dropTrailing0s ('0':xs) =
564 case dropTrailing0s xs of
566 Just ls -> Just ('0':ls)
567 dropTrailing0s (x:xs) =
568 case dropTrailing0s xs of
570 Just ls -> Just (x:ls)
574 [In response to a request for documentation of how fromRational works,
575 Joe Fasel writes:] A quite reasonable request! This code was added to
576 the Prelude just before the 1.2 release, when Lennart, working with an
577 early version of hbi, noticed that (read . show) was not the identity
578 for floating-point numbers. (There was a one-bit error about half the
579 time.) The original version of the conversion function was in fact
580 simply a floating-point divide, as you suggest above. The new version
581 is, I grant you, somewhat denser.
583 Unfortunately, Joe's code doesn't work! Here's an example:
585 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
590 1.8217369128763981e-300
592 Lennart's code follows, and it works...
595 {-# SPECIALISE fromRat ::
597 Rational -> Float #-}
598 fromRat :: (RealFloat a) => Rational -> a
602 -- If the exponent of the nearest floating-point number to x
603 -- is e, then the significand is the integer nearest xb^(-e),
604 -- where b is the floating-point radix. We start with a good
605 -- guess for e, and if it is correct, the exponent of the
606 -- floating-point number we construct will again be e. If
607 -- not, one more iteration is needed.
609 f e = if e' == e then y else f e'
610 where y = encodeFloat (round (x * (1 % b)^^e)) e
611 (_,e') = decodeFloat y
614 -- We obtain a trial exponent by doing a floating-point
615 -- division of x's numerator by its denominator. The
616 -- result of this division may not itself be the ultimate
617 -- result, because of an accumulation of three rounding
620 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
621 / fromInteger (denominator x))
624 Now, here's Lennart's code.
627 fromRat :: (RealFloat a) => Rational -> a
629 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
630 | x < 0 = - fromRat' (-x) -- first.
631 | otherwise = fromRat' x
633 -- Conversion process:
634 -- Scale the rational number by the RealFloat base until
635 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
636 -- Then round the rational to an Integer and encode it with the exponent
637 -- that we got from the scaling.
638 -- To speed up the scaling process we compute the log2 of the number to get
639 -- a first guess of the exponent.
641 fromRat' :: (RealFloat a) => Rational -> a
643 where b = floatRadix r
645 (minExp0, _) = floatRange r
646 minExp = minExp0 - p -- the real minimum exponent
647 xMin = toRational (expt b (p-1))
648 xMax = toRational (expt b p)
649 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
650 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
651 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
652 r = encodeFloat (round x') p'
654 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
655 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
656 scaleRat b minExp xMin xMax p x
657 | p <= minExp = (x, p)
658 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
659 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
662 -- Exponentiation with a cache for the most common numbers.
663 minExpt, maxExpt :: Int
667 expt :: Integer -> Int -> Integer
669 if base == 2 && n >= minExpt && n <= maxExpt then
674 expts :: Array Int Integer
675 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
677 -- Compute the (floor of the) log of i in base b.
678 -- Simplest way would be just divide i by b until it's smaller then b, but that would
679 -- be very slow! We are just slightly more clever.
680 integerLogBase :: Integer -> Integer -> Int
683 | otherwise = doDiv (i `div` (b^l)) l
685 -- Try squaring the base first to cut down the number of divisions.
686 l = 2 * integerLogBase (b*b) i
688 doDiv :: Integer -> Int -> Int
691 | otherwise = doDiv (x `div` b) (y+1)
695 %*********************************************************
697 \subsection{Printing out numbers}
699 %*********************************************************
702 --Exported from std library Numeric, defined here to
703 --avoid mut. rec. between PrelNum and Numeric.
704 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
705 showSigned showPos p x
706 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
707 | otherwise = showPos x
709 showFloat :: (RealFloat a) => a -> ShowS
710 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
712 -- These are the format types. This type is not exported.
714 data FFFormat = FFExponent | FFFixed | FFGeneric --no need: deriving (Eq, Ord, Show)
716 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
717 formatRealFloat fmt decs x
719 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
720 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
721 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
725 doFmt format (is, e) =
726 let ds = map intToDigit is in
729 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
734 let e' = if e==0 then 0 else e-1 in
737 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
739 let dec' = max dec 1 in
741 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
744 (ei,is') = roundTo base (dec'+1) is
745 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
747 d:'.':ds' ++ 'e':show (e-1+ei)
750 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
755 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
756 f n s "" = f (n-1) ('0':s) ""
757 f n s (r:rs) = f (n-1) (r:s) rs
761 let dec' = max dec 0 in
764 (ei,is') = roundTo base (dec' + e) is
765 (ls,rs) = splitAt (e+ei) (map intToDigit is')
767 mk0 ls ++ (if null rs then "" else '.':rs)
770 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
771 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
776 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
784 f n [] = (0, replicate n 0)
785 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
787 | i' == base = (1,0:ds)
788 | otherwise = (0,i':ds)
794 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
795 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
796 -- This version uses a much slower logarithm estimator. It should be improved.
798 -- This function returns a list of digits (Ints in [0..base-1]) and an
801 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
802 floatToDigits _ 0 = ([0], 0)
803 floatToDigits base x =
805 (f0, e0) = decodeFloat x
806 (minExp0, _) = floatRange x
809 minExp = minExp0 - p -- the real minimum exponent
810 -- Haskell requires that f be adjusted so denormalized numbers
811 -- will have an impossibly low exponent. Adjust for this.
813 let n = minExp - e0 in
814 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
819 (f*be*b*2, 2*b, be*b, b)
823 if e > minExp && f == b^(p-1) then
824 (f*b*2, b^(-e+1)*2, b, 1)
826 (f*2, b^(-e)*2, 1, 1)
830 if b == 2 && base == 10 then
831 -- logBase 10 2 is slightly bigger than 3/10 so
832 -- the following will err on the low side. Ignoring
833 -- the fraction will make it err even more.
834 -- Haskell promises that p-1 <= logBase b f < p.
835 (p - 1 + e0) * 3 `div` 10
837 ceiling ((log (fromInteger (f+1)) +
838 fromInt e * log (fromInteger b)) /
839 log (fromInteger base))
840 --WAS: fromInt e * log (fromInteger b))
844 if r + mUp <= expt base n * s then n else fixup (n+1)
846 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
850 gen ds rn sN mUpN mDnN =
852 (dn, rn') = (rn * base) `divMod` sN
856 case (rn' < mDnN', rn' + mUpN' > sN) of
857 (True, False) -> dn : ds
858 (False, True) -> dn+1 : ds
859 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
860 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
864 gen [] r (s * expt base k) mUp mDn
866 let bk = expt base (-k) in
867 gen [] (r * bk) s (mUp * bk) (mDn * bk)
869 (map toInt (reverse rds), k)
873 %*********************************************************
875 \subsection{Numeric primops}
877 %*********************************************************
879 Definitions of the boxed PrimOps; these will be
880 used in the case of partial applications, etc.
883 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
884 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
885 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
886 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
887 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
889 negateFloat :: Float -> Float
890 negateFloat (F# x) = F# (negateFloat# x)
892 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
893 gtFloat (F# x) (F# y) = gtFloat# x y
894 geFloat (F# x) (F# y) = geFloat# x y
895 eqFloat (F# x) (F# y) = eqFloat# x y
896 neFloat (F# x) (F# y) = neFloat# x y
897 ltFloat (F# x) (F# y) = ltFloat# x y
898 leFloat (F# x) (F# y) = leFloat# x y
900 float2Int :: Float -> Int
901 float2Int (F# x) = I# (float2Int# x)
903 int2Float :: Int -> Float
904 int2Float (I# x) = F# (int2Float# x)
906 expFloat, logFloat, sqrtFloat :: Float -> Float
907 sinFloat, cosFloat, tanFloat :: Float -> Float
908 asinFloat, acosFloat, atanFloat :: Float -> Float
909 sinhFloat, coshFloat, tanhFloat :: Float -> Float
910 expFloat (F# x) = F# (expFloat# x)
911 logFloat (F# x) = F# (logFloat# x)
912 sqrtFloat (F# x) = F# (sqrtFloat# x)
913 sinFloat (F# x) = F# (sinFloat# x)
914 cosFloat (F# x) = F# (cosFloat# x)
915 tanFloat (F# x) = F# (tanFloat# x)
916 asinFloat (F# x) = F# (asinFloat# x)
917 acosFloat (F# x) = F# (acosFloat# x)
918 atanFloat (F# x) = F# (atanFloat# x)
919 sinhFloat (F# x) = F# (sinhFloat# x)
920 coshFloat (F# x) = F# (coshFloat# x)
921 tanhFloat (F# x) = F# (tanhFloat# x)
923 powerFloat :: Float -> Float -> Float
924 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
926 -- definitions of the boxed PrimOps; these will be
927 -- used in the case of partial applications, etc.
929 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
930 plusDouble (D# x) (D# y) = D# (x +## y)
931 minusDouble (D# x) (D# y) = D# (x -## y)
932 timesDouble (D# x) (D# y) = D# (x *## y)
933 divideDouble (D# x) (D# y) = D# (x /## y)
935 negateDouble :: Double -> Double
936 negateDouble (D# x) = D# (negateDouble# x)
938 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
939 gtDouble (D# x) (D# y) = x >## y
940 geDouble (D# x) (D# y) = x >=## y
941 eqDouble (D# x) (D# y) = x ==## y
942 neDouble (D# x) (D# y) = x /=## y
943 ltDouble (D# x) (D# y) = x <## y
944 leDouble (D# x) (D# y) = x <=## y
946 double2Int :: Double -> Int
947 double2Int (D# x) = I# (double2Int# x)
949 int2Double :: Int -> Double
950 int2Double (I# x) = D# (int2Double# x)
952 double2Float :: Double -> Float
953 double2Float (D# x) = F# (double2Float# x)
954 float2Double :: Float -> Double
955 float2Double (F# x) = D# (float2Double# x)
957 expDouble, logDouble, sqrtDouble :: Double -> Double
958 sinDouble, cosDouble, tanDouble :: Double -> Double
959 asinDouble, acosDouble, atanDouble :: Double -> Double
960 sinhDouble, coshDouble, tanhDouble :: Double -> Double
961 expDouble (D# x) = D# (expDouble# x)
962 logDouble (D# x) = D# (logDouble# x)
963 sqrtDouble (D# x) = D# (sqrtDouble# x)
964 sinDouble (D# x) = D# (sinDouble# x)
965 cosDouble (D# x) = D# (cosDouble# x)
966 tanDouble (D# x) = D# (tanDouble# x)
967 asinDouble (D# x) = D# (asinDouble# x)
968 acosDouble (D# x) = D# (acosDouble# x)
969 atanDouble (D# x) = D# (atanDouble# x)
970 sinhDouble (D# x) = D# (sinhDouble# x)
971 coshDouble (D# x) = D# (coshDouble# x)
972 tanhDouble (D# x) = D# (tanhDouble# x)
974 powerDouble :: Double -> Double -> Double
975 powerDouble (D# x) (D# y) = D# (x **## y)