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