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 takeN (I# n#) ls rs = takeUInt_append n# ls rs
574 drop0 :: String -> String -> String
576 drop0 (c:cs) rs = c : fromMaybe rs (dropTrailing0s cs) --WAS (yuck): reverse (dropWhile (=='0') (reverse cs))
578 dropTrailing0s [] = Nothing
579 dropTrailing0s ('0':xs) =
580 case dropTrailing0s xs of
582 Just ls -> Just ('0':ls)
583 dropTrailing0s (x:xs) =
584 case dropTrailing0s xs of
586 Just ls -> Just (x:ls)
590 [In response to a request for documentation of how fromRational works,
591 Joe Fasel writes:] A quite reasonable request! This code was added to
592 the Prelude just before the 1.2 release, when Lennart, working with an
593 early version of hbi, noticed that (read . show) was not the identity
594 for floating-point numbers. (There was a one-bit error about half the
595 time.) The original version of the conversion function was in fact
596 simply a floating-point divide, as you suggest above. The new version
597 is, I grant you, somewhat denser.
599 Unfortunately, Joe's code doesn't work! Here's an example:
601 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
606 1.8217369128763981e-300
608 Lennart's code follows, and it works...
611 fromRat :: (RealFloat a) => Rational -> a
615 -- If the exponent of the nearest floating-point number to x
616 -- is e, then the significand is the integer nearest xb^(-e),
617 -- where b is the floating-point radix. We start with a good
618 -- guess for e, and if it is correct, the exponent of the
619 -- floating-point number we construct will again be e. If
620 -- not, one more iteration is needed.
622 f e = if e' == e then y else f e'
623 where y = encodeFloat (round (x * (1 % b)^^e)) e
624 (_,e') = decodeFloat y
627 -- We obtain a trial exponent by doing a floating-point
628 -- division of x's numerator by its denominator. The
629 -- result of this division may not itself be the ultimate
630 -- result, because of an accumulation of three rounding
633 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
634 / fromInteger (denominator x))
637 Now, here's Lennart's code.
640 {-# SPECIALISE fromRat ::
642 Rational -> Float #-}
643 fromRat :: (RealFloat a) => Rational -> a
645 | x == 0 = encodeFloat 0 0 -- Handle exceptional cases
646 | x < 0 = - fromRat' (-x) -- first.
647 | otherwise = fromRat' x
649 -- Conversion process:
650 -- Scale the rational number by the RealFloat base until
651 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
652 -- Then round the rational to an Integer and encode it with the exponent
653 -- that we got from the scaling.
654 -- To speed up the scaling process we compute the log2 of the number to get
655 -- a first guess of the exponent.
657 fromRat' :: (RealFloat a) => Rational -> a
659 where b = floatRadix r
661 (minExp0, _) = floatRange r
662 minExp = minExp0 - p -- the real minimum exponent
663 xMin = toRational (expt b (p-1))
664 xMax = toRational (expt b p)
665 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
666 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
667 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
668 r = encodeFloat (round x') p'
670 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
671 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
672 scaleRat b minExp xMin xMax p x
673 | p <= minExp = (x, p)
674 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
675 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
678 -- Exponentiation with a cache for the most common numbers.
679 minExpt, maxExpt :: Int
683 expt :: Integer -> Int -> Integer
685 if base == 2 && n >= minExpt && n <= maxExpt then
690 expts :: Array Int Integer
691 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
693 -- Compute the (floor of the) log of i in base b.
694 -- Simplest way would be just divide i by b until it's smaller then b, but that would
695 -- be very slow! We are just slightly more clever.
696 integerLogBase :: Integer -> Integer -> Int
699 | otherwise = doDiv (i `div` (b^l)) l
701 -- Try squaring the base first to cut down the number of divisions.
702 l = 2 * integerLogBase (b*b) i
704 doDiv :: Integer -> Int -> Int
707 | otherwise = doDiv (x `div` b) (y+1)
711 %*********************************************************
713 \subsection{Printing out numbers}
715 %*********************************************************
718 --Exported from std library Numeric, defined here to
719 --avoid mut. rec. between PrelNum and Numeric.
720 showSigned :: (Real a) => (a -> ShowS) -> Int -> a -> ShowS
721 showSigned showPos p x
722 | x < 0 = showParen (p > 6) (showChar '-' . showPos (-x))
723 | otherwise = showPos x
725 showFloat :: (RealFloat a) => a -> ShowS
726 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
728 -- These are the format types. This type is not exported.
730 data FFFormat = FFExponent | FFFixed | FFGeneric
732 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
733 formatRealFloat fmt decs x
735 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
736 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
737 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
741 doFmt format (is, e) =
742 let ds = map intToDigit is in
745 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
750 let e' = if e==0 then 0 else e-1 in
753 (d:ds') -> d : '.' : ds' ++ "e") ++ show e'
755 let dec' = max dec 1 in
757 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
760 (ei,is') = roundTo base (dec'+1) is
761 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
763 d:'.':ds' ++ 'e':show (e-1+ei)
766 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
771 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
772 f n s "" = f (n-1) ('0':s) ""
773 f n s (r:rs) = f (n-1) (r:s) rs
777 let dec' = max dec 0 in
780 (ei,is') = roundTo base (dec' + e) is
781 (ls,rs) = splitAt (e+ei) (map intToDigit is')
783 mk0 ls ++ (if null rs then "" else '.':rs)
786 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
787 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
792 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
800 f n [] = (0, replicate n 0)
801 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
803 | i' == base = (1,0:ds)
804 | otherwise = (0,i':ds)
810 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
811 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
812 -- This version uses a much slower logarithm estimator. It should be improved.
814 -- This function returns a list of digits (Ints in [0..base-1]) and an
817 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
818 floatToDigits _ 0 = ([0], 0)
819 floatToDigits base x =
821 (f0, e0) = decodeFloat x
822 (minExp0, _) = floatRange x
825 minExp = minExp0 - p -- the real minimum exponent
826 -- Haskell requires that f be adjusted so denormalized numbers
827 -- will have an impossibly low exponent. Adjust for this.
829 let n = minExp - e0 in
830 if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
835 (f*be*b*2, 2*b, be*b, b)
839 if e > minExp && f == b^(p-1) then
840 (f*b*2, b^(-e+1)*2, b, 1)
842 (f*2, b^(-e)*2, 1, 1)
846 if b == 2 && base == 10 then
847 -- logBase 10 2 is slightly bigger than 3/10 so
848 -- the following will err on the low side. Ignoring
849 -- the fraction will make it err even more.
850 -- Haskell promises that p-1 <= logBase b f < p.
851 (p - 1 + e0) * 3 `div` 10
853 ceiling ((log (fromInteger (f+1)) +
854 fromInt e * log (fromInteger b)) /
855 log (fromInteger base))
856 --WAS: fromInt e * log (fromInteger b))
860 if r + mUp <= expt base n * s then n else fixup (n+1)
862 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
866 gen ds rn sN mUpN mDnN =
868 (dn, rn') = (rn * base) `divMod` sN
872 case (rn' < mDnN', rn' + mUpN' > sN) of
873 (True, False) -> dn : ds
874 (False, True) -> dn+1 : ds
875 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
876 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
880 gen [] r (s * expt base k) mUp mDn
882 let bk = expt base (-k) in
883 gen [] (r * bk) s (mUp * bk) (mDn * bk)
885 (map toInt (reverse rds), k)
889 %*********************************************************
891 \subsection{Numeric primops}
893 %*********************************************************
895 Definitions of the boxed PrimOps; these will be
896 used in the case of partial applications, etc.
899 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
900 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
901 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
902 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
903 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
905 negateFloat :: Float -> Float
906 negateFloat (F# x) = F# (negateFloat# x)
908 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
909 gtFloat (F# x) (F# y) = gtFloat# x y
910 geFloat (F# x) (F# y) = geFloat# x y
911 eqFloat (F# x) (F# y) = eqFloat# x y
912 neFloat (F# x) (F# y) = neFloat# x y
913 ltFloat (F# x) (F# y) = ltFloat# x y
914 leFloat (F# x) (F# y) = leFloat# x y
916 float2Int :: Float -> Int
917 float2Int (F# x) = I# (float2Int# x)
919 int2Float :: Int -> Float
920 int2Float (I# x) = F# (int2Float# x)
922 expFloat, logFloat, sqrtFloat :: Float -> Float
923 sinFloat, cosFloat, tanFloat :: Float -> Float
924 asinFloat, acosFloat, atanFloat :: Float -> Float
925 sinhFloat, coshFloat, tanhFloat :: Float -> Float
926 expFloat (F# x) = F# (expFloat# x)
927 logFloat (F# x) = F# (logFloat# x)
928 sqrtFloat (F# x) = F# (sqrtFloat# x)
929 sinFloat (F# x) = F# (sinFloat# x)
930 cosFloat (F# x) = F# (cosFloat# x)
931 tanFloat (F# x) = F# (tanFloat# x)
932 asinFloat (F# x) = F# (asinFloat# x)
933 acosFloat (F# x) = F# (acosFloat# x)
934 atanFloat (F# x) = F# (atanFloat# x)
935 sinhFloat (F# x) = F# (sinhFloat# x)
936 coshFloat (F# x) = F# (coshFloat# x)
937 tanhFloat (F# x) = F# (tanhFloat# x)
939 powerFloat :: Float -> Float -> Float
940 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
942 -- definitions of the boxed PrimOps; these will be
943 -- used in the case of partial applications, etc.
945 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
946 plusDouble (D# x) (D# y) = D# (x +## y)
947 minusDouble (D# x) (D# y) = D# (x -## y)
948 timesDouble (D# x) (D# y) = D# (x *## y)
949 divideDouble (D# x) (D# y) = D# (x /## y)
951 negateDouble :: Double -> Double
952 negateDouble (D# x) = D# (negateDouble# x)
954 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
955 gtDouble (D# x) (D# y) = x >## y
956 geDouble (D# x) (D# y) = x >=## y
957 eqDouble (D# x) (D# y) = x ==## y
958 neDouble (D# x) (D# y) = x /=## y
959 ltDouble (D# x) (D# y) = x <## y
960 leDouble (D# x) (D# y) = x <=## y
962 double2Int :: Double -> Int
963 double2Int (D# x) = I# (double2Int# x)
965 int2Double :: Int -> Double
966 int2Double (I# x) = D# (int2Double# x)
968 double2Float :: Double -> Float
969 double2Float (D# x) = F# (double2Float# x)
970 float2Double :: Float -> Double
971 float2Double (F# x) = D# (float2Double# x)
973 expDouble, logDouble, sqrtDouble :: Double -> Double
974 sinDouble, cosDouble, tanDouble :: Double -> Double
975 asinDouble, acosDouble, atanDouble :: Double -> Double
976 sinhDouble, coshDouble, tanhDouble :: Double -> Double
977 expDouble (D# x) = D# (expDouble# x)
978 logDouble (D# x) = D# (logDouble# x)
979 sqrtDouble (D# x) = D# (sqrtDouble# x)
980 sinDouble (D# x) = D# (sinDouble# x)
981 cosDouble (D# x) = D# (cosDouble# x)
982 tanDouble (D# x) = D# (tanDouble# x)
983 asinDouble (D# x) = D# (asinDouble# x)
984 acosDouble (D# x) = D# (acosDouble# x)
985 atanDouble (D# x) = D# (atanDouble# x)
986 sinhDouble (D# x) = D# (sinhDouble# x)
987 coshDouble (D# x) = D# (coshDouble# x)
988 tanhDouble (D# x) = D# (tanhDouble# x)
990 powerDouble :: Double -> Double -> Double
991 powerDouble (D# x) (D# y) = D# (x **## y)