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