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