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