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