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