6 , ForeignFunctionInterface
8 -- We believe we could deorphan this module, by moving lots of things
9 -- around, but we haven't got there yet:
10 {-# OPTIONS_GHC -fno-warn-orphans #-}
11 {-# OPTIONS_HADDOCK hide #-}
13 -----------------------------------------------------------------------------
16 -- Copyright : (c) The University of Glasgow 1994-2002
17 -- License : see libraries/base/LICENSE
19 -- Maintainer : cvs-ghc@haskell.org
20 -- Stability : internal
21 -- Portability : non-portable (GHC Extensions)
23 -- The types 'Float' and 'Double', and the classes 'Floating' and 'RealFloat'.
25 -----------------------------------------------------------------------------
27 #include "ieee-flpt.h"
30 module GHC.Float( module GHC.Float, Float(..), Double(..), Float#, Double#
31 , double2Int, int2Double, float2Int, int2Float )
44 import GHC.Float.RealFracMethods
45 import GHC.Float.ConversionUtils
46 import GHC.Integer.Logarithms ( integerLogBase# )
47 import GHC.Integer.Logarithms.Internals
52 %*********************************************************
54 \subsection{Standard numeric classes}
56 %*********************************************************
59 -- | Trigonometric and hyperbolic functions and related functions.
61 -- Minimal complete definition:
62 -- 'pi', 'exp', 'log', 'sin', 'cos', 'sinh', 'cosh',
63 -- 'asin', 'acos', 'atan', 'asinh', 'acosh' and 'atanh'
64 class (Fractional a) => Floating a where
66 exp, log, sqrt :: a -> a
67 (**), logBase :: a -> a -> a
68 sin, cos, tan :: a -> a
69 asin, acos, atan :: a -> a
70 sinh, cosh, tanh :: a -> a
71 asinh, acosh, atanh :: a -> a
74 {-# INLINE logBase #-}
78 x ** y = exp (log x * y)
79 logBase x y = log y / log x
82 tanh x = sinh x / cosh x
84 -- | Efficient, machine-independent access to the components of a
85 -- floating-point number.
87 -- Minimal complete definition:
88 -- all except 'exponent', 'significand', 'scaleFloat' and 'atan2'
89 class (RealFrac a, Floating a) => RealFloat a where
90 -- | a constant function, returning the radix of the representation
92 floatRadix :: a -> Integer
93 -- | a constant function, returning the number of digits of
94 -- 'floatRadix' in the significand
95 floatDigits :: a -> Int
96 -- | a constant function, returning the lowest and highest values
97 -- the exponent may assume
98 floatRange :: a -> (Int,Int)
99 -- | The function 'decodeFloat' applied to a real floating-point
100 -- number returns the significand expressed as an 'Integer' and an
101 -- appropriately scaled exponent (an 'Int'). If @'decodeFloat' x@
102 -- yields @(m,n)@, then @x@ is equal in value to @m*b^^n@, where @b@
103 -- is the floating-point radix, and furthermore, either @m@ and @n@
104 -- are both zero or else @b^(d-1) <= m < b^d@, where @d@ is the value
105 -- of @'floatDigits' x@. In particular, @'decodeFloat' 0 = (0,0)@.
106 decodeFloat :: a -> (Integer,Int)
107 -- | 'encodeFloat' performs the inverse of 'decodeFloat'
108 encodeFloat :: Integer -> Int -> a
109 -- | the second component of 'decodeFloat'.
111 -- | the first component of 'decodeFloat', scaled to lie in the open
112 -- interval (@-1@,@1@)
113 significand :: a -> a
114 -- | multiplies a floating-point number by an integer power of the radix
115 scaleFloat :: Int -> a -> a
116 -- | 'True' if the argument is an IEEE \"not-a-number\" (NaN) value
118 -- | 'True' if the argument is an IEEE infinity or negative infinity
119 isInfinite :: a -> Bool
120 -- | 'True' if the argument is too small to be represented in
122 isDenormalized :: a -> Bool
123 -- | 'True' if the argument is an IEEE negative zero
124 isNegativeZero :: a -> Bool
125 -- | 'True' if the argument is an IEEE floating point number
127 -- | a version of arctangent taking two real floating-point arguments.
128 -- For real floating @x@ and @y@, @'atan2' y x@ computes the angle
129 -- (from the positive x-axis) of the vector from the origin to the
130 -- point @(x,y)@. @'atan2' y x@ returns a value in the range [@-pi@,
131 -- @pi@]. It follows the Common Lisp semantics for the origin when
132 -- signed zeroes are supported. @'atan2' y 1@, with @y@ in a type
133 -- that is 'RealFloat', should return the same value as @'atan' y@.
134 -- A default definition of 'atan2' is provided, but implementors
135 -- can provide a more accurate implementation.
139 exponent x = if m == 0 then 0 else n + floatDigits x
140 where (m,n) = decodeFloat x
142 significand x = encodeFloat m (negate (floatDigits x))
143 where (m,_) = decodeFloat x
145 scaleFloat k x = encodeFloat m (n + clamp b k)
146 where (m,n) = decodeFloat x
150 -- n+k may overflow, which would lead
151 -- to wrong results, hence we clamp the
152 -- scaling parameter.
153 -- If n + k would be larger than h,
154 -- n + clamp b k must be too, simliar
155 -- for smaller than l - d.
156 -- Add a little extra to keep clear
157 -- from the boundary cases.
161 | x == 0 && y > 0 = pi/2
162 | x < 0 && y > 0 = pi + atan (y/x)
163 |(x <= 0 && y < 0) ||
164 (x < 0 && isNegativeZero y) ||
165 (isNegativeZero x && isNegativeZero y)
167 | y == 0 && (x < 0 || isNegativeZero x)
168 = pi -- must be after the previous test on zero y
169 | x==0 && y==0 = y -- must be after the other double zero tests
170 | otherwise = x + y -- x or y is a NaN, return a NaN (via +)
174 %*********************************************************
176 \subsection{Type @Float@}
178 %*********************************************************
181 instance Num Float where
182 (+) x y = plusFloat x y
183 (-) x y = minusFloat x y
184 negate x = negateFloat x
185 (*) x y = timesFloat x y
187 | otherwise = negateFloat x
188 signum x | x == 0.0 = 0
190 | otherwise = negate 1
192 {-# INLINE fromInteger #-}
193 fromInteger i = F# (floatFromInteger i)
195 instance Real Float where
197 case decodeFloat_Int# x# of
200 (smallInteger m# `shiftLInteger` e#) :% 1
201 | (int2Word# m# `and#` 1##) `eqWord#` 0## ->
202 case elimZerosInt# m# (negateInt# e#) of
203 (# n, d# #) -> n :% shiftLInteger 1 d#
205 smallInteger m# :% shiftLInteger 1 (negateInt# e#)
207 instance Fractional Float where
208 (/) x y = divideFloat x y
214 | n == 0 = encodeFloat 0 0
215 | n < 0 = -(fromRat'' minEx mantDigs (-n) d)
216 | otherwise = fromRat'' minEx mantDigs n d
219 mantDigs = FLT_MANT_DIG
222 -- RULES for Integer and Int
224 "properFraction/Float->Integer" properFraction = properFractionFloatInteger
225 "truncate/Float->Integer" truncate = truncateFloatInteger
226 "floor/Float->Integer" floor = floorFloatInteger
227 "ceiling/Float->Integer" ceiling = ceilingFloatInteger
228 "round/Float->Integer" round = roundFloatInteger
229 "properFraction/Float->Int" properFraction = properFractionFloatInt
230 "truncate/Float->Int" truncate = float2Int
231 "floor/Float->Int" floor = floorFloatInt
232 "ceiling/Float->Int" ceiling = ceilingFloatInt
233 "round/Float->Int" round = roundFloatInt
235 instance RealFrac Float where
237 -- ceiling, floor, and truncate are all small
238 {-# INLINE [1] ceiling #-}
239 {-# INLINE [1] floor #-}
240 {-# INLINE [1] truncate #-}
242 -- We assume that FLT_RADIX is 2 so that we can use more efficient code
244 #error FLT_RADIX must be 2
246 properFraction (F# x#)
247 = case decodeFloat_Int# x# of
253 then (fromIntegral m * (2 ^ n), 0.0)
254 else let i = if m >= 0 then m `shiftR` negate n
255 else negate (negate m `shiftR` negate n)
256 f = m - (i `shiftL` negate n)
257 in (fromIntegral i, encodeFloat (fromIntegral f) n)
259 truncate x = case properFraction x of
262 round x = case properFraction x of
264 m = if r < 0.0 then n - 1 else n + 1
265 half_down = abs r - 0.5
267 case (compare half_down 0.0) of
269 EQ -> if even n then n else m
272 ceiling x = case properFraction x of
273 (n,r) -> if r > 0.0 then n + 1 else n
275 floor x = case properFraction x of
276 (n,r) -> if r < 0.0 then n - 1 else n
278 instance Floating Float where
279 pi = 3.141592653589793238
292 (**) x y = powerFloat x y
293 logBase x y = log y / log x
295 asinh x = log (x + sqrt (1.0+x*x))
296 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
297 atanh x = 0.5 * log ((1.0+x) / (1.0-x))
299 instance RealFloat Float where
300 floatRadix _ = FLT_RADIX -- from float.h
301 floatDigits _ = FLT_MANT_DIG -- ditto
302 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
304 decodeFloat (F# f#) = case decodeFloat_Int# f# of
305 (# i, e #) -> (smallInteger i, I# e)
307 encodeFloat i (I# e) = F# (encodeFloatInteger i e)
309 exponent x = case decodeFloat x of
310 (m,n) -> if m == 0 then 0 else n + floatDigits x
312 significand x = case decodeFloat x of
313 (m,_) -> encodeFloat m (negate (floatDigits x))
315 scaleFloat k x = case decodeFloat x of
316 (m,n) -> encodeFloat m (n + clamp bf k)
317 where bf = FLT_MAX_EXP - (FLT_MIN_EXP) + 4*FLT_MANT_DIG
319 isNaN x = 0 /= isFloatNaN x
320 isInfinite x = 0 /= isFloatInfinite x
321 isDenormalized x = 0 /= isFloatDenormalized x
322 isNegativeZero x = 0 /= isFloatNegativeZero x
325 instance Show Float where
326 showsPrec x = showSignedFloat showFloat x
327 showList = showList__ (showsPrec 0)
330 %*********************************************************
332 \subsection{Type @Double@}
334 %*********************************************************
337 instance Num Double where
338 (+) x y = plusDouble x y
339 (-) x y = minusDouble x y
340 negate x = negateDouble x
341 (*) x y = timesDouble x y
343 | otherwise = negateDouble x
344 signum x | x == 0.0 = 0
346 | otherwise = negate 1
348 {-# INLINE fromInteger #-}
349 fromInteger i = D# (doubleFromInteger i)
352 instance Real Double where
354 case decodeDoubleInteger x# of
357 shiftLInteger m e# :% 1
358 | (int2Word# (toInt# m) `and#` 1##) `eqWord#` 0## ->
359 case elimZerosInteger m (negateInt# e#) of
360 (# n, d# #) -> n :% shiftLInteger 1 d#
362 m :% shiftLInteger 1 (negateInt# e#)
364 instance Fractional Double where
365 (/) x y = divideDouble x y
371 | n == 0 = encodeFloat 0 0
372 | n < 0 = -(fromRat'' minEx mantDigs (-n) d)
373 | otherwise = fromRat'' minEx mantDigs n d
376 mantDigs = DBL_MANT_DIG
379 instance Floating Double where
380 pi = 3.141592653589793238
383 sqrt x = sqrtDouble x
387 asin x = asinDouble x
388 acos x = acosDouble x
389 atan x = atanDouble x
390 sinh x = sinhDouble x
391 cosh x = coshDouble x
392 tanh x = tanhDouble x
393 (**) x y = powerDouble x y
394 logBase x y = log y / log x
396 asinh x = log (x + sqrt (1.0+x*x))
397 acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
398 atanh x = 0.5 * log ((1.0+x) / (1.0-x))
400 -- RULES for Integer and Int
402 "properFraction/Double->Integer" properFraction = properFractionDoubleInteger
403 "truncate/Double->Integer" truncate = truncateDoubleInteger
404 "floor/Double->Integer" floor = floorDoubleInteger
405 "ceiling/Double->Integer" ceiling = ceilingDoubleInteger
406 "round/Double->Integer" round = roundDoubleInteger
407 "properFraction/Double->Int" properFraction = properFractionDoubleInt
408 "truncate/Double->Int" truncate = double2Int
409 "floor/Double->Int" floor = floorDoubleInt
410 "ceiling/Double->Int" ceiling = ceilingDoubleInt
411 "round/Double->Int" round = roundDoubleInt
413 instance RealFrac Double where
415 -- ceiling, floor, and truncate are all small
416 {-# INLINE [1] ceiling #-}
417 {-# INLINE [1] floor #-}
418 {-# INLINE [1] truncate #-}
421 = case (decodeFloat x) of { (m,n) ->
423 (fromInteger m * 2 ^ n, 0.0)
425 case (quotRem m (2^(negate n))) of { (w,r) ->
426 (fromInteger w, encodeFloat r n)
430 truncate x = case properFraction x of
433 round x = case properFraction x of
435 m = if r < 0.0 then n - 1 else n + 1
436 half_down = abs r - 0.5
438 case (compare half_down 0.0) of
440 EQ -> if even n then n else m
443 ceiling x = case properFraction x of
444 (n,r) -> if r > 0.0 then n + 1 else n
446 floor x = case properFraction x of
447 (n,r) -> if r < 0.0 then n - 1 else n
449 instance RealFloat Double where
450 floatRadix _ = FLT_RADIX -- from float.h
451 floatDigits _ = DBL_MANT_DIG -- ditto
452 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
455 = case decodeDoubleInteger x# of
456 (# i, j #) -> (i, I# j)
458 encodeFloat i (I# j) = D# (encodeDoubleInteger i j)
460 exponent x = case decodeFloat x of
461 (m,n) -> if m == 0 then 0 else n + floatDigits x
463 significand x = case decodeFloat x of
464 (m,_) -> encodeFloat m (negate (floatDigits x))
466 scaleFloat k x = case decodeFloat x of
467 (m,n) -> encodeFloat m (n + clamp bd k)
468 where bd = DBL_MAX_EXP - (DBL_MIN_EXP) + 4*DBL_MANT_DIG
470 isNaN x = 0 /= isDoubleNaN x
471 isInfinite x = 0 /= isDoubleInfinite x
472 isDenormalized x = 0 /= isDoubleDenormalized x
473 isNegativeZero x = 0 /= isDoubleNegativeZero x
476 instance Show Double where
477 showsPrec x = showSignedFloat showFloat x
478 showList = showList__ (showsPrec 0)
481 %*********************************************************
483 \subsection{@Enum@ instances}
485 %*********************************************************
487 The @Enum@ instances for Floats and Doubles are slightly unusual.
488 The @toEnum@ function truncates numbers to Int. The definitions
489 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
490 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
491 dubious. This example may have either 10 or 11 elements, depending on
492 how 0.1 is represented.
494 NOTE: The instances for Float and Double do not make use of the default
495 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
496 a `non-lossy' conversion to and from Ints. Instead we make use of the
497 1.2 default methods (back in the days when Enum had Ord as a superclass)
498 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
501 instance Enum Float where
505 fromEnum = fromInteger . truncate -- may overflow
506 enumFrom = numericEnumFrom
507 enumFromTo = numericEnumFromTo
508 enumFromThen = numericEnumFromThen
509 enumFromThenTo = numericEnumFromThenTo
511 instance Enum Double where
515 fromEnum = fromInteger . truncate -- may overflow
516 enumFrom = numericEnumFrom
517 enumFromTo = numericEnumFromTo
518 enumFromThen = numericEnumFromThen
519 enumFromThenTo = numericEnumFromThenTo
523 %*********************************************************
525 \subsection{Printing floating point}
527 %*********************************************************
531 -- | Show a signed 'RealFloat' value to full precision
532 -- using standard decimal notation for arguments whose absolute value lies
533 -- between @0.1@ and @9,999,999@, and scientific notation otherwise.
534 showFloat :: (RealFloat a) => a -> ShowS
535 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
537 -- These are the format types. This type is not exported.
539 data FFFormat = FFExponent | FFFixed | FFGeneric
541 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
542 formatRealFloat fmt decs x
544 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
545 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
546 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
550 doFmt format (is, e) =
551 let ds = map intToDigit is in
554 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
559 let show_e' = show (e-1) in
562 [d] -> d : ".0e" ++ show_e'
563 (d:ds') -> d : '.' : ds' ++ "e" ++ show_e'
564 [] -> error "formatRealFloat/doFmt/FFExponent: []"
566 let dec' = max dec 1 in
568 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
571 (ei,is') = roundTo base (dec'+1) is
572 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
574 d:'.':ds' ++ 'e':show (e-1+ei)
577 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
581 | e <= 0 -> "0." ++ replicate (-e) '0' ++ ds
584 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
585 f n s "" = f (n-1) ('0':s) ""
586 f n s (r:rs) = f (n-1) (r:s) rs
590 let dec' = max dec 0 in
593 (ei,is') = roundTo base (dec' + e) is
594 (ls,rs) = splitAt (e+ei) (map intToDigit is')
596 mk0 ls ++ (if null rs then "" else '.':rs)
599 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
600 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
602 d : (if null ds' then "" else '.':ds')
605 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
610 _ -> error "roundTo: bad Value"
614 f n [] = (0, replicate n 0)
615 f 0 (x:_) = (if x >= b2 then 1 else 0, [])
617 | i' == base = (1,0:ds)
618 | otherwise = (0,i':ds)
623 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
624 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
625 -- This version uses a much slower logarithm estimator. It should be improved.
627 -- | 'floatToDigits' takes a base and a non-negative 'RealFloat' number,
628 -- and returns a list of digits and an exponent.
629 -- In particular, if @x>=0@, and
631 -- > floatToDigits base x = ([d1,d2,...,dn], e)
637 -- (2) @x = 0.d1d2...dn * (base**e)@
639 -- (3) @0 <= di <= base-1@
641 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
642 floatToDigits _ 0 = ([0], 0)
643 floatToDigits base x =
645 (f0, e0) = decodeFloat x
646 (minExp0, _) = floatRange x
649 minExp = minExp0 - p -- the real minimum exponent
650 -- Haskell requires that f be adjusted so denormalized numbers
651 -- will have an impossibly low exponent. Adjust for this.
653 let n = minExp - e0 in
654 if n > 0 then (f0 `quot` (expt b n), e0+n) else (f0, e0)
658 if f == expt b (p-1) then
659 (f*be*b*2, 2*b, be*b, be) -- according to Burger and Dybvig
663 if e > minExp && f == expt b (p-1) then
664 (f*b*2, expt b (-e+1)*2, b, 1)
666 (f*2, expt b (-e)*2, 1, 1)
672 if b == 2 && base == 10 then
673 -- logBase 10 2 is very slightly larger than 8651/28738
674 -- (about 5.3558e-10), so if log x >= 0, the approximation
675 -- k1 is too small, hence we add one and need one fixup step less.
676 -- If log x < 0, the approximation errs rather on the high side.
677 -- That is usually more than compensated for by ignoring the
678 -- fractional part of logBase 2 x, but when x is a power of 1/2
679 -- or slightly larger and the exponent is a multiple of the
680 -- denominator of the rational approximation to logBase 10 2,
681 -- k1 is larger than logBase 10 x. If k1 > 1 + logBase 10 x,
682 -- we get a leading zero-digit we don't want.
683 -- With the approximation 3/10, this happened for
684 -- 0.5^1030, 0.5^1040, ..., 0.5^1070 and values close above.
685 -- The approximation 8651/28738 guarantees k1 < 1 + logBase 10 x
686 -- for IEEE-ish floating point types with exponent fields
687 -- <= 17 bits and mantissae of several thousand bits, earlier
688 -- convergents to logBase 10 2 would fail for long double.
689 -- Using quot instead of div is a little faster and requires
690 -- fewer fixup steps for negative lx.
692 k1 = (lx * 8651) `quot` 28738
693 in if lx >= 0 then k1 + 1 else k1
695 -- f :: Integer, log :: Float -> Float,
696 -- ceiling :: Float -> Int
697 ceiling ((log (fromInteger (f+1) :: Float) +
698 fromIntegral e * log (fromInteger b)) /
699 log (fromInteger base))
700 --WAS: fromInt e * log (fromInteger b))
704 if r + mUp <= expt base n * s then n else fixup (n+1)
706 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
710 gen ds rn sN mUpN mDnN =
712 (dn, rn') = (rn * base) `quotRem` sN
716 case (rn' < mDnN', rn' + mUpN' > sN) of
717 (True, False) -> dn : ds
718 (False, True) -> dn+1 : ds
719 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
720 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
724 gen [] r (s * expt base k) mUp mDn
726 let bk = expt base (-k) in
727 gen [] (r * bk) s (mUp * bk) (mDn * bk)
729 (map fromIntegral (reverse rds), k)
734 %*********************************************************
736 \subsection{Converting from a Rational to a RealFloat
738 %*********************************************************
740 [In response to a request for documentation of how fromRational works,
741 Joe Fasel writes:] A quite reasonable request! This code was added to
742 the Prelude just before the 1.2 release, when Lennart, working with an
743 early version of hbi, noticed that (read . show) was not the identity
744 for floating-point numbers. (There was a one-bit error about half the
745 time.) The original version of the conversion function was in fact
746 simply a floating-point divide, as you suggest above. The new version
747 is, I grant you, somewhat denser.
749 Unfortunately, Joe's code doesn't work! Here's an example:
751 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
756 1.8217369128763981e-300
761 fromRat :: (RealFloat a) => Rational -> a
765 -- If the exponent of the nearest floating-point number to x
766 -- is e, then the significand is the integer nearest xb^(-e),
767 -- where b is the floating-point radix. We start with a good
768 -- guess for e, and if it is correct, the exponent of the
769 -- floating-point number we construct will again be e. If
770 -- not, one more iteration is needed.
772 f e = if e' == e then y else f e'
773 where y = encodeFloat (round (x * (1 % b)^^e)) e
774 (_,e') = decodeFloat y
777 -- We obtain a trial exponent by doing a floating-point
778 -- division of x's numerator by its denominator. The
779 -- result of this division may not itself be the ultimate
780 -- result, because of an accumulation of three rounding
783 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
784 / fromInteger (denominator x))
787 Now, here's Lennart's code (which works)
790 -- | Converts a 'Rational' value into any type in class 'RealFloat'.
792 "fromRat/Float" fromRat = (fromRational :: Rational -> Float)
793 "fromRat/Double" fromRat = (fromRational :: Rational -> Double)
795 fromRat :: (RealFloat a) => Rational -> a
797 -- Deal with special cases first, delegating the real work to fromRat'
798 fromRat (n :% 0) | n > 0 = 1/0 -- +Infinity
799 | n < 0 = -1/0 -- -Infinity
800 | otherwise = 0/0 -- NaN
802 fromRat (n :% d) | n > 0 = fromRat' (n :% d)
803 | n < 0 = - fromRat' ((-n) :% d)
804 | otherwise = encodeFloat 0 0 -- Zero
806 -- Conversion process:
807 -- Scale the rational number by the RealFloat base until
808 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
809 -- Then round the rational to an Integer and encode it with the exponent
810 -- that we got from the scaling.
811 -- To speed up the scaling process we compute the log2 of the number to get
812 -- a first guess of the exponent.
814 fromRat' :: (RealFloat a) => Rational -> a
815 -- Invariant: argument is strictly positive
817 where b = floatRadix r
819 (minExp0, _) = floatRange r
820 minExp = minExp0 - p -- the real minimum exponent
821 xMin = toRational (expt b (p-1))
822 xMax = toRational (expt b p)
823 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
824 f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1
825 (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)
826 r = encodeFloat (round x') p'
828 -- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.
829 scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)
830 scaleRat b minExp xMin xMax p x
831 | p <= minExp = (x, p)
832 | x >= xMax = scaleRat b minExp xMin xMax (p+1) (x/b)
833 | x < xMin = scaleRat b minExp xMin xMax (p-1) (x*b)
836 -- Exponentiation with a cache for the most common numbers.
837 minExpt, maxExpt :: Int
841 expt :: Integer -> Int -> Integer
843 if base == 2 && n >= minExpt && n <= maxExpt then
846 if base == 10 && n <= maxExpt10 then
851 expts :: Array Int Integer
852 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
857 expts10 :: Array Int Integer
858 expts10 = array (minExpt,maxExpt10) [(n,10^n) | n <- [minExpt .. maxExpt10]]
860 -- Compute the (floor of the) log of i in base b.
861 -- Simplest way would be just divide i by b until it's smaller then b, but that would
862 -- be very slow! We are just slightly more clever, except for base 2, where
863 -- we take advantage of the representation of Integers.
864 -- The general case could be improved by a lookup table for
865 -- approximating the result by integerLog2 i / integerLog2 b.
866 integerLogBase :: Integer -> Integer -> Int
869 | b == 2 = I# (integerLog2# i)
870 | otherwise = I# (integerLogBase# b i)
874 Unfortunately, the old conversion code was awfully slow due to
875 a) a slow integer logarithm
876 b) repeated calculation of gcd's
878 For the case of Rational's coming from a Float or Double via toRational,
879 we can exploit the fact that the denominator is a power of two, which for
880 these brings a huge speedup since we need only shift and add instead
883 The below is an adaption of fromRat' for the conversion to
884 Float or Double exploiting the know floatRadix and avoiding
885 divisions as much as possible.
888 {-# SPECIALISE fromRat'' :: Int -> Int -> Integer -> Integer -> Float,
889 Int -> Int -> Integer -> Integer -> Double #-}
890 fromRat'' :: RealFloat a => Int -> Int -> Integer -> Integer -> a
891 fromRat'' minEx@(I# me#) mantDigs@(I# md#) n d =
892 case integerLog2IsPowerOf2# d of
895 case integerLog2# n of
896 ln# | ln# ># (ld# +# me#) ->
898 then encodeFloat (n `shiftL` (I# (md# -# 1# -# ln#)))
899 (I# (ln# +# 1# -# ld# -# md#))
900 else let n' = n `shiftR` (I# (ln# +# 1# -# md#))
901 n'' = case roundingMode# n (ln# -# md#) of
904 _ -> case fromInteger n' .&. (1 :: Int) of
907 in encodeFloat n'' (I# (ln# -# ld# +# 1# -# md#))
909 case ld# +# (me# -# md#) of
910 ld'# | ld'# ># (ln# +# 1#) -> encodeFloat 0 0
911 | ld'# ==# (ln# +# 1#) ->
912 case integerLog2IsPowerOf2# n of
913 (# _, 0# #) -> encodeFloat 0 0
914 (# _, _ #) -> encodeFloat 1 (minEx - mantDigs)
916 encodeFloat n (I# ((me# -# md#) -# ld'#))
918 let n' = n `shiftR` (I# ld'#)
919 in case roundingMode# n (ld'# -# 1#) of
920 0# -> encodeFloat n' (minEx - mantDigs)
921 1# -> if fromInteger n' .&. (1 :: Int) == 0
922 then encodeFloat n' (minEx-mantDigs)
923 else encodeFloat (n' + 1) (minEx-mantDigs)
924 _ -> encodeFloat (n' + 1) (minEx-mantDigs)
926 let ln = I# (integerLog2# n)
928 p0 = max minEx (ln - ld)
930 | p0 < mantDigs = (n `shiftL` (mantDigs - p0), d)
931 | p0 == mantDigs = (n, d)
932 | otherwise = (n, d `shiftL` (p0 - mantDigs))
934 | p <= minEx-mantDigs = (p,a,b)
935 | a < (b `shiftL` (mantDigs-1)) = (p-1, a `shiftL` 1, b)
936 | (b `shiftL` mantDigs) <= a = (p+1, a, b `shiftL` 1)
937 | otherwise = (p, a, b)
938 (p', n'', d'') = scale (p0-mantDigs) n' d'
939 rdq = case n'' `quotRem` d'' of
940 (q,r) -> case compare (r `shiftL` 1) d'' of
942 EQ -> if fromInteger q .&. (1 :: Int) == 0
945 in encodeFloat rdq p'
949 %*********************************************************
951 \subsection{Floating point numeric primops}
953 %*********************************************************
955 Definitions of the boxed PrimOps; these will be
956 used in the case of partial applications, etc.
959 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
960 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
961 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
962 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
963 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
965 negateFloat :: Float -> Float
966 negateFloat (F# x) = F# (negateFloat# x)
968 gtFloat, geFloat, eqFloat, neFloat, ltFloat, leFloat :: Float -> Float -> Bool
969 gtFloat (F# x) (F# y) = gtFloat# x y
970 geFloat (F# x) (F# y) = geFloat# x y
971 eqFloat (F# x) (F# y) = eqFloat# x y
972 neFloat (F# x) (F# y) = neFloat# x y
973 ltFloat (F# x) (F# y) = ltFloat# x y
974 leFloat (F# x) (F# y) = leFloat# x y
976 expFloat, logFloat, sqrtFloat :: Float -> Float
977 sinFloat, cosFloat, tanFloat :: Float -> Float
978 asinFloat, acosFloat, atanFloat :: Float -> Float
979 sinhFloat, coshFloat, tanhFloat :: Float -> Float
980 expFloat (F# x) = F# (expFloat# x)
981 logFloat (F# x) = F# (logFloat# x)
982 sqrtFloat (F# x) = F# (sqrtFloat# x)
983 sinFloat (F# x) = F# (sinFloat# x)
984 cosFloat (F# x) = F# (cosFloat# x)
985 tanFloat (F# x) = F# (tanFloat# x)
986 asinFloat (F# x) = F# (asinFloat# x)
987 acosFloat (F# x) = F# (acosFloat# x)
988 atanFloat (F# x) = F# (atanFloat# x)
989 sinhFloat (F# x) = F# (sinhFloat# x)
990 coshFloat (F# x) = F# (coshFloat# x)
991 tanhFloat (F# x) = F# (tanhFloat# x)
993 powerFloat :: Float -> Float -> Float
994 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
996 -- definitions of the boxed PrimOps; these will be
997 -- used in the case of partial applications, etc.
999 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
1000 plusDouble (D# x) (D# y) = D# (x +## y)
1001 minusDouble (D# x) (D# y) = D# (x -## y)
1002 timesDouble (D# x) (D# y) = D# (x *## y)
1003 divideDouble (D# x) (D# y) = D# (x /## y)
1005 negateDouble :: Double -> Double
1006 negateDouble (D# x) = D# (negateDouble# x)
1008 gtDouble, geDouble, eqDouble, neDouble, leDouble, ltDouble :: Double -> Double -> Bool
1009 gtDouble (D# x) (D# y) = x >## y
1010 geDouble (D# x) (D# y) = x >=## y
1011 eqDouble (D# x) (D# y) = x ==## y
1012 neDouble (D# x) (D# y) = x /=## y
1013 ltDouble (D# x) (D# y) = x <## y
1014 leDouble (D# x) (D# y) = x <=## y
1016 double2Float :: Double -> Float
1017 double2Float (D# x) = F# (double2Float# x)
1019 float2Double :: Float -> Double
1020 float2Double (F# x) = D# (float2Double# x)
1022 expDouble, logDouble, sqrtDouble :: Double -> Double
1023 sinDouble, cosDouble, tanDouble :: Double -> Double
1024 asinDouble, acosDouble, atanDouble :: Double -> Double
1025 sinhDouble, coshDouble, tanhDouble :: Double -> Double
1026 expDouble (D# x) = D# (expDouble# x)
1027 logDouble (D# x) = D# (logDouble# x)
1028 sqrtDouble (D# x) = D# (sqrtDouble# x)
1029 sinDouble (D# x) = D# (sinDouble# x)
1030 cosDouble (D# x) = D# (cosDouble# x)
1031 tanDouble (D# x) = D# (tanDouble# x)
1032 asinDouble (D# x) = D# (asinDouble# x)
1033 acosDouble (D# x) = D# (acosDouble# x)
1034 atanDouble (D# x) = D# (atanDouble# x)
1035 sinhDouble (D# x) = D# (sinhDouble# x)
1036 coshDouble (D# x) = D# (coshDouble# x)
1037 tanhDouble (D# x) = D# (tanhDouble# x)
1039 powerDouble :: Double -> Double -> Double
1040 powerDouble (D# x) (D# y) = D# (x **## y)
1044 foreign import ccall unsafe "isFloatNaN" isFloatNaN :: Float -> Int
1045 foreign import ccall unsafe "isFloatInfinite" isFloatInfinite :: Float -> Int
1046 foreign import ccall unsafe "isFloatDenormalized" isFloatDenormalized :: Float -> Int
1047 foreign import ccall unsafe "isFloatNegativeZero" isFloatNegativeZero :: Float -> Int
1050 foreign import ccall unsafe "isDoubleNaN" isDoubleNaN :: Double -> Int
1051 foreign import ccall unsafe "isDoubleInfinite" isDoubleInfinite :: Double -> Int
1052 foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Double -> Int
1053 foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
1056 %*********************************************************
1058 \subsection{Coercion rules}
1060 %*********************************************************
1064 "fromIntegral/Int->Float" fromIntegral = int2Float
1065 "fromIntegral/Int->Double" fromIntegral = int2Double
1066 "realToFrac/Float->Float" realToFrac = id :: Float -> Float
1067 "realToFrac/Float->Double" realToFrac = float2Double
1068 "realToFrac/Double->Float" realToFrac = double2Float
1069 "realToFrac/Double->Double" realToFrac = id :: Double -> Double
1070 "realToFrac/Int->Double" realToFrac = int2Double -- See Note [realToFrac int-to-float]
1071 "realToFrac/Int->Float" realToFrac = int2Float -- ..ditto
1075 Note [realToFrac int-to-float]
1076 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1077 Don found that the RULES for realToFrac/Int->Double and simliarly
1078 Float made a huge difference to some stream-fusion programs. Here's
1081 import Data.Array.Vector
1086 let c = replicateU n (2::Double)
1087 a = mapU realToFrac (enumFromToU 0 (n-1) ) :: UArr Double
1088 print (sumU (zipWithU (*) c a))
1090 Without the RULE we get this loop body:
1092 case $wtoRational sc_sY4 of ww_aM7 { (# ww1_aM9, ww2_aMa #) ->
1093 case $wfromRat ww1_aM9 ww2_aMa of tpl_X1P { D# ipv_sW3 ->
1097 (+## sc2_sY6 (*## 2.0 ipv_sW3))
1104 (+## sc2_sXV (*## 2.0 (int2Double# sc_sXT)))
1106 The running time of the program goes from 120 seconds to 0.198 seconds
1107 with the native backend, and 0.143 seconds with the C backend.
1109 A few more details in Trac #2251, and the patch message
1110 "Add RULES for realToFrac from Int".
1112 %*********************************************************
1116 %*********************************************************
1119 showSignedFloat :: (RealFloat a)
1120 => (a -> ShowS) -- ^ a function that can show unsigned values
1121 -> Int -- ^ the precedence of the enclosing context
1122 -> a -- ^ the value to show
1124 showSignedFloat showPos p x
1125 | x < 0 || isNegativeZero x
1126 = showParen (p > 6) (showChar '-' . showPos (-x))
1127 | otherwise = showPos x
1130 We need to prevent over/underflow of the exponent in encodeFloat when
1131 called from scaleFloat, hence we clamp the scaling parameter.
1132 We must have a large enough range to cover the maximum difference of
1133 exponents returned by decodeFloat.
1135 clamp :: Int -> Int -> Int
1136 clamp bd k = max (-bd) (min bd k)