~ubuntu-branches/ubuntu/precise/ghc/precise

« back to all changes in this revision

Viewing changes to libraries/base/GHC/Float.lhs

  • Committer: Bazaar Package Importer
  • Author(s): Joachim Breitner
  • Date: 2011-01-17 12:49:24 UTC
  • Revision ID: james.westby@ubuntu.com-20110117124924-do1pym1jlf5o636m
Tags: upstream-7.0.1
ImportĀ upstreamĀ versionĀ 7.0.1

Show diffs side-by-side

added added

removed removed

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