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
23 import PrelErr ( error )
26 import Maybe ( fromMaybe )
28 import PrelArr ( Array, array, (!) )
29 import PrelIOBase ( unsafePerformIO )
30 import PrelCCall () -- we need the definitions of CCallable and
31 -- CReturnable for the _ccall_s herein.
34 %*********************************************************
36 \subsection{Type @Float@}
38 %*********************************************************
41 instance Eq Float where
42 (F# x) == (F# y) = x `eqFloat#` y
44 instance Ord Float where
45 (F# x) `compare` (F# y) | x `ltFloat#` y = LT
49 (F# x) < (F# y) = x `ltFloat#` y
50 (F# x) <= (F# y) = x `leFloat#` y
51 (F# x) >= (F# y) = x `geFloat#` y
52 (F# x) > (F# y) = x `gtFloat#` y
54 instance Num Float where
55 (+) x y = plusFloat x y
56 (-) x y = minusFloat x y
57 negate x = negateFloat x
58 (*) x y = timesFloat x y
60 | otherwise = negateFloat x
61 signum x | x == 0.0 = 0
63 | otherwise = negate 1
65 {-# INLINE fromInteger #-}
66 fromInteger n = encodeFloat n 0
67 -- It's important that encodeFloat inlines here, and that
68 -- fromInteger in turn inlines,
69 -- so that if fromInteger is applied to an (S# i) the right thing happens
71 {-# INLINE fromInt #-}
72 fromInt i = int2Float i
74 instance Real Float where
75 toRational x = (m%1)*(b%1)^^n
76 where (m,n) = decodeFloat x
79 instance Fractional Float where
80 (/) x y = divideFloat x y
81 fromRational x = fromRat x
84 instance Floating Float where
85 pi = 3.141592653589793238
98 (**) x y = powerFloat x y
99 logBase x y = log y / log x
101 asinh x = log (x + sqrt (1.0+x*x))
102 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
103 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
105 instance RealFrac Float where
107 {-# SPECIALIZE properFraction :: Float -> (Int, Float) #-}
108 {-# SPECIALIZE truncate :: Float -> Int #-}
109 {-# SPECIALIZE round :: Float -> Int #-}
110 {-# SPECIALIZE ceiling :: Float -> Int #-}
111 {-# SPECIALIZE floor :: Float -> Int #-}
113 {-# SPECIALIZE properFraction :: Float -> (Integer, Float) #-}
114 {-# SPECIALIZE truncate :: Float -> Integer #-}
115 {-# SPECIALIZE round :: Float -> Integer #-}
116 {-# SPECIALIZE ceiling :: Float -> Integer #-}
117 {-# SPECIALIZE floor :: Float -> Integer #-}
120 = case (decodeFloat x) of { (m,n) ->
121 let b = floatRadix x in
123 (fromInteger m * fromInteger b ^ n, 0.0)
125 case (quotRem m (b^(negate n))) of { (w,r) ->
126 (fromInteger w, encodeFloat r n)
130 truncate x = case properFraction x of
133 round x = case properFraction x of
135 m = if r < 0.0 then n - 1 else n + 1
136 half_down = abs r - 0.5
138 case (compare half_down 0.0) of
140 EQ -> if even n then n else m
143 ceiling x = case properFraction x of
144 (n,r) -> if r > 0.0 then n + 1 else n
146 floor x = case properFraction x of
147 (n,r) -> if r < 0.0 then n - 1 else n
149 foreign import ccall "__encodeFloat" unsafe
150 encodeFloat# :: Int# -> ByteArray# -> Int -> Float
151 foreign import ccall "__int_encodeFloat" unsafe
152 int_encodeFloat# :: Int# -> Int -> Float
155 foreign import ccall "isFloatNaN" unsafe isFloatNaN :: Float -> Int
156 foreign import ccall "isFloatInfinite" unsafe isFloatInfinite :: Float -> Int
157 foreign import ccall "isFloatDenormalized" unsafe isFloatDenormalized :: Float -> Int
158 foreign import ccall "isFloatNegativeZero" unsafe isFloatNegativeZero :: Float -> Int
160 instance RealFloat Float where
161 floatRadix _ = FLT_RADIX -- from float.h
162 floatDigits _ = FLT_MANT_DIG -- ditto
163 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
166 = case decodeFloat# f# of
167 (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
169 encodeFloat (S# i) j = int_encodeFloat# i j
170 encodeFloat (J# s# d#) e = encodeFloat# s# d# e
172 exponent x = case decodeFloat x of
173 (m,n) -> if m == 0 then 0 else n + floatDigits x
175 significand x = case decodeFloat x of
176 (m,_) -> encodeFloat m (negate (floatDigits x))
178 scaleFloat k x = case decodeFloat x of
179 (m,n) -> encodeFloat m (n+k)
180 isNaN x = 0 /= isFloatNaN x
181 isInfinite x = 0 /= isFloatInfinite x
182 isDenormalized x = 0 /= isFloatDenormalized x
183 isNegativeZero x = 0 /= isFloatNegativeZero x
187 %*********************************************************
189 \subsection{Type @Double@}
191 %*********************************************************
194 instance Show Float where
195 showsPrec x = showSigned showFloat x
196 showList = showList__ (showsPrec 0)
198 instance Eq Double where
199 (D# x) == (D# y) = x ==## y
201 instance Ord Double where
202 (D# x) `compare` (D# y) | x <## y = LT
206 (D# x) < (D# y) = x <## y
207 (D# x) <= (D# y) = x <=## y
208 (D# x) >= (D# y) = x >=## y
209 (D# x) > (D# y) = x >## y
211 instance Num Double where
212 (+) x y = plusDouble x y
213 (-) x y = minusDouble x y
214 negate x = negateDouble x
215 (*) x y = timesDouble x y
217 | otherwise = negateDouble x
218 signum x | x == 0.0 = 0
220 | otherwise = negate 1
222 {-# INLINE fromInteger #-}
223 -- See comments with Num Float
224 fromInteger n = encodeFloat n 0
225 fromInt (I# n#) = case (int2Double# n#) of { d# -> D# d# }
227 instance Real Double where
228 toRational x = (m%1)*(b%1)^^n
229 where (m,n) = decodeFloat x
232 instance Fractional Double where
233 (/) x y = divideDouble x y
234 fromRational x = fromRat x
237 instance Floating Double where
238 pi = 3.141592653589793238
241 sqrt x = sqrtDouble x
245 asin x = asinDouble x
246 acos x = acosDouble x
247 atan x = atanDouble x
248 sinh x = sinhDouble x
249 cosh x = coshDouble x
250 tanh x = tanhDouble x
251 (**) x y = powerDouble x y
252 logBase x y = log y / log x
254 asinh x = log (x + sqrt (1.0+x*x))
255 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
256 atanh x = log ((x+1.0) / sqrt (1.0-x*x))
258 instance RealFrac Double where
260 {-# SPECIALIZE properFraction :: Double -> (Int, Double) #-}
261 {-# SPECIALIZE truncate :: Double -> Int #-}
262 {-# SPECIALIZE round :: Double -> Int #-}
263 {-# SPECIALIZE ceiling :: Double -> Int #-}
264 {-# SPECIALIZE floor :: Double -> Int #-}
266 {-# SPECIALIZE properFraction :: Double -> (Integer, Double) #-}
267 {-# SPECIALIZE truncate :: Double -> Integer #-}
268 {-# SPECIALIZE round :: Double -> Integer #-}
269 {-# SPECIALIZE ceiling :: Double -> Integer #-}
270 {-# SPECIALIZE floor :: Double -> Integer #-}
273 = case (decodeFloat x) of { (m,n) ->
274 let b = floatRadix x in
276 (fromInteger m * fromInteger b ^ n, 0.0)
278 case (quotRem m (b^(negate n))) of { (w,r) ->
279 (fromInteger w, encodeFloat r n)
283 truncate x = case properFraction x of
286 round x = case properFraction x of
288 m = if r < 0.0 then n - 1 else n + 1
289 half_down = abs r - 0.5
291 case (compare half_down 0.0) of
293 EQ -> if even n then n else m
296 ceiling x = case properFraction x of
297 (n,r) -> if r > 0.0 then n + 1 else n
299 floor x = case properFraction x of
300 (n,r) -> if r < 0.0 then n - 1 else n
302 foreign import ccall "__encodeDouble" unsafe
303 encodeDouble# :: Int# -> ByteArray# -> Int -> Double
304 foreign import ccall "__int_encodeDouble" unsafe
305 int_encodeDouble# :: Int# -> Int -> Double
307 foreign import ccall "isDoubleNaN" unsafe isDoubleNaN :: Double -> Int
308 foreign import ccall "isDoubleInfinite" unsafe isDoubleInfinite :: Double -> Int
309 foreign import ccall "isDoubleDenormalized" unsafe isDoubleDenormalized :: Double -> Int
310 foreign import ccall "isDoubleNegativeZero" unsafe isDoubleNegativeZero :: Double -> Int
312 instance RealFloat Double where
313 floatRadix _ = FLT_RADIX -- from float.h
314 floatDigits _ = DBL_MANT_DIG -- ditto
315 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
318 = case decodeDouble# x# of
319 (# exp#, s#, d# #) -> (J# s# d#, I# exp#)
321 encodeFloat (S# i) j = int_encodeDouble# i j
322 encodeFloat (J# s# d#) e = encodeDouble# s# d# e
324 exponent x = case decodeFloat x of
325 (m,n) -> if m == 0 then 0 else n + floatDigits x
327 significand x = case decodeFloat x of
328 (m,_) -> encodeFloat m (negate (floatDigits x))
330 scaleFloat k x = case decodeFloat x of
331 (m,n) -> encodeFloat m (n+k)
333 isNaN x = 0 /= isDoubleNaN x
334 isInfinite x = 0 /= isDoubleInfinite x
335 isDenormalized x = 0 /= isDoubleDenormalized x
336 isNegativeZero x = 0 /= isDoubleNegativeZero x
339 instance Show Double where
340 showsPrec x = showSigned showFloat x
341 showList = showList__ (showsPrec 0)
344 %*********************************************************
346 \subsection{Coercions}
348 %*********************************************************
351 {-# SPECIALIZE fromIntegral ::
361 Integer -> Double #-}
362 fromIntegral :: (Integral a, Num b) => a -> b
363 fromIntegral = fromInteger . toInteger
365 {-# SPECIALIZE realToFrac ::
370 Rational -> Rational,
375 realToFrac :: (Real a, Fractional b) => a -> b
376 realToFrac = fromRational . toRational
379 %*********************************************************
381 \subsection{Common code for @Float@ and @Double@}
383 %*********************************************************
385 The @Enum@ instances for Floats and Doubles are slightly unusual.
386 The @toEnum@ function truncates numbers to Int. The definitions
387 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
388 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
389 dubious. This example may have either 10 or 11 elements, depending on
390 how 0.1 is represented.
392 NOTE: The instances for Float and Double do not make use of the default
393 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
394 a `non-lossy' conversion to and from Ints. Instead we make use of the
395 1.2 default methods (back in the days when Enum had Ord as a superclass)
396 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
399 instance Enum Float where
402 toEnum = fromIntegral
403 fromEnum = fromInteger . truncate -- may overflow
404 enumFrom = numericEnumFrom
405 enumFromTo = numericEnumFromTo
406 enumFromThen = numericEnumFromThen
407 enumFromThenTo = numericEnumFromThenTo
409 instance Enum Double where
412 toEnum = fromIntegral
413 fromEnum = fromInteger . truncate -- may overflow
414 enumFrom = numericEnumFrom
415 enumFromTo = numericEnumFromTo
416 enumFromThen = numericEnumFromThen
417 enumFromThenTo = numericEnumFromThenTo
419 numericEnumFrom :: (Fractional a) => a -> [a]
420 numericEnumFrom = iterate (+1)
422 numericEnumFromThen :: (Fractional a) => a -> a -> [a]
423 numericEnumFromThen n m = iterate (+(m-n)) n
425 numericEnumFromTo :: (Ord a, Fractional a) => a -> a -> [a]
426 numericEnumFromTo n m = takeWhile (<= m + 1/2) (numericEnumFrom n)
428 numericEnumFromThenTo :: (Ord a, Fractional a) => a -> a -> a -> [a]
429 numericEnumFromThenTo e1 e2 e3 = takeWhile pred (numericEnumFromThen e1 e2)
432 pred | e2 > e1 = (<= e3 + mid)
433 | otherwise = (>= e3 + mid)
437 @approxRational@, applied to two real fractional numbers x and epsilon,
438 returns the simplest rational number within epsilon of x. A rational
439 number n%d in reduced form is said to be simpler than another n'%d' if
440 abs n <= abs n' && d <= d'. Any real interval contains a unique
441 simplest rational; here, for simplicity, we assume a closed rational
442 interval. If such an interval includes at least one whole number, then
443 the simplest rational is the absolutely least whole number. Otherwise,
444 the bounds are of the form q%1 + r%d and q%1 + r'%d', where abs r < d
445 and abs r' < d', and the simplest rational is q%1 + the reciprocal of
446 the simplest rational between d'%r' and d%r.
449 approxRational :: (RealFrac a) => a -> a -> Rational
450 approxRational rat eps = simplest (rat-eps) (rat+eps)
451 where simplest x y | y < x = simplest y x
453 | x > 0 = simplest' n d n' d'
454 | y < 0 = - simplest' (-n') d' (-n) d
456 where xr = toRational x
463 simplest' n d n' d' -- assumes 0 < n%d < n'%d'
465 | q /= q' = (q+1) :% 1
466 | otherwise = (q*n''+d'') :% n''
467 where (q,r) = quotRem n d
468 (q',r') = quotRem n' d'
469 nd'' = simplest' d' r' d r
471 d'' = denominator nd''
476 instance (Integral a) => Ord (Ratio a) where
477 (x:%y) <= (x':%y') = x * y' <= x' * y
478 (x:%y) < (x':%y') = x * y' < x' * y
480 instance (Integral a) => Num (Ratio a) where
481 (x:%y) + (x':%y') = reduce (x*y' + x'*y) (y*y')
482 (x:%y) - (x':%y') = reduce (x*y' - x'*y) (y*y')
483 (x:%y) * (x':%y') = reduce (x * x') (y * y')
484 negate (x:%y) = (-x) :% y
485 abs (x:%y) = abs x :% y
486 signum (x:%_) = signum x :% 1
487 fromInteger x = fromInteger x :% 1
489 instance (Integral a) => Real (Ratio a) where
490 toRational (x:%y) = toInteger x :% toInteger y
492 instance (Integral a) => Fractional (Ratio a) where
493 (x:%y) / (x':%y') = (x*y') % (y*x')
494 recip (x:%y) = if x < 0 then (-y) :% (-x) else y :% x
495 fromRational (x:%y) = fromInteger x :% fromInteger y
497 instance (Integral a) => RealFrac (Ratio a) where
498 properFraction (x:%y) = (fromIntegral q, r:%y)
499 where (q,r) = quotRem x y
501 instance (Integral a) => Enum (Ratio a) where
505 toEnum n = fromIntegral n :% 1
506 fromEnum = fromInteger . truncate
508 enumFrom = bounded_iterator True (1)
509 enumFromThen n m = bounded_iterator (diff >= 0) diff n
513 bounded_iterator :: (Ord a, Num a) => Bool -> a -> a -> [a]
514 bounded_iterator inc step v
515 | inc && v > new_v = [v] -- oflow
516 | not inc && v < new_v = [v] -- uflow
517 | otherwise = v : bounded_iterator inc step new_v
524 instance (Integral a) => Show (Ratio a) where
525 showsPrec p (x:%y) = showParen (p > ratio_prec)
526 (shows x . showString " % " . shows y)
529 @showRational@ converts a Rational to a string that looks like a
530 floating point number, but without converting to any floating type
531 (because of the possible overflow).
533 From/by Lennart, 94/09/26
536 showRational :: Int -> Rational -> String
541 let (r', e) = normalize r
547 -- make sure 1 <= r < 10
548 normalize :: Rational -> (Rational, Int)
549 normalize r = if r < 1 then
550 case norm startExpExp (1 / r) 0 of (r', e) -> (10 / r', -e-1)
553 where norm :: Int -> Rational -> Int -> (Rational, Int)
554 -- Invariant: x*10^e == original r
559 in if x >= tn then norm ee (x/tn) (e+n) else norm (ee-1) x e
561 prR :: Int -> Rational -> Int -> String
562 prR n r e | r < 1 = prR n (r*10) (e-1) -- final adjustment
563 prR n r e | r >= 10 = prR n (r/10) (e+1)
565 | e > 0 && e < 8 = takeN e s ('.' : drop0 (drop e s) [])
566 | e <= 0 && e > -3 = '0': '.' : takeN (-e) (repeat '0') (drop0 s [])
567 | otherwise = h : '.' : drop0 t ('e':show e0)
569 s@(h:t) = show ((round (r * 10^n))::Integer)
572 #ifdef USE_REPORT_PRELUDE
573 takeN n ls rs = take n ls ++ rs
575 takeN (I# n#) ls rs = takeUInt_append n# ls rs
578 drop0 :: String -> String -> String
580 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
582 dropTrailing0s [] = Nothing
583 dropTrailing0s ('0':xs) =
584 case dropTrailing0s xs of
586 Just ls -> Just ('0':ls)
587 dropTrailing0s (x:xs) =
588 case dropTrailing0s xs of
590 Just ls -> Just (x:ls)
594 [In response to a request for documentation of how fromRational works,
595 Joe Fasel writes:] A quite reasonable request! This code was added to
596 the Prelude just before the 1.2 release, when Lennart, working with an
597 early version of hbi, noticed that (read . show) was not the identity
598 for floating-point numbers. (There was a one-bit error about half the
599 time.) The original version of the conversion function was in fact
600 simply a floating-point divide, as you suggest above. The new version
601 is, I grant you, somewhat denser.
603 Unfortunately, Joe's code doesn't work! Here's an example:
605 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
610 1.8217369128763981e-300
612 Lennart's code follows, and it works...
615 fromRat :: (RealFloat a) => Rational -> a
619 -- If the exponent of the nearest floating-point number to x
620 -- is e, then the significand is the integer nearest xb^(-e),
621 -- where b is the floating-point radix. We start with a good
622 -- guess for e, and if it is correct, the exponent of the
623 -- floating-point number we construct will again be e. If
624 -- not, one more iteration is needed.
626 f e = if e' == e then y else f e'
627 where y = encodeFloat (round (x * (1 % b)^^e)) e
628 (_,e') = decodeFloat y
631 -- We obtain a trial exponent by doing a floating-point
632 -- division of x's numerator by its denominator. The
633 -- result of this division may not itself be the ultimate
634 -- result, because of an accumulation of three rounding
637 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
638 / fromInteger (denominator x))
641 Now, here's Lennart's code.
644 {-# SPECIALISE fromRat ::
646 Rational -> Float #-}
647 fromRat :: (RealFloat a) => Rational -> a
649 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
650 | x < 0 = - fromRat' (-x) -- first.
651 | otherwise = fromRat' x
653 -- Conversion process:
654 -- Scale the rational number by the RealFloat base until
655 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
656 -- Then round the rational to an Integer and encode it with the exponent
657 -- that we got from the scaling.
658 -- To speed up the scaling process we compute the log2 of the number to get
659 -- a first guess of the exponent.
661 fromRat' :: (RealFloat a) => Rational -> a
663 where b = floatRadix r
665 (minExp0, _) = floatRange r
666 minExp = minExp0 - p -- the real minimum exponent
667 xMin = toRational (expt b (p-1))
668 xMax = toRational (expt b p)
669 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
670 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
671 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
672 r = encodeFloat (round x') p'
674 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
675 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
676 scaleRat b minExp xMin xMax p x
677 | p <= minExp = (x, p)
678 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
679 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
682 -- Exponentiation with a cache for the most common numbers.
683 minExpt, maxExpt :: Int
687 expt :: Integer -> Int -> Integer
689 if base == 2 && n >= minExpt && n <= maxExpt then
694 expts :: Array Int Integer
695 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
697 -- Compute the (floor of the) log of i in base b.
698 -- Simplest way would be just divide i by b until it's smaller then b, but that would
699 -- be very slow! We are just slightly more clever.
700 integerLogBase :: Integer -> Integer -> Int
703 | otherwise = doDiv (i `div` (b^l)) l
705 -- Try squaring the base first to cut down the number of divisions.
706 l = 2 * integerLogBase (b*b) i
708 doDiv :: Integer -> Int -> Int
711 | otherwise = doDiv (x `div` b) (y+1)
715 %*********************************************************
717 \subsection{Printing out numbers}
719 %*********************************************************
722 --Exported from std library Numeric, defined here to
723 --avoid mut. rec. between PrelNum and Numeric.
724 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
725 showSigned showPos p x
726 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
727 | otherwise = showPos x
729 showFloat :: (RealFloat a) => a -> ShowS
730 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
732 -- These are the format types. This type is not exported.
734 data FFFormat = FFExponent | FFFixed | FFGeneric
736 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
737 formatRealFloat fmt decs x
739 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
740 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
741 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
745 doFmt format (is, e) =
746 let ds = map intToDigit is in
749 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
754 let e' = if e==0 then 0 else e-1 in
757 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
759 let dec' = max dec 1 in
761 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
764 (ei,is') = roundTo base (dec'+1) is
765 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
767 d:'.':ds' ++ 'e':show (e-1+ei)
770 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
775 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
776 f n s "" = f (n-1) ('0':s) ""
777 f n s (r:rs) = f (n-1) (r:s) rs
781 let dec' = max dec 0 in
784 (ei,is') = roundTo base (dec' + e) is
785 (ls,rs) = splitAt (e+ei) (map intToDigit is')
787 mk0 ls ++ (if null rs then "" else '.':rs)
790 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
791 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
796 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
804 f n [] = (0, replicate n 0)
805 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
807 | i' == base = (1,0:ds)
808 | otherwise = (0,i':ds)
814 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
815 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
816 -- This version uses a much slower logarithm estimator. It should be improved.
818 -- This function returns a list of digits (Ints in [0..base-1]) and an
821 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
822 floatToDigits _ 0 = ([0], 0)
823 floatToDigits base x =
825 (f0, e0) = decodeFloat x
826 (minExp0, _) = floatRange x
829 minExp = minExp0 - p -- the real minimum exponent
830 -- Haskell requires that f be adjusted so denormalized numbers
831 -- will have an impossibly low exponent. Adjust for this.
833 let n = minExp - e0 in
834 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
839 (f*be*b*2, 2*b, be*b, b)
843 if e > minExp && f == b^(p-1) then
844 (f*b*2, b^(-e+1)*2, b, 1)
846 (f*2, b^(-e)*2, 1, 1)
850 if b == 2 && base == 10 then
851 -- logBase 10 2 is slightly bigger than 3/10 so
852 -- the following will err on the low side. Ignoring
853 -- the fraction will make it err even more.
854 -- Haskell promises that p-1 <= logBase b f < p.
855 (p - 1 + e0) * 3 `div` 10
857 ceiling ((log (fromInteger (f+1)) +
858 fromInt e * log (fromInteger b)) /
859 log (fromInteger base))
860 --WAS: fromInt e * log (fromInteger b))
864 if r + mUp <= expt base n * s then n else fixup (n+1)
866 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
870 gen ds rn sN mUpN mDnN =
872 (dn, rn') = (rn * base) `divMod` sN
876 case (rn' < mDnN', rn' + mUpN' > sN) of
877 (True, False) -> dn : ds
878 (False, True) -> dn+1 : ds
879 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
880 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
884 gen [] r (s * expt base k) mUp mDn
886 let bk = expt base (-k) in
887 gen [] (r * bk) s (mUp * bk) (mDn * bk)
889 (map toInt (reverse rds), k)
893 %*********************************************************
895 \subsection{Numeric primops}
897 %*********************************************************
899 Definitions of the boxed PrimOps; these will be
900 used in the case of partial applications, etc.
903 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
904 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
905 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
906 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
907 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
909 negateFloat :: Float -> Float
910 negateFloat (F# x) = F# (negateFloat# x)
912 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
913 gtFloat (F# x) (F# y) = gtFloat# x y
914 geFloat (F# x) (F# y) = geFloat# x y
915 eqFloat (F# x) (F# y) = eqFloat# x y
916 neFloat (F# x) (F# y) = neFloat# x y
917 ltFloat (F# x) (F# y) = ltFloat# x y
918 leFloat (F# x) (F# y) = leFloat# x y
920 float2Int :: Float -> Int
921 float2Int (F# x) = I# (float2Int# x)
923 int2Float :: Int -> Float
924 int2Float (I# x) = F# (int2Float# x)
926 expFloat, logFloat, sqrtFloat :: Float -> Float
927 sinFloat, cosFloat, tanFloat :: Float -> Float
928 asinFloat, acosFloat, atanFloat :: Float -> Float
929 sinhFloat, coshFloat, tanhFloat :: Float -> Float
930 expFloat (F# x) = F# (expFloat# x)
931 logFloat (F# x) = F# (logFloat# x)
932 sqrtFloat (F# x) = F# (sqrtFloat# x)
933 sinFloat (F# x) = F# (sinFloat# x)
934 cosFloat (F# x) = F# (cosFloat# x)
935 tanFloat (F# x) = F# (tanFloat# x)
936 asinFloat (F# x) = F# (asinFloat# x)
937 acosFloat (F# x) = F# (acosFloat# x)
938 atanFloat (F# x) = F# (atanFloat# x)
939 sinhFloat (F# x) = F# (sinhFloat# x)
940 coshFloat (F# x) = F# (coshFloat# x)
941 tanhFloat (F# x) = F# (tanhFloat# x)
943 powerFloat :: Float -> Float -> Float
944 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
946 -- definitions of the boxed PrimOps; these will be
947 -- used in the case of partial applications, etc.
949 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
950 plusDouble (D# x) (D# y) = D# (x +## y)
951 minusDouble (D# x) (D# y) = D# (x -## y)
952 timesDouble (D# x) (D# y) = D# (x *## y)
953 divideDouble (D# x) (D# y) = D# (x /## y)
955 negateDouble :: Double -> Double
956 negateDouble (D# x) = D# (negateDouble# x)
958 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
959 gtDouble (D# x) (D# y) = x >## y
960 geDouble (D# x) (D# y) = x >=## y
961 eqDouble (D# x) (D# y) = x ==## y
962 neDouble (D# x) (D# y) = x /=## y
963 ltDouble (D# x) (D# y) = x <## y
964 leDouble (D# x) (D# y) = x <=## y
966 double2Int :: Double -> Int
967 double2Int (D# x) = I# (double2Int# x)
969 int2Double :: Int -> Double
970 int2Double (I# x) = D# (int2Double# x)
972 double2Float :: Double -> Float
973 double2Float (D# x) = F# (double2Float# x)
974 float2Double :: Float -> Double
975 float2Double (F# x) = D# (float2Double# x)
977 expDouble, logDouble, sqrtDouble :: Double -> Double
978 sinDouble, cosDouble, tanDouble :: Double -> Double
979 asinDouble, acosDouble, atanDouble :: Double -> Double
980 sinhDouble, coshDouble, tanhDouble :: Double -> Double
981 expDouble (D# x) = D# (expDouble# x)
982 logDouble (D# x) = D# (logDouble# x)
983 sqrtDouble (D# x) = D# (sqrtDouble# x)
984 sinDouble (D# x) = D# (sinDouble# x)
985 cosDouble (D# x) = D# (cosDouble# x)
986 tanDouble (D# x) = D# (tanDouble# x)
987 asinDouble (D# x) = D# (asinDouble# x)
988 acosDouble (D# x) = D# (acosDouble# x)
989 atanDouble (D# x) = D# (atanDouble# x)
990 sinhDouble (D# x) = D# (sinhDouble# x)
991 coshDouble (D# x) = D# (coshDouble# x)
992 tanhDouble (D# x) = D# (tanhDouble# x)
994 powerDouble :: Double -> Double -> Double
995 powerDouble (D# x) (D# y) = D# (x **## y)