1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
-- !! fromRational woes
import Data.Ratio -- 1.3
main = putStr (
shows tinyFloat ( '\n'
: shows t_f ( '\n'
: shows hugeFloat ( '\n'
: shows h_f ( '\n'
: shows tinyDouble ( '\n'
: shows t_d ( '\n'
: shows hugeDouble ( '\n'
: shows h_d ( '\n'
: shows x_f ( '\n'
: shows x_d ( '\n'
: shows y_f ( '\n'
: shows y_d ( "\n"
)))))))))))))
where
t_f :: Float
t_d :: Double
h_f :: Float
h_d :: Double
x_f :: Float
x_d :: Double
y_f :: Float
y_d :: Double
t_f = fromRationalX (toRational tinyFloat)
t_d = fromRationalX (toRational tinyDouble)
h_f = fromRationalX (toRational hugeFloat)
h_d = fromRationalX (toRational hugeDouble)
x_f = fromRationalX (1.82173691287639817263897126389712638972163e-300 :: Rational)
x_d = fromRationalX (1.82173691287639817263897126389712638972163e-300 :: Rational)
y_f = 1.82173691287639817263897126389712638972163e-300
y_d = 1.82173691287639817263897126389712638972163e-300
fromRationalX :: (RealFloat a) => Rational -> a
fromRationalX r =
let
h = ceiling (huge `asTypeOf` x)
b = toInteger (floatRadix x)
x = fromRat 0 r
fromRat e0 r' =
let d = denominator r'
n = numerator r'
in if d > h then
let e = integerLogBase b (d `div` h) + 1
in fromRat (e0-e) (n % (d `div` (b^e)))
else if abs n > h then
let e = integerLogBase b (abs n `div` h) + 1
in fromRat (e0+e) ((n `div` (b^e)) % d)
else
scaleFloat e0 (rationalToRealFloat {-fromRational-} r')
in x
{-
fromRationalX r =
rationalToRealFloat r
{- Hmmm...
let
h = ceiling (huge `asTypeOf` x)
b = toInteger (floatRadix x)
x = fromRat 0 r
fromRat e0 r' =
{--} trace (shows e0 ('/' : shows r' ('/' : shows h "\n"))) (
let d = denominator r'
n = numerator r'
in if d > h then
let e = integerLogBase b (d `div` h) + 1
in fromRat (e0-e) (n % (d `div` (b^e)))
else if abs n > h then
let e = integerLogBase b (abs n `div` h) + 1
in fromRat (e0+e) ((n `div` (b^e)) % d)
else
scaleFloat e0 (rationalToRealFloat r')
-- now that we know things are in-bounds,
-- we use the "old" Prelude code.
{--} )
in x
-}
-}
-- Compute the discrete log of i in base b.
-- Simplest way would be just divide i by b until it's smaller then b, but that would
-- be very slow! We are just slightly more clever.
integerLogBase :: Integer -> Integer -> Int
integerLogBase b i =
if i < b then
0
else
-- Try squaring the base first to cut down the number of divisions.
let l = 2 * integerLogBase (b*b) i
doDiv :: Integer -> Int -> Int
doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
in doDiv (i `div` (b^l)) l
------------
-- Compute smallest and largest floating point values.
tiny :: (RealFloat a) => a
tiny =
let (l, _) = floatRange x
x = encodeFloat 1 (l-1)
in x
huge :: (RealFloat a) => a
huge =
let (_, u) = floatRange x
d = floatDigits x
x = encodeFloat (floatRadix x ^ d - 1) (u - d)
in x
tinyDouble = tiny :: Double
tinyFloat = tiny :: Float
hugeDouble = huge :: Double
hugeFloat = huge :: Float
{-
[In response to a request by simonpj, Joe Fasel writes:]
A quite reasonable request! This code was added to the Prelude just
before the 1.2 release, when Lennart, working with an early version
of hbi, noticed that (read . show) was not the identity for
floating-point numbers. (There was a one-bit error about half the time.)
The original version of the conversion function was in fact simply
a floating-point divide, as you suggest above. The new version is,
I grant you, somewhat denser.
How's this?
--Joe
-}
rationalToRealFloat :: (RealFloat a) => Rational -> a
rationalToRealFloat x = x'
where x' = f e
-- If the exponent of the nearest floating-point number to x
-- is e, then the significand is the integer nearest xb^(-e),
-- where b is the floating-point radix. We start with a good
-- guess for e, and if it is correct, the exponent of the
-- floating-point number we construct will again be e. If
-- not, one more iteration is needed.
f e = if e' == e then y else f e'
where y = encodeFloat (round (x * (1%b)^^e)) e
(_,e') = decodeFloat y
b = floatRadix x'
-- We obtain a trial exponent by doing a floating-point
-- division of x's numerator by its denominator. The
-- result of this division may not itself be the ultimate
-- result, because of an accumulation of three rounding
-- errors.
(s,e) = decodeFloat (fromInteger (numerator x) `asTypeOf` x'
/ fromInteger (denominator x))
|
