Update Trac ticket URLs to point to GitLab
[ghc.git] / libraries / base / GHC / Float.hs
1 {-# LANGUAGE Trustworthy #-}
2 {-# LANGUAGE CPP
3 , GHCForeignImportPrim
4 , NoImplicitPrelude
5 , MagicHash
6 , UnboxedTuples
7 , UnliftedFFITypes
8 #-}
9 {-# LANGUAGE CApiFFI #-}
10 -- We believe we could deorphan this module, by moving lots of things
11 -- around, but we haven't got there yet:
12 {-# OPTIONS_GHC -Wno-orphans #-}
13 {-# OPTIONS_HADDOCK not-home #-}
14
15 -----------------------------------------------------------------------------
16 -- |
17 -- Module : GHC.Float
18 -- Copyright : (c) The University of Glasgow 1994-2002
19 -- Portions obtained from hbc (c) Lennart Augusstson
20 -- License : see libraries/base/LICENSE
21 --
22 -- Maintainer : cvs-ghc@haskell.org
23 -- Stability : internal
24 -- Portability : non-portable (GHC Extensions)
25 --
26 -- The types 'Float' and 'Double', the classes 'Floating' and 'RealFloat' and
27 -- casting between Word32 and Float and Word64 and Double.
28 --
29 -----------------------------------------------------------------------------
30
31 #include "ieee-flpt.h"
32 #include "MachDeps.h"
33
34 module GHC.Float
35 ( module GHC.Float
36 , Float(..), Double(..), Float#, Double#
37 , double2Int, int2Double, float2Int, int2Float
38
39 -- * Monomorphic equality operators
40 -- | See GHC.Classes#matching_overloaded_methods_in_rules
41 , eqFloat, eqDouble
42 ) where
43
44 import Data.Maybe
45
46 import Data.Bits
47 import GHC.Base
48 import GHC.List
49 import GHC.Enum
50 import GHC.Show
51 import GHC.Num
52 import GHC.Real
53 import GHC.Word
54 import GHC.Arr
55 import GHC.Float.RealFracMethods
56 import GHC.Float.ConversionUtils
57 import GHC.Integer.Logarithms ( integerLogBase# )
58 import GHC.Integer.Logarithms.Internals
59
60 infixr 8 **
61
62 ------------------------------------------------------------------------
63 -- Standard numeric classes
64 ------------------------------------------------------------------------
65
66 -- | Trigonometric and hyperbolic functions and related functions.
67 --
68 -- The Haskell Report defines no laws for 'Floating'. However, @('+')@, @('*')@
69 -- and 'exp' are customarily expected to define an exponential field and have
70 -- the following properties:
71 --
72 -- * @exp (a + b)@ = @exp a * exp b@
73 -- * @exp (fromInteger 0)@ = @fromInteger 1@
74 --
75 class (Fractional a) => Floating a where
76 pi :: a
77 exp, log, sqrt :: a -> a
78 (**), logBase :: a -> a -> a
79 sin, cos, tan :: a -> a
80 asin, acos, atan :: a -> a
81 sinh, cosh, tanh :: a -> a
82 asinh, acosh, atanh :: a -> a
83
84 -- | @'log1p' x@ computes @'log' (1 + x)@, but provides more precise
85 -- results for small (absolute) values of @x@ if possible.
86 --
87 -- @since 4.9.0.0
88 log1p :: a -> a
89
90 -- | @'expm1' x@ computes @'exp' x - 1@, but provides more precise
91 -- results for small (absolute) values of @x@ if possible.
92 --
93 -- @since 4.9.0.0
94 expm1 :: a -> a
95
96 -- | @'log1pexp' x@ computes @'log' (1 + 'exp' x)@, but provides more
97 -- precise results if possible.
98 --
99 -- Examples:
100 --
101 -- * if @x@ is a large negative number, @'log' (1 + 'exp' x)@ will be
102 -- imprecise for the reasons given in 'log1p'.
103 --
104 -- * if @'exp' x@ is close to @-1@, @'log' (1 + 'exp' x)@ will be
105 -- imprecise for the reasons given in 'expm1'.
106 --
107 -- @since 4.9.0.0
108 log1pexp :: a -> a
109
110 -- | @'log1mexp' x@ computes @'log' (1 - 'exp' x)@, but provides more
111 -- precise results if possible.
112 --
113 -- Examples:
114 --
115 -- * if @x@ is a large negative number, @'log' (1 - 'exp' x)@ will be
116 -- imprecise for the reasons given in 'log1p'.
117 --
118 -- * if @'exp' x@ is close to @1@, @'log' (1 - 'exp' x)@ will be
119 -- imprecise for the reasons given in 'expm1'.
120 --
121 -- @since 4.9.0.0
122 log1mexp :: a -> a
123
124 {-# INLINE (**) #-}
125 {-# INLINE logBase #-}
126 {-# INLINE sqrt #-}
127 {-# INLINE tan #-}
128 {-# INLINE tanh #-}
129 x ** y = exp (log x * y)
130 logBase x y = log y / log x
131 sqrt x = x ** 0.5
132 tan x = sin x / cos x
133 tanh x = sinh x / cosh x
134
135 {-# INLINE log1p #-}
136 {-# INLINE expm1 #-}
137 {-# INLINE log1pexp #-}
138 {-# INLINE log1mexp #-}
139 log1p x = log (1 + x)
140 expm1 x = exp x - 1
141 log1pexp x = log1p (exp x)
142 log1mexp x = log1p (negate (exp x))
143
144 -- | Efficient, machine-independent access to the components of a
145 -- floating-point number.
146 class (RealFrac a, Floating a) => RealFloat a where
147 -- | a constant function, returning the radix of the representation
148 -- (often @2@)
149 floatRadix :: a -> Integer
150 -- | a constant function, returning the number of digits of
151 -- 'floatRadix' in the significand
152 floatDigits :: a -> Int
153 -- | a constant function, returning the lowest and highest values
154 -- the exponent may assume
155 floatRange :: a -> (Int,Int)
156 -- | The function 'decodeFloat' applied to a real floating-point
157 -- number returns the significand expressed as an 'Integer' and an
158 -- appropriately scaled exponent (an 'Int'). If @'decodeFloat' x@
159 -- yields @(m,n)@, then @x@ is equal in value to @m*b^^n@, where @b@
160 -- is the floating-point radix, and furthermore, either @m@ and @n@
161 -- are both zero or else @b^(d-1) <= 'abs' m < b^d@, where @d@ is
162 -- the value of @'floatDigits' x@.
163 -- In particular, @'decodeFloat' 0 = (0,0)@. If the type
164 -- contains a negative zero, also @'decodeFloat' (-0.0) = (0,0)@.
165 -- /The result of/ @'decodeFloat' x@ /is unspecified if either of/
166 -- @'isNaN' x@ /or/ @'isInfinite' x@ /is/ 'True'.
167 decodeFloat :: a -> (Integer,Int)
168 -- | 'encodeFloat' performs the inverse of 'decodeFloat' in the
169 -- sense that for finite @x@ with the exception of @-0.0@,
170 -- @'Prelude.uncurry' 'encodeFloat' ('decodeFloat' x) = x@.
171 -- @'encodeFloat' m n@ is one of the two closest representable
172 -- floating-point numbers to @m*b^^n@ (or @&#177;Infinity@ if overflow
173 -- occurs); usually the closer, but if @m@ contains too many bits,
174 -- the result may be rounded in the wrong direction.
175 encodeFloat :: Integer -> Int -> a
176 -- | 'exponent' corresponds to the second component of 'decodeFloat'.
177 -- @'exponent' 0 = 0@ and for finite nonzero @x@,
178 -- @'exponent' x = snd ('decodeFloat' x) + 'floatDigits' x@.
179 -- If @x@ is a finite floating-point number, it is equal in value to
180 -- @'significand' x * b ^^ 'exponent' x@, where @b@ is the
181 -- floating-point radix.
182 -- The behaviour is unspecified on infinite or @NaN@ values.
183 exponent :: a -> Int
184 -- | The first component of 'decodeFloat', scaled to lie in the open
185 -- interval (@-1@,@1@), either @0.0@ or of absolute value @>= 1\/b@,
186 -- where @b@ is the floating-point radix.
187 -- The behaviour is unspecified on infinite or @NaN@ values.
188 significand :: a -> a
189 -- | multiplies a floating-point number by an integer power of the radix
190 scaleFloat :: Int -> a -> a
191 -- | 'True' if the argument is an IEEE \"not-a-number\" (NaN) value
192 isNaN :: a -> Bool
193 -- | 'True' if the argument is an IEEE infinity or negative infinity
194 isInfinite :: a -> Bool
195 -- | 'True' if the argument is too small to be represented in
196 -- normalized format
197 isDenormalized :: a -> Bool
198 -- | 'True' if the argument is an IEEE negative zero
199 isNegativeZero :: a -> Bool
200 -- | 'True' if the argument is an IEEE floating point number
201 isIEEE :: a -> Bool
202 -- | a version of arctangent taking two real floating-point arguments.
203 -- For real floating @x@ and @y@, @'atan2' y x@ computes the angle
204 -- (from the positive x-axis) of the vector from the origin to the
205 -- point @(x,y)@. @'atan2' y x@ returns a value in the range [@-pi@,
206 -- @pi@]. It follows the Common Lisp semantics for the origin when
207 -- signed zeroes are supported. @'atan2' y 1@, with @y@ in a type
208 -- that is 'RealFloat', should return the same value as @'atan' y@.
209 -- A default definition of 'atan2' is provided, but implementors
210 -- can provide a more accurate implementation.
211 atan2 :: a -> a -> a
212
213
214 exponent x = if m == 0 then 0 else n + floatDigits x
215 where (m,n) = decodeFloat x
216
217 significand x = encodeFloat m (negate (floatDigits x))
218 where (m,_) = decodeFloat x
219
220 scaleFloat 0 x = x
221 scaleFloat k x
222 | isFix = x
223 | otherwise = encodeFloat m (n + clamp b k)
224 where (m,n) = decodeFloat x
225 (l,h) = floatRange x
226 d = floatDigits x
227 b = h - l + 4*d
228 -- n+k may overflow, which would lead
229 -- to wrong results, hence we clamp the
230 -- scaling parameter.
231 -- If n + k would be larger than h,
232 -- n + clamp b k must be too, simliar
233 -- for smaller than l - d.
234 -- Add a little extra to keep clear
235 -- from the boundary cases.
236 isFix = x == 0 || isNaN x || isInfinite x
237
238 atan2 y x
239 | x > 0 = atan (y/x)
240 | x == 0 && y > 0 = pi/2
241 | x < 0 && y > 0 = pi + atan (y/x)
242 |(x <= 0 && y < 0) ||
243 (x < 0 && isNegativeZero y) ||
244 (isNegativeZero x && isNegativeZero y)
245 = -atan2 (-y) x
246 | y == 0 && (x < 0 || isNegativeZero x)
247 = pi -- must be after the previous test on zero y
248 | x==0 && y==0 = y -- must be after the other double zero tests
249 | otherwise = x + y -- x or y is a NaN, return a NaN (via +)
250
251 ------------------------------------------------------------------------
252 -- Float
253 ------------------------------------------------------------------------
254
255 -- | @since 2.01
256 -- Note that due to the presence of @NaN@, not all elements of 'Float' have an
257 -- additive inverse.
258 --
259 -- >>> 0/0 + (negate 0/0 :: Float)
260 -- NaN
261 --
262 -- Also note that due to the presence of -0, `Float`'s 'Num' instance doesn't
263 -- have an additive identity
264 --
265 -- >>> 0 + (-0 :: Float)
266 -- 0.0
267 instance Num Float where
268 (+) x y = plusFloat x y
269 (-) x y = minusFloat x y
270 negate x = negateFloat x
271 (*) x y = timesFloat x y
272 abs x = fabsFloat x
273 signum x | x > 0 = 1
274 | x < 0 = negateFloat 1
275 | otherwise = x -- handles 0.0, (-0.0), and NaN
276
277 {-# INLINE fromInteger #-}
278 fromInteger i = F# (floatFromInteger i)
279
280 -- | @since 2.01
281 instance Real Float where
282 toRational (F# x#) =
283 case decodeFloat_Int# x# of
284 (# m#, e# #)
285 | isTrue# (e# >=# 0#) ->
286 (smallInteger m# `shiftLInteger` e#) :% 1
287 | isTrue# ((int2Word# m# `and#` 1##) `eqWord#` 0##) ->
288 case elimZerosInt# m# (negateInt# e#) of
289 (# n, d# #) -> n :% shiftLInteger 1 d#
290 | otherwise ->
291 smallInteger m# :% shiftLInteger 1 (negateInt# e#)
292
293 -- | @since 2.01
294 -- Note that due to the presence of @NaN@, not all elements of 'Float' have an
295 -- multiplicative inverse.
296 --
297 -- >>> 0/0 * (recip 0/0 :: Float)
298 -- NaN
299 instance Fractional Float where
300 (/) x y = divideFloat x y
301 {-# INLINE fromRational #-}
302 fromRational (n:%d) = rationalToFloat n d
303 recip x = 1.0 / x
304
305 rationalToFloat :: Integer -> Integer -> Float
306 {-# NOINLINE [1] rationalToFloat #-}
307 rationalToFloat n 0
308 | n == 0 = 0/0
309 | n < 0 = (-1)/0
310 | otherwise = 1/0
311 rationalToFloat n d
312 | n == 0 = encodeFloat 0 0
313 | n < 0 = -(fromRat'' minEx mantDigs (-n) d)
314 | otherwise = fromRat'' minEx mantDigs n d
315 where
316 minEx = FLT_MIN_EXP
317 mantDigs = FLT_MANT_DIG
318
319 -- RULES for Integer and Int
320 {-# RULES
321 "properFraction/Float->Integer" properFraction = properFractionFloatInteger
322 "truncate/Float->Integer" truncate = truncateFloatInteger
323 "floor/Float->Integer" floor = floorFloatInteger
324 "ceiling/Float->Integer" ceiling = ceilingFloatInteger
325 "round/Float->Integer" round = roundFloatInteger
326 "properFraction/Float->Int" properFraction = properFractionFloatInt
327 "truncate/Float->Int" truncate = float2Int
328 "floor/Float->Int" floor = floorFloatInt
329 "ceiling/Float->Int" ceiling = ceilingFloatInt
330 "round/Float->Int" round = roundFloatInt
331 #-}
332 -- | @since 2.01
333 instance RealFrac Float where
334
335 -- ceiling, floor, and truncate are all small
336 {-# INLINE [1] ceiling #-}
337 {-# INLINE [1] floor #-}
338 {-# INLINE [1] truncate #-}
339
340 -- We assume that FLT_RADIX is 2 so that we can use more efficient code
341 #if FLT_RADIX != 2
342 #error FLT_RADIX must be 2
343 #endif
344 properFraction (F# x#)
345 = case decodeFloat_Int# x# of
346 (# m#, n# #) ->
347 let m = I# m#
348 n = I# n#
349 in
350 if n >= 0
351 then (fromIntegral m * (2 ^ n), 0.0)
352 else let i = if m >= 0 then m `shiftR` negate n
353 else negate (negate m `shiftR` negate n)
354 f = m - (i `shiftL` negate n)
355 in (fromIntegral i, encodeFloat (fromIntegral f) n)
356
357 truncate x = case properFraction x of
358 (n,_) -> n
359
360 round x = case properFraction x of
361 (n,r) -> let
362 m = if r < 0.0 then n - 1 else n + 1
363 half_down = abs r - 0.5
364 in
365 case (compare half_down 0.0) of
366 LT -> n
367 EQ -> if even n then n else m
368 GT -> m
369
370 ceiling x = case properFraction x of
371 (n,r) -> if r > 0.0 then n + 1 else n
372
373 floor x = case properFraction x of
374 (n,r) -> if r < 0.0 then n - 1 else n
375
376 -- | @since 2.01
377 instance Floating Float where
378 pi = 3.141592653589793238
379 exp x = expFloat x
380 log x = logFloat x
381 sqrt x = sqrtFloat x
382 sin x = sinFloat x
383 cos x = cosFloat x
384 tan x = tanFloat x
385 asin x = asinFloat x
386 acos x = acosFloat x
387 atan x = atanFloat x
388 sinh x = sinhFloat x
389 cosh x = coshFloat x
390 tanh x = tanhFloat x
391 (**) x y = powerFloat x y
392 logBase x y = log y / log x
393
394 asinh x = asinhFloat x
395 acosh x = acoshFloat x
396 atanh x = atanhFloat x
397
398 log1p = log1pFloat
399 expm1 = expm1Float
400
401 log1mexp a
402 | a <= log 2 = log (negate (expm1Float a))
403 | otherwise = log1pFloat (negate (exp a))
404 {-# INLINE log1mexp #-}
405 log1pexp a
406 | a <= 18 = log1pFloat (exp a)
407 | a <= 100 = a + exp (negate a)
408 | otherwise = a
409 {-# INLINE log1pexp #-}
410
411 -- | @since 2.01
412 instance RealFloat Float where
413 floatRadix _ = FLT_RADIX -- from float.h
414 floatDigits _ = FLT_MANT_DIG -- ditto
415 floatRange _ = (FLT_MIN_EXP, FLT_MAX_EXP) -- ditto
416
417 decodeFloat (F# f#) = case decodeFloat_Int# f# of
418 (# i, e #) -> (smallInteger i, I# e)
419
420 encodeFloat i (I# e) = F# (encodeFloatInteger i e)
421
422 exponent x = case decodeFloat x of
423 (m,n) -> if m == 0 then 0 else n + floatDigits x
424
425 significand x = case decodeFloat x of
426 (m,_) -> encodeFloat m (negate (floatDigits x))
427
428 scaleFloat 0 x = x
429 scaleFloat k x
430 | isFix = x
431 | otherwise = case decodeFloat x of
432 (m,n) -> encodeFloat m (n + clamp bf k)
433 where bf = FLT_MAX_EXP - (FLT_MIN_EXP) + 4*FLT_MANT_DIG
434 isFix = x == 0 || isFloatFinite x == 0
435
436 isNaN x = 0 /= isFloatNaN x
437 isInfinite x = 0 /= isFloatInfinite x
438 isDenormalized x = 0 /= isFloatDenormalized x
439 isNegativeZero x = 0 /= isFloatNegativeZero x
440 isIEEE _ = True
441
442 -- | @since 2.01
443 instance Show Float where
444 showsPrec x = showSignedFloat showFloat x
445 showList = showList__ (showsPrec 0)
446
447 ------------------------------------------------------------------------
448 -- Double
449 ------------------------------------------------------------------------
450
451 -- | @since 2.01
452 -- Note that due to the presence of @NaN@, not all elements of 'Double' have an
453 -- additive inverse.
454 --
455 -- >>> 0/0 + (negate 0/0 :: Double)
456 -- NaN
457 --
458 -- Also note that due to the presence of -0, `Double`'s 'Num' instance doesn't
459 -- have an additive identity
460 --
461 -- >>> 0 + (-0 :: Double)
462 -- 0.0
463 instance Num Double where
464 (+) x y = plusDouble x y
465 (-) x y = minusDouble x y
466 negate x = negateDouble x
467 (*) x y = timesDouble x y
468 abs x = fabsDouble x
469 signum x | x > 0 = 1
470 | x < 0 = negateDouble 1
471 | otherwise = x -- handles 0.0, (-0.0), and NaN
472
473
474 {-# INLINE fromInteger #-}
475 fromInteger i = D# (doubleFromInteger i)
476
477
478 -- | @since 2.01
479 instance Real Double where
480 toRational (D# x#) =
481 case decodeDoubleInteger x# of
482 (# m, e# #)
483 | isTrue# (e# >=# 0#) ->
484 shiftLInteger m e# :% 1
485 | isTrue# ((integerToWord m `and#` 1##) `eqWord#` 0##) ->
486 case elimZerosInteger m (negateInt# e#) of
487 (# n, d# #) -> n :% shiftLInteger 1 d#
488 | otherwise ->
489 m :% shiftLInteger 1 (negateInt# e#)
490
491 -- | @since 2.01
492 -- Note that due to the presence of @NaN@, not all elements of 'Double' have an
493 -- multiplicative inverse.
494 --
495 -- >>> 0/0 * (recip 0/0 :: Double)
496 -- NaN
497 instance Fractional Double where
498 (/) x y = divideDouble x y
499 {-# INLINE fromRational #-}
500 fromRational (n:%d) = rationalToDouble n d
501 recip x = 1.0 / x
502
503 rationalToDouble :: Integer -> Integer -> Double
504 {-# NOINLINE [1] rationalToDouble #-}
505 rationalToDouble n 0
506 | n == 0 = 0/0
507 | n < 0 = (-1)/0
508 | otherwise = 1/0
509 rationalToDouble n d
510 | n == 0 = encodeFloat 0 0
511 | n < 0 = -(fromRat'' minEx mantDigs (-n) d)
512 | otherwise = fromRat'' minEx mantDigs n d
513 where
514 minEx = DBL_MIN_EXP
515 mantDigs = DBL_MANT_DIG
516
517 -- | @since 2.01
518 instance Floating Double where
519 pi = 3.141592653589793238
520 exp x = expDouble x
521 log x = logDouble x
522 sqrt x = sqrtDouble x
523 sin x = sinDouble x
524 cos x = cosDouble x
525 tan x = tanDouble x
526 asin x = asinDouble x
527 acos x = acosDouble x
528 atan x = atanDouble x
529 sinh x = sinhDouble x
530 cosh x = coshDouble x
531 tanh x = tanhDouble x
532 (**) x y = powerDouble x y
533 logBase x y = log y / log x
534
535 asinh x = asinhDouble x
536 acosh x = acoshDouble x
537 atanh x = atanhDouble x
538
539 log1p = log1pDouble
540 expm1 = expm1Double
541
542 log1mexp a
543 | a <= log 2 = log (negate (expm1Double a))
544 | otherwise = log1pDouble (negate (exp a))
545 {-# INLINE log1mexp #-}
546 log1pexp a
547 | a <= 18 = log1pDouble (exp a)
548 | a <= 100 = a + exp (negate a)
549 | otherwise = a
550 {-# INLINE log1pexp #-}
551
552 -- RULES for Integer and Int
553 {-# RULES
554 "properFraction/Double->Integer" properFraction = properFractionDoubleInteger
555 "truncate/Double->Integer" truncate = truncateDoubleInteger
556 "floor/Double->Integer" floor = floorDoubleInteger
557 "ceiling/Double->Integer" ceiling = ceilingDoubleInteger
558 "round/Double->Integer" round = roundDoubleInteger
559 "properFraction/Double->Int" properFraction = properFractionDoubleInt
560 "truncate/Double->Int" truncate = double2Int
561 "floor/Double->Int" floor = floorDoubleInt
562 "ceiling/Double->Int" ceiling = ceilingDoubleInt
563 "round/Double->Int" round = roundDoubleInt
564 #-}
565 -- | @since 2.01
566 instance RealFrac Double where
567
568 -- ceiling, floor, and truncate are all small
569 {-# INLINE [1] ceiling #-}
570 {-# INLINE [1] floor #-}
571 {-# INLINE [1] truncate #-}
572
573 properFraction x
574 = case (decodeFloat x) of { (m,n) ->
575 if n >= 0 then
576 (fromInteger m * 2 ^ n, 0.0)
577 else
578 case (quotRem m (2^(negate n))) of { (w,r) ->
579 (fromInteger w, encodeFloat r n)
580 }
581 }
582
583 truncate x = case properFraction x of
584 (n,_) -> n
585
586 round x = case properFraction x of
587 (n,r) -> let
588 m = if r < 0.0 then n - 1 else n + 1
589 half_down = abs r - 0.5
590 in
591 case (compare half_down 0.0) of
592 LT -> n
593 EQ -> if even n then n else m
594 GT -> m
595
596 ceiling x = case properFraction x of
597 (n,r) -> if r > 0.0 then n + 1 else n
598
599 floor x = case properFraction x of
600 (n,r) -> if r < 0.0 then n - 1 else n
601
602 -- | @since 2.01
603 instance RealFloat Double where
604 floatRadix _ = FLT_RADIX -- from float.h
605 floatDigits _ = DBL_MANT_DIG -- ditto
606 floatRange _ = (DBL_MIN_EXP, DBL_MAX_EXP) -- ditto
607
608 decodeFloat (D# x#)
609 = case decodeDoubleInteger x# of
610 (# i, j #) -> (i, I# j)
611
612 encodeFloat i (I# j) = D# (encodeDoubleInteger i j)
613
614 exponent x = case decodeFloat x of
615 (m,n) -> if m == 0 then 0 else n + floatDigits x
616
617 significand x = case decodeFloat x of
618 (m,_) -> encodeFloat m (negate (floatDigits x))
619
620 scaleFloat 0 x = x
621 scaleFloat k x
622 | isFix = x
623 | otherwise = case decodeFloat x of
624 (m,n) -> encodeFloat m (n + clamp bd k)
625 where bd = DBL_MAX_EXP - (DBL_MIN_EXP) + 4*DBL_MANT_DIG
626 isFix = x == 0 || isDoubleFinite x == 0
627
628 isNaN x = 0 /= isDoubleNaN x
629 isInfinite x = 0 /= isDoubleInfinite x
630 isDenormalized x = 0 /= isDoubleDenormalized x
631 isNegativeZero x = 0 /= isDoubleNegativeZero x
632 isIEEE _ = True
633
634 -- | @since 2.01
635 instance Show Double where
636 showsPrec x = showSignedFloat showFloat x
637 showList = showList__ (showsPrec 0)
638
639
640 ------------------------------------------------------------------------
641 -- Enum instances
642 ------------------------------------------------------------------------
643
644 {-
645 The @Enum@ instances for Floats and Doubles are slightly unusual.
646 The @toEnum@ function truncates numbers to Int. The definitions
647 of @enumFrom@ and @enumFromThen@ allow floats to be used in arithmetic
648 series: [0,0.1 .. 1.0]. However, roundoff errors make these somewhat
649 dubious. This example may have either 10 or 11 elements, depending on
650 how 0.1 is represented.
651
652 NOTE: The instances for Float and Double do not make use of the default
653 methods for @enumFromTo@ and @enumFromThenTo@, as these rely on there being
654 a `non-lossy' conversion to and from Ints. Instead we make use of the
655 1.2 default methods (back in the days when Enum had Ord as a superclass)
656 for these (@numericEnumFromTo@ and @numericEnumFromThenTo@ below.)
657 -}
658
659 -- | @since 2.01
660 instance Enum Float where
661 succ x = x + 1
662 pred x = x - 1
663 toEnum = int2Float
664 fromEnum = fromInteger . truncate -- may overflow
665 enumFrom = numericEnumFrom
666 enumFromTo = numericEnumFromTo
667 enumFromThen = numericEnumFromThen
668 enumFromThenTo = numericEnumFromThenTo
669
670 -- | @since 2.01
671 instance Enum Double where
672 succ x = x + 1
673 pred x = x - 1
674 toEnum = int2Double
675 fromEnum = fromInteger . truncate -- may overflow
676 enumFrom = numericEnumFrom
677 enumFromTo = numericEnumFromTo
678 enumFromThen = numericEnumFromThen
679 enumFromThenTo = numericEnumFromThenTo
680
681 ------------------------------------------------------------------------
682 -- Printing floating point
683 ------------------------------------------------------------------------
684
685 -- | Show a signed 'RealFloat' value to full precision
686 -- using standard decimal notation for arguments whose absolute value lies
687 -- between @0.1@ and @9,999,999@, and scientific notation otherwise.
688 showFloat :: (RealFloat a) => a -> ShowS
689 showFloat x = showString (formatRealFloat FFGeneric Nothing x)
690
691 -- These are the format types. This type is not exported.
692
693 data FFFormat = FFExponent | FFFixed | FFGeneric
694
695 -- This is just a compatibility stub, as the "alt" argument formerly
696 -- didn't exist.
697 formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
698 formatRealFloat fmt decs x = formatRealFloatAlt fmt decs False x
699
700 formatRealFloatAlt :: (RealFloat a) => FFFormat -> Maybe Int -> Bool -> a
701 -> String
702 formatRealFloatAlt fmt decs alt x
703 | isNaN x = "NaN"
704 | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
705 | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
706 | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
707 where
708 base = 10
709
710 doFmt format (is, e) =
711 let ds = map intToDigit is in
712 case format of
713 FFGeneric ->
714 doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
715 (is,e)
716 FFExponent ->
717 case decs of
718 Nothing ->
719 let show_e' = show (e-1) in
720 case ds of
721 "0" -> "0.0e0"
722 [d] -> d : ".0e" ++ show_e'
723 (d:ds') -> d : '.' : ds' ++ "e" ++ show_e'
724 [] -> errorWithoutStackTrace "formatRealFloat/doFmt/FFExponent: []"
725 Just d | d <= 0 ->
726 -- handle this case specifically since we need to omit the
727 -- decimal point as well (#15115).
728 -- Note that this handles negative precisions as well for consistency
729 -- (see #15509).
730 case is of
731 [0] -> "0e0"
732 _ ->
733 let
734 (ei,is') = roundTo base 1 is
735 n:_ = map intToDigit (if ei > 0 then init is' else is')
736 in n : 'e' : show (e-1+ei)
737 Just dec ->
738 let dec' = max dec 1 in
739 case is of
740 [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
741 _ ->
742 let
743 (ei,is') = roundTo base (dec'+1) is
744 (d:ds') = map intToDigit (if ei > 0 then init is' else is')
745 in
746 d:'.':ds' ++ 'e':show (e-1+ei)
747 FFFixed ->
748 let
749 mk0 ls = case ls of { "" -> "0" ; _ -> ls}
750 in
751 case decs of
752 Nothing
753 | e <= 0 -> "0." ++ replicate (-e) '0' ++ ds
754 | otherwise ->
755 let
756 f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
757 f n s "" = f (n-1) ('0':s) ""
758 f n s (r:rs) = f (n-1) (r:s) rs
759 in
760 f e "" ds
761 Just dec ->
762 let dec' = max dec 0 in
763 if e >= 0 then
764 let
765 (ei,is') = roundTo base (dec' + e) is
766 (ls,rs) = splitAt (e+ei) (map intToDigit is')
767 in
768 mk0 ls ++ (if null rs && not alt then "" else '.':rs)
769 else
770 let
771 (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
772 d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
773 in
774 d : (if null ds' && not alt then "" else '.':ds')
775
776
777 roundTo :: Int -> Int -> [Int] -> (Int,[Int])
778 roundTo base d is =
779 case f d True is of
780 x@(0,_) -> x
781 (1,xs) -> (1, 1:xs)
782 _ -> errorWithoutStackTrace "roundTo: bad Value"
783 where
784 b2 = base `quot` 2
785
786 f n _ [] = (0, replicate n 0)
787 f 0 e (x:xs) | x == b2 && e && all (== 0) xs = (0, []) -- Round to even when at exactly half the base
788 | otherwise = (if x >= b2 then 1 else 0, [])
789 f n _ (i:xs)
790 | i' == base = (1,0:ds)
791 | otherwise = (0,i':ds)
792 where
793 (c,ds) = f (n-1) (even i) xs
794 i' = c + i
795
796 -- Based on "Printing Floating-Point Numbers Quickly and Accurately"
797 -- by R.G. Burger and R.K. Dybvig in PLDI 96.
798 -- This version uses a much slower logarithm estimator. It should be improved.
799
800 -- | 'floatToDigits' takes a base and a non-negative 'RealFloat' number,
801 -- and returns a list of digits and an exponent.
802 -- In particular, if @x>=0@, and
803 --
804 -- > floatToDigits base x = ([d1,d2,...,dn], e)
805 --
806 -- then
807 --
808 -- (1) @n >= 1@
809 --
810 -- (2) @x = 0.d1d2...dn * (base**e)@
811 --
812 -- (3) @0 <= di <= base-1@
813
814 floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
815 floatToDigits _ 0 = ([0], 0)
816 floatToDigits base x =
817 let
818 (f0, e0) = decodeFloat x
819 (minExp0, _) = floatRange x
820 p = floatDigits x
821 b = floatRadix x
822 minExp = minExp0 - p -- the real minimum exponent
823 -- Haskell requires that f be adjusted so denormalized numbers
824 -- will have an impossibly low exponent. Adjust for this.
825 (f, e) =
826 let n = minExp - e0 in
827 if n > 0 then (f0 `quot` (expt b n), e0+n) else (f0, e0)
828 (r, s, mUp, mDn) =
829 if e >= 0 then
830 let be = expt b e in
831 if f == expt b (p-1) then
832 (f*be*b*2, 2*b, be*b, be) -- according to Burger and Dybvig
833 else
834 (f*be*2, 2, be, be)
835 else
836 if e > minExp && f == expt b (p-1) then
837 (f*b*2, expt b (-e+1)*2, b, 1)
838 else
839 (f*2, expt b (-e)*2, 1, 1)
840 k :: Int
841 k =
842 let
843 k0 :: Int
844 k0 =
845 if b == 2 && base == 10 then
846 -- logBase 10 2 is very slightly larger than 8651/28738
847 -- (about 5.3558e-10), so if log x >= 0, the approximation
848 -- k1 is too small, hence we add one and need one fixup step less.
849 -- If log x < 0, the approximation errs rather on the high side.
850 -- That is usually more than compensated for by ignoring the
851 -- fractional part of logBase 2 x, but when x is a power of 1/2
852 -- or slightly larger and the exponent is a multiple of the
853 -- denominator of the rational approximation to logBase 10 2,
854 -- k1 is larger than logBase 10 x. If k1 > 1 + logBase 10 x,
855 -- we get a leading zero-digit we don't want.
856 -- With the approximation 3/10, this happened for
857 -- 0.5^1030, 0.5^1040, ..., 0.5^1070 and values close above.
858 -- The approximation 8651/28738 guarantees k1 < 1 + logBase 10 x
859 -- for IEEE-ish floating point types with exponent fields
860 -- <= 17 bits and mantissae of several thousand bits, earlier
861 -- convergents to logBase 10 2 would fail for long double.
862 -- Using quot instead of div is a little faster and requires
863 -- fewer fixup steps for negative lx.
864 let lx = p - 1 + e0
865 k1 = (lx * 8651) `quot` 28738
866 in if lx >= 0 then k1 + 1 else k1
867 else
868 -- f :: Integer, log :: Float -> Float,
869 -- ceiling :: Float -> Int
870 ceiling ((log (fromInteger (f+1) :: Float) +
871 fromIntegral e * log (fromInteger b)) /
872 log (fromInteger base))
873 --WAS: fromInt e * log (fromInteger b))
874
875 fixup n =
876 if n >= 0 then
877 if r + mUp <= expt base n * s then n else fixup (n+1)
878 else
879 if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
880 in
881 fixup k0
882
883 gen ds rn sN mUpN mDnN =
884 let
885 (dn, rn') = (rn * base) `quotRem` sN
886 mUpN' = mUpN * base
887 mDnN' = mDnN * base
888 in
889 case (rn' < mDnN', rn' + mUpN' > sN) of
890 (True, False) -> dn : ds
891 (False, True) -> dn+1 : ds
892 (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
893 (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
894
895 rds =
896 if k >= 0 then
897 gen [] r (s * expt base k) mUp mDn
898 else
899 let bk = expt base (-k) in
900 gen [] (r * bk) s (mUp * bk) (mDn * bk)
901 in
902 (map fromIntegral (reverse rds), k)
903
904 ------------------------------------------------------------------------
905 -- Converting from a Rational to a RealFloa
906 ------------------------------------------------------------------------
907
908 {-
909 [In response to a request for documentation of how fromRational works,
910 Joe Fasel writes:] A quite reasonable request! This code was added to
911 the Prelude just before the 1.2 release, when Lennart, working with an
912 early version of hbi, noticed that (read . show) was not the identity
913 for floating-point numbers. (There was a one-bit error about half the
914 time.) The original version of the conversion function was in fact
915 simply a floating-point divide, as you suggest above. The new version
916 is, I grant you, somewhat denser.
917
918 Unfortunately, Joe's code doesn't work! Here's an example:
919
920 main = putStr (shows (1.82173691287639817263897126389712638972163e-300::Double) "\n")
921
922 This program prints
923 0.0000000000000000
924 instead of
925 1.8217369128763981e-300
926
927 Here's Joe's code:
928
929 \begin{pseudocode}
930 fromRat :: (RealFloat a) => Rational -> a
931 fromRat x = x'
932 where x' = f e
933
934 -- If the exponent of the nearest floating-point number to x
935 -- is e, then the significand is the integer nearest xb^(-e),
936 -- where b is the floating-point radix. We start with a good
937 -- guess for e, and if it is correct, the exponent of the
938 -- floating-point number we construct will again be e. If
939 -- not, one more iteration is needed.
940
941 f e = if e' == e then y else f e'
942 where y = encodeFloat (round (x * (1 % b)^^e)) e
943 (_,e') = decodeFloat y
944 b = floatRadix x'
945
946 -- We obtain a trial exponent by doing a floating-point
947 -- division of x's numerator by its denominator. The
948 -- result of this division may not itself be the ultimate
949 -- result, because of an accumulation of three rounding
950 -- errors.
951
952 (s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
953 / fromInteger (denominator x))
954 \end{pseudocode}
955
956 Now, here's Lennart's code (which works):
957 -}
958
959 -- | Converts a 'Rational' value into any type in class 'RealFloat'.
960 {-# RULES
961 "fromRat/Float" fromRat = (fromRational :: Rational -> Float)
962 "fromRat/Double" fromRat = (fromRational :: Rational -> Double)
963 #-}
964
965 {-# NOINLINE [1] fromRat #-}
966 fromRat :: (RealFloat a) => Rational -> a
967
968 -- Deal with special cases first, delegating the real work to fromRat'
969 fromRat (n :% 0) | n > 0 = 1/0 -- +Infinity
970 | n < 0 = -1/0 -- -Infinity
971 | otherwise = 0/0 -- NaN
972
973 fromRat (n :% d) | n > 0 = fromRat' (n :% d)
974 | n < 0 = - fromRat' ((-n) :% d)
975 | otherwise = encodeFloat 0 0 -- Zero
976
977 -- Conversion process:
978 -- Scale the rational number by the RealFloat base until
979 -- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).
980 -- Then round the rational to an Integer and encode it with the exponent
981 -- that we got from the scaling.
982 -- To speed up the scaling process we compute the log2 of the number to get
983 -- a first guess of the exponent.
984
985 fromRat' :: (RealFloat a) => Rational -> a
986 -- Invariant: argument is strictly positive
987 fromRat' x = r
988 where b = floatRadix r
989 p = floatDigits r
990 (minExp0, _) = floatRange r
991 minExp = minExp0 - p -- the real minimum exponent
992 xMax = toRational (expt b p)
993 p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp
994 -- if x = n/d and ln = integerLogBase b n, ld = integerLogBase b d,
995 -- then b^(ln-ld-1) < x < b^(ln-ld+1)
996 f = if p0 < 0 then 1 :% expt b (-p0) else expt b p0 :% 1
997 x0 = x / f
998 -- if ln - ld >= minExp0, then b^(p-1) < x0 < b^(p+1), so there's at most
999 -- one scaling step needed, otherwise, x0 < b^p and no scaling is needed
1000 (x', p') = if x0 >= xMax then (x0 / toRational b, p0+1) else (x0, p0)
1001 r = encodeFloat (round x') p'
1002
1003 -- Exponentiation with a cache for the most common numbers.
1004 minExpt, maxExpt :: Int
1005 minExpt = 0
1006 maxExpt = 1100
1007
1008 expt :: Integer -> Int -> Integer
1009 expt base n =
1010 if base == 2 && n >= minExpt && n <= maxExpt then
1011 expts!n
1012 else
1013 if base == 10 && n <= maxExpt10 then
1014 expts10!n
1015 else
1016 base^n
1017
1018 expts :: Array Int Integer
1019 expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
1020
1021 maxExpt10 :: Int
1022 maxExpt10 = 324
1023
1024 expts10 :: Array Int Integer
1025 expts10 = array (minExpt,maxExpt10) [(n,10^n) | n <- [minExpt .. maxExpt10]]
1026
1027 -- Compute the (floor of the) log of i in base b.
1028 -- Simplest way would be just divide i by b until it's smaller then b, but that would
1029 -- be very slow! We are just slightly more clever, except for base 2, where
1030 -- we take advantage of the representation of Integers.
1031 -- The general case could be improved by a lookup table for
1032 -- approximating the result by integerLog2 i / integerLog2 b.
1033 integerLogBase :: Integer -> Integer -> Int
1034 integerLogBase b i
1035 | i < b = 0
1036 | b == 2 = I# (integerLog2# i)
1037 | otherwise = I# (integerLogBase# b i)
1038
1039 {-
1040 Unfortunately, the old conversion code was awfully slow due to
1041 a) a slow integer logarithm
1042 b) repeated calculation of gcd's
1043
1044 For the case of Rational's coming from a Float or Double via toRational,
1045 we can exploit the fact that the denominator is a power of two, which for
1046 these brings a huge speedup since we need only shift and add instead
1047 of division.
1048
1049 The below is an adaption of fromRat' for the conversion to
1050 Float or Double exploiting the known floatRadix and avoiding
1051 divisions as much as possible.
1052 -}
1053
1054 {-# SPECIALISE fromRat'' :: Int -> Int -> Integer -> Integer -> Float,
1055 Int -> Int -> Integer -> Integer -> Double #-}
1056 fromRat'' :: RealFloat a => Int -> Int -> Integer -> Integer -> a
1057 -- Invariant: n and d strictly positive
1058 fromRat'' minEx@(I# me#) mantDigs@(I# md#) n d =
1059 case integerLog2IsPowerOf2# d of
1060 (# ld#, pw# #)
1061 | isTrue# (pw# ==# 0#) ->
1062 case integerLog2# n of
1063 ln# | isTrue# (ln# >=# (ld# +# me# -# 1#)) ->
1064 -- this means n/d >= 2^(minEx-1), i.e. we are guaranteed to get
1065 -- a normalised number, round to mantDigs bits
1066 if isTrue# (ln# <# md#)
1067 then encodeFloat n (I# (negateInt# ld#))
1068 else let n' = n `shiftR` (I# (ln# +# 1# -# md#))
1069 n'' = case roundingMode# n (ln# -# md#) of
1070 0# -> n'
1071 2# -> n' + 1
1072 _ -> case fromInteger n' .&. (1 :: Int) of
1073 0 -> n'
1074 _ -> n' + 1
1075 in encodeFloat n'' (I# (ln# -# ld# +# 1# -# md#))
1076 | otherwise ->
1077 -- n/d < 2^(minEx-1), a denorm or rounded to 2^(minEx-1)
1078 -- the exponent for encoding is always minEx-mantDigs
1079 -- so we must shift right by (minEx-mantDigs) - (-ld)
1080 case ld# +# (me# -# md#) of
1081 ld'# | isTrue# (ld'# <=# 0#) -> -- we would shift left, so we don't shift
1082 encodeFloat n (I# ((me# -# md#) -# ld'#))
1083 | isTrue# (ld'# <=# ln#) ->
1084 let n' = n `shiftR` (I# ld'#)
1085 in case roundingMode# n (ld'# -# 1#) of
1086 0# -> encodeFloat n' (minEx - mantDigs)
1087 1# -> if fromInteger n' .&. (1 :: Int) == 0
1088 then encodeFloat n' (minEx-mantDigs)
1089 else encodeFloat (n' + 1) (minEx-mantDigs)
1090 _ -> encodeFloat (n' + 1) (minEx-mantDigs)
1091 | isTrue# (ld'# ># (ln# +# 1#)) -> encodeFloat 0 0 -- result of shift < 0.5
1092 | otherwise -> -- first bit of n shifted to 0.5 place
1093 case integerLog2IsPowerOf2# n of
1094 (# _, 0# #) -> encodeFloat 0 0 -- round to even
1095 (# _, _ #) -> encodeFloat 1 (minEx - mantDigs)
1096 | otherwise ->
1097 let ln = I# (integerLog2# n)
1098 ld = I# ld#
1099 -- 2^(ln-ld-1) < n/d < 2^(ln-ld+1)
1100 p0 = max minEx (ln - ld)
1101 (n', d')
1102 | p0 < mantDigs = (n `shiftL` (mantDigs - p0), d)
1103 | p0 == mantDigs = (n, d)
1104 | otherwise = (n, d `shiftL` (p0 - mantDigs))
1105 -- if ln-ld < minEx, then n'/d' < 2^mantDigs, else
1106 -- 2^(mantDigs-1) < n'/d' < 2^(mantDigs+1) and we
1107 -- may need one scaling step
1108 scale p a b
1109 | (b `shiftL` mantDigs) <= a = (p+1, a, b `shiftL` 1)
1110 | otherwise = (p, a, b)
1111 (p', n'', d'') = scale (p0-mantDigs) n' d'
1112 -- n''/d'' < 2^mantDigs and p' == minEx-mantDigs or n''/d'' >= 2^(mantDigs-1)
1113 rdq = case n'' `quotRem` d'' of
1114 (q,r) -> case compare (r `shiftL` 1) d'' of
1115 LT -> q
1116 EQ -> if fromInteger q .&. (1 :: Int) == 0
1117 then q else q+1
1118 GT -> q+1
1119 in encodeFloat rdq p'
1120
1121 ------------------------------------------------------------------------
1122 -- Floating point numeric primops
1123 ------------------------------------------------------------------------
1124
1125 -- Definitions of the boxed PrimOps; these will be
1126 -- used in the case of partial applications, etc.
1127
1128 plusFloat, minusFloat, timesFloat, divideFloat :: Float -> Float -> Float
1129 plusFloat (F# x) (F# y) = F# (plusFloat# x y)
1130 minusFloat (F# x) (F# y) = F# (minusFloat# x y)
1131 timesFloat (F# x) (F# y) = F# (timesFloat# x y)
1132 divideFloat (F# x) (F# y) = F# (divideFloat# x y)
1133
1134 negateFloat :: Float -> Float
1135 negateFloat (F# x) = F# (negateFloat# x)
1136
1137 gtFloat, geFloat, ltFloat, leFloat :: Float -> Float -> Bool
1138 gtFloat (F# x) (F# y) = isTrue# (gtFloat# x y)
1139 geFloat (F# x) (F# y) = isTrue# (geFloat# x y)
1140 ltFloat (F# x) (F# y) = isTrue# (ltFloat# x y)
1141 leFloat (F# x) (F# y) = isTrue# (leFloat# x y)
1142
1143 expFloat, logFloat, sqrtFloat, fabsFloat :: Float -> Float
1144 sinFloat, cosFloat, tanFloat :: Float -> Float
1145 asinFloat, acosFloat, atanFloat :: Float -> Float
1146 sinhFloat, coshFloat, tanhFloat :: Float -> Float
1147 asinhFloat, acoshFloat, atanhFloat :: Float -> Float
1148 expFloat (F# x) = F# (expFloat# x)
1149 logFloat (F# x) = F# (logFloat# x)
1150 sqrtFloat (F# x) = F# (sqrtFloat# x)
1151 fabsFloat (F# x) = F# (fabsFloat# x)
1152 sinFloat (F# x) = F# (sinFloat# x)
1153 cosFloat (F# x) = F# (cosFloat# x)
1154 tanFloat (F# x) = F# (tanFloat# x)
1155 asinFloat (F# x) = F# (asinFloat# x)
1156 acosFloat (F# x) = F# (acosFloat# x)
1157 atanFloat (F# x) = F# (atanFloat# x)
1158 sinhFloat (F# x) = F# (sinhFloat# x)
1159 coshFloat (F# x) = F# (coshFloat# x)
1160 tanhFloat (F# x) = F# (tanhFloat# x)
1161 asinhFloat (F# x) = F# (asinhFloat# x)
1162 acoshFloat (F# x) = F# (acoshFloat# x)
1163 atanhFloat (F# x) = F# (atanhFloat# x)
1164
1165 powerFloat :: Float -> Float -> Float
1166 powerFloat (F# x) (F# y) = F# (powerFloat# x y)
1167
1168 -- definitions of the boxed PrimOps; these will be
1169 -- used in the case of partial applications, etc.
1170
1171 plusDouble, minusDouble, timesDouble, divideDouble :: Double -> Double -> Double
1172 plusDouble (D# x) (D# y) = D# (x +## y)
1173 minusDouble (D# x) (D# y) = D# (x -## y)
1174 timesDouble (D# x) (D# y) = D# (x *## y)
1175 divideDouble (D# x) (D# y) = D# (x /## y)
1176
1177 negateDouble :: Double -> Double
1178 negateDouble (D# x) = D# (negateDouble# x)
1179
1180 gtDouble, geDouble, leDouble, ltDouble :: Double -> Double -> Bool
1181 gtDouble (D# x) (D# y) = isTrue# (x >## y)
1182 geDouble (D# x) (D# y) = isTrue# (x >=## y)
1183 ltDouble (D# x) (D# y) = isTrue# (x <## y)
1184 leDouble (D# x) (D# y) = isTrue# (x <=## y)
1185
1186 double2Float :: Double -> Float
1187 double2Float (D# x) = F# (double2Float# x)
1188
1189 float2Double :: Float -> Double
1190 float2Double (F# x) = D# (float2Double# x)
1191
1192 expDouble, logDouble, sqrtDouble, fabsDouble :: Double -> Double
1193 sinDouble, cosDouble, tanDouble :: Double -> Double
1194 asinDouble, acosDouble, atanDouble :: Double -> Double
1195 sinhDouble, coshDouble, tanhDouble :: Double -> Double
1196 asinhDouble, acoshDouble, atanhDouble :: Double -> Double
1197 expDouble (D# x) = D# (expDouble# x)
1198 logDouble (D# x) = D# (logDouble# x)
1199 sqrtDouble (D# x) = D# (sqrtDouble# x)
1200 fabsDouble (D# x) = D# (fabsDouble# x)
1201 sinDouble (D# x) = D# (sinDouble# x)
1202 cosDouble (D# x) = D# (cosDouble# x)
1203 tanDouble (D# x) = D# (tanDouble# x)
1204 asinDouble (D# x) = D# (asinDouble# x)
1205 acosDouble (D# x) = D# (acosDouble# x)
1206 atanDouble (D# x) = D# (atanDouble# x)
1207 sinhDouble (D# x) = D# (sinhDouble# x)
1208 coshDouble (D# x) = D# (coshDouble# x)
1209 tanhDouble (D# x) = D# (tanhDouble# x)
1210 asinhDouble (D# x) = D# (asinhDouble# x)
1211 acoshDouble (D# x) = D# (acoshDouble# x)
1212 atanhDouble (D# x) = D# (atanhDouble# x)
1213
1214 powerDouble :: Double -> Double -> Double
1215 powerDouble (D# x) (D# y) = D# (x **## y)
1216
1217 foreign import ccall unsafe "isFloatNaN" isFloatNaN :: Float -> Int
1218 foreign import ccall unsafe "isFloatInfinite" isFloatInfinite :: Float -> Int
1219 foreign import ccall unsafe "isFloatDenormalized" isFloatDenormalized :: Float -> Int
1220 foreign import ccall unsafe "isFloatNegativeZero" isFloatNegativeZero :: Float -> Int
1221 foreign import ccall unsafe "isFloatFinite" isFloatFinite :: Float -> Int
1222
1223 foreign import ccall unsafe "isDoubleNaN" isDoubleNaN :: Double -> Int
1224 foreign import ccall unsafe "isDoubleInfinite" isDoubleInfinite :: Double -> Int
1225 foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Double -> Int
1226 foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
1227 foreign import ccall unsafe "isDoubleFinite" isDoubleFinite :: Double -> Int
1228
1229
1230 ------------------------------------------------------------------------
1231 -- libm imports for extended floating
1232 ------------------------------------------------------------------------
1233 foreign import capi unsafe "math.h log1p" log1pDouble :: Double -> Double
1234 foreign import capi unsafe "math.h expm1" expm1Double :: Double -> Double
1235 foreign import capi unsafe "math.h log1pf" log1pFloat :: Float -> Float
1236 foreign import capi unsafe "math.h expm1f" expm1Float :: Float -> Float
1237
1238
1239 ------------------------------------------------------------------------
1240 -- Coercion rules
1241 ------------------------------------------------------------------------
1242
1243 word2Double :: Word -> Double
1244 word2Double (W# w) = D# (word2Double# w)
1245
1246 word2Float :: Word -> Float
1247 word2Float (W# w) = F# (word2Float# w)
1248
1249 {-# RULES
1250 "fromIntegral/Int->Float" fromIntegral = int2Float
1251 "fromIntegral/Int->Double" fromIntegral = int2Double
1252 "fromIntegral/Word->Float" fromIntegral = word2Float
1253 "fromIntegral/Word->Double" fromIntegral = word2Double
1254 "realToFrac/Float->Float" realToFrac = id :: Float -> Float
1255 "realToFrac/Float->Double" realToFrac = float2Double
1256 "realToFrac/Double->Float" realToFrac = double2Float
1257 "realToFrac/Double->Double" realToFrac = id :: Double -> Double
1258 "realToFrac/Int->Double" realToFrac = int2Double -- See Note [realToFrac int-to-float]
1259 "realToFrac/Int->Float" realToFrac = int2Float -- ..ditto
1260 #-}
1261
1262 {-
1263 Note [realToFrac int-to-float]
1264 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1265 Don found that the RULES for realToFrac/Int->Double and simliarly
1266 Float made a huge difference to some stream-fusion programs. Here's
1267 an example
1268
1269 import Data.Array.Vector
1270
1271 n = 40000000
1272
1273 main = do
1274 let c = replicateU n (2::Double)
1275 a = mapU realToFrac (enumFromToU 0 (n-1) ) :: UArr Double
1276 print (sumU (zipWithU (*) c a))
1277
1278 Without the RULE we get this loop body:
1279
1280 case $wtoRational sc_sY4 of ww_aM7 { (# ww1_aM9, ww2_aMa #) ->
1281 case $wfromRat ww1_aM9 ww2_aMa of tpl_X1P { D# ipv_sW3 ->
1282 Main.$s$wfold
1283 (+# sc_sY4 1)
1284 (+# wild_X1i 1)
1285 (+## sc2_sY6 (*## 2.0 ipv_sW3))
1286
1287 And with the rule:
1288
1289 Main.$s$wfold
1290 (+# sc_sXT 1)
1291 (+# wild_X1h 1)
1292 (+## sc2_sXV (*## 2.0 (int2Double# sc_sXT)))
1293
1294 The running time of the program goes from 120 seconds to 0.198 seconds
1295 with the native backend, and 0.143 seconds with the C backend.
1296
1297 A few more details in #2251, and the patch message
1298 "Add RULES for realToFrac from Int".
1299 -}
1300
1301 -- Utils
1302
1303 showSignedFloat :: (RealFloat a)
1304 => (a -> ShowS) -- ^ a function that can show unsigned values
1305 -> Int -- ^ the precedence of the enclosing context
1306 -> a -- ^ the value to show
1307 -> ShowS
1308 showSignedFloat showPos p x
1309 | x < 0 || isNegativeZero x
1310 = showParen (p > 6) (showChar '-' . showPos (-x))
1311 | otherwise = showPos x
1312
1313 {-
1314 We need to prevent over/underflow of the exponent in encodeFloat when
1315 called from scaleFloat, hence we clamp the scaling parameter.
1316 We must have a large enough range to cover the maximum difference of
1317 exponents returned by decodeFloat.
1318 -}
1319 clamp :: Int -> Int -> Int
1320 clamp bd k = max (-bd) (min bd k)
1321
1322
1323 {-
1324 Note [Casting from integral to floating point types]
1325 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1326 To implement something like `reinterpret_cast` from C++ to go from a
1327 floating-point type to an integral type one might niavely think that the
1328 following should work:
1329
1330 cast :: Float -> Word32
1331 cast (F# f#) = W32# (unsafeCoerce# f#)
1332
1333 Unfortunately that is not the case, because all the `unsafeCoerce#` does is tell
1334 the compiler that the types have changed. When one does the above cast and
1335 tries to operate on the resulting `Word32` the code generator will generate code
1336 that performs an integer/word operation on a floating-point register, which
1337 results in a compile error.
1338
1339 The correct way of implementing `reinterpret_cast` to implement a primpop, but
1340 that requires a unique implementation for all supported archetectures. The next
1341 best solution is to write the value from the source register to memory and then
1342 read it from memory into the destination register and the best way to do that
1343 is using CMM.
1344 -}
1345
1346 -- | @'castWord32ToFloat' w@ does a bit-for-bit copy from an integral value
1347 -- to a floating-point value.
1348 --
1349 -- @since 4.10.0.0
1350
1351 {-# INLINE castWord32ToFloat #-}
1352 castWord32ToFloat :: Word32 -> Float
1353 castWord32ToFloat (W32# w#) = F# (stgWord32ToFloat w#)
1354
1355 foreign import prim "stg_word32ToFloatzh"
1356 stgWord32ToFloat :: Word# -> Float#
1357
1358
1359 -- | @'castFloatToWord32' f@ does a bit-for-bit copy from a floating-point value
1360 -- to an integral value.
1361 --
1362 -- @since 4.10.0.0
1363
1364 {-# INLINE castFloatToWord32 #-}
1365 castFloatToWord32 :: Float -> Word32
1366 castFloatToWord32 (F# f#) = W32# (stgFloatToWord32 f#)
1367
1368 foreign import prim "stg_floatToWord32zh"
1369 stgFloatToWord32 :: Float# -> Word#
1370
1371
1372
1373 -- | @'castWord64ToDouble' w@ does a bit-for-bit copy from an integral value
1374 -- to a floating-point value.
1375 --
1376 -- @since 4.10.0.0
1377
1378 {-# INLINE castWord64ToDouble #-}
1379 castWord64ToDouble :: Word64 -> Double
1380 castWord64ToDouble (W64# w) = D# (stgWord64ToDouble w)
1381
1382 foreign import prim "stg_word64ToDoublezh"
1383 #if WORD_SIZE_IN_BITS == 64
1384 stgWord64ToDouble :: Word# -> Double#
1385 #else
1386 stgWord64ToDouble :: Word64# -> Double#
1387 #endif
1388
1389
1390 -- | @'castFloatToWord32' f@ does a bit-for-bit copy from a floating-point value
1391 -- to an integral value.
1392 --
1393 -- @since 4.10.0.0
1394
1395 {-# INLINE castDoubleToWord64 #-}
1396 castDoubleToWord64 :: Double -> Word64
1397 castDoubleToWord64 (D# d#) = W64# (stgDoubleToWord64 d#)
1398
1399 foreign import prim "stg_doubleToWord64zh"
1400 #if WORD_SIZE_IN_BITS == 64
1401 stgDoubleToWord64 :: Double# -> Word#
1402 #else
1403 stgDoubleToWord64 :: Double# -> Word64#
1404 #endif