22964
|
1 |
(* Title: HOL/Real/Float.ML
|
|
2 |
ID: $Id$
|
|
3 |
Author: Steven Obua
|
|
4 |
*)
|
|
5 |
|
|
6 |
signature FLOAT_ARITH =
|
|
7 |
sig
|
|
8 |
exception Destruct_floatstr of string
|
|
9 |
val destruct_floatstr: (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string
|
|
10 |
|
|
11 |
exception Floating_point of string
|
23261
|
12 |
val approx_dec_by_bin: integer -> Float.float -> Float.float * Float.float
|
22964
|
13 |
val approx_decstr_by_bin: int -> string -> Float.float * Float.float
|
|
14 |
|
|
15 |
val mk_float: Float.float -> term
|
|
16 |
val dest_float: term -> Float.float
|
|
17 |
|
|
18 |
val approx_float: int -> (Float.float * Float.float -> Float.float * Float.float)
|
|
19 |
-> string -> term * term
|
|
20 |
end;
|
|
21 |
|
|
22 |
structure FloatArith : FLOAT_ARITH =
|
|
23 |
struct
|
|
24 |
|
|
25 |
exception Destruct_floatstr of string;
|
|
26 |
|
|
27 |
fun destruct_floatstr isDigit isExp number =
|
|
28 |
let
|
|
29 |
val numlist = filter (not o Char.isSpace) (String.explode number)
|
|
30 |
|
|
31 |
fun countsigns ((#"+")::cs) = countsigns cs
|
|
32 |
| countsigns ((#"-")::cs) =
|
|
33 |
let
|
|
34 |
val (positive, rest) = countsigns cs
|
|
35 |
in
|
|
36 |
(not positive, rest)
|
|
37 |
end
|
|
38 |
| countsigns cs = (true, cs)
|
|
39 |
|
|
40 |
fun readdigits [] = ([], [])
|
|
41 |
| readdigits (q as c::cs) =
|
|
42 |
if (isDigit c) then
|
|
43 |
let
|
|
44 |
val (digits, rest) = readdigits cs
|
|
45 |
in
|
|
46 |
(c::digits, rest)
|
|
47 |
end
|
|
48 |
else
|
|
49 |
([], q)
|
|
50 |
|
|
51 |
fun readfromexp_helper cs =
|
|
52 |
let
|
|
53 |
val (positive, rest) = countsigns cs
|
|
54 |
val (digits, rest') = readdigits rest
|
|
55 |
in
|
|
56 |
case rest' of
|
|
57 |
[] => (positive, digits)
|
|
58 |
| _ => raise (Destruct_floatstr number)
|
|
59 |
end
|
|
60 |
|
|
61 |
fun readfromexp [] = (true, [])
|
|
62 |
| readfromexp (c::cs) =
|
|
63 |
if isExp c then
|
|
64 |
readfromexp_helper cs
|
|
65 |
else
|
|
66 |
raise (Destruct_floatstr number)
|
|
67 |
|
|
68 |
fun readfromdot [] = ([], readfromexp [])
|
|
69 |
| readfromdot ((#".")::cs) =
|
|
70 |
let
|
|
71 |
val (digits, rest) = readdigits cs
|
|
72 |
val exp = readfromexp rest
|
|
73 |
in
|
|
74 |
(digits, exp)
|
|
75 |
end
|
|
76 |
| readfromdot cs = readfromdot ((#".")::cs)
|
|
77 |
|
|
78 |
val (positive, numlist) = countsigns numlist
|
|
79 |
val (digits1, numlist) = readdigits numlist
|
|
80 |
val (digits2, exp) = readfromdot numlist
|
|
81 |
in
|
|
82 |
(positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp))
|
|
83 |
end
|
|
84 |
|
|
85 |
exception Floating_point of string;
|
|
86 |
|
|
87 |
val ln2_10 = Math.ln 10.0 / Math.ln 2.0;
|
23297
|
88 |
val exp5 = Integer.pow 5;
|
|
89 |
val exp10 = Integer.pow 10;
|
22964
|
90 |
|
|
91 |
fun find_most_significant q r =
|
|
92 |
let
|
|
93 |
fun int2real i =
|
23261
|
94 |
case (Real.fromString o Integer.string_of_int) i of
|
22964
|
95 |
SOME r => r
|
|
96 |
| NONE => raise (Floating_point "int2real")
|
|
97 |
fun subtract (q, r) (q', r') =
|
23297
|
98 |
if r <=% r' then
|
|
99 |
(q - q' * exp10 (r' - r), r)
|
22964
|
100 |
else
|
23297
|
101 |
(q * exp10 (r - r') - q', r')
|
22964
|
102 |
fun bin2dec d =
|
23297
|
103 |
if 0 <=% d then
|
|
104 |
(Integer.exp d, 0)
|
22964
|
105 |
else
|
23297
|
106 |
(exp5 (~ d), d)
|
22964
|
107 |
|
23261
|
108 |
val L = Integer.int (Real.floor (int2real (Integer.log q) + int2real r * ln2_10))
|
23297
|
109 |
val L1 = L + 1
|
22964
|
110 |
|
|
111 |
val (q1, r1) = subtract (q, r) (bin2dec L1)
|
|
112 |
in
|
23297
|
113 |
if 0 <=% q1 then
|
22964
|
114 |
let
|
23297
|
115 |
val (q2, r2) = subtract (q, r) (bin2dec (L1 + 1))
|
22964
|
116 |
in
|
23297
|
117 |
if 0 <=% q2 then
|
22964
|
118 |
raise (Floating_point "find_most_significant")
|
|
119 |
else
|
|
120 |
(L1, (q1, r1))
|
|
121 |
end
|
|
122 |
else
|
|
123 |
let
|
|
124 |
val (q0, r0) = subtract (q, r) (bin2dec L)
|
|
125 |
in
|
23297
|
126 |
if 0 <=% q0 then
|
22964
|
127 |
(L, (q0, r0))
|
|
128 |
else
|
|
129 |
raise (Floating_point "find_most_significant")
|
|
130 |
end
|
|
131 |
end
|
|
132 |
|
|
133 |
fun approx_dec_by_bin n (q,r) =
|
|
134 |
let
|
|
135 |
fun addseq acc d' [] = acc
|
23297
|
136 |
| addseq acc d' (d::ds) = addseq (acc +% Integer.exp (d - d')) d' ds
|
22964
|
137 |
|
23297
|
138 |
fun seq2bin [] = (0, 0)
|
|
139 |
| seq2bin (d::ds) = (addseq 0 d ds +% 1, d)
|
22964
|
140 |
|
|
141 |
fun approx d_seq d0 precision (q,r) =
|
23297
|
142 |
if q = 0 then
|
22964
|
143 |
let val x = seq2bin d_seq in
|
|
144 |
(x, x)
|
|
145 |
end
|
|
146 |
else
|
|
147 |
let
|
|
148 |
val (d, (q', r')) = find_most_significant q r
|
|
149 |
in
|
23297
|
150 |
if precision <% d0 - d then
|
22964
|
151 |
let
|
23297
|
152 |
val d' = d0 - precision
|
22964
|
153 |
val x1 = seq2bin (d_seq)
|
23297
|
154 |
val x2 = (fst x1 * Integer.exp (snd x1 - d') + 1, d') (* = seq2bin (d'::d_seq) *)
|
22964
|
155 |
in
|
|
156 |
(x1, x2)
|
|
157 |
end
|
|
158 |
else
|
|
159 |
approx (d::d_seq) d0 precision (q', r')
|
|
160 |
end
|
|
161 |
|
|
162 |
fun approx_start precision (q, r) =
|
23297
|
163 |
if q =% 0 then
|
|
164 |
((0, 0), (0, 0))
|
22964
|
165 |
else
|
|
166 |
let
|
|
167 |
val (d, (q', r')) = find_most_significant q r
|
|
168 |
in
|
23297
|
169 |
if precision <=% 0 then
|
22964
|
170 |
let
|
|
171 |
val x1 = seq2bin [d]
|
|
172 |
in
|
23297
|
173 |
if q' = 0 then
|
22964
|
174 |
(x1, x1)
|
|
175 |
else
|
23297
|
176 |
(x1, seq2bin [d +% 1])
|
22964
|
177 |
end
|
|
178 |
else
|
|
179 |
approx [d] d precision (q', r')
|
|
180 |
end
|
|
181 |
in
|
23297
|
182 |
if 0 <=% q then
|
22964
|
183 |
approx_start n (q,r)
|
|
184 |
else
|
|
185 |
let
|
23297
|
186 |
val ((a1,b1), (a2, b2)) = approx_start n (~ q, r)
|
22964
|
187 |
in
|
23297
|
188 |
((~ a2, b2), (~ a1, b1))
|
22964
|
189 |
end
|
|
190 |
end
|
|
191 |
|
|
192 |
fun approx_decstr_by_bin n decstr =
|
|
193 |
let
|
23297
|
194 |
fun str2int s = the_default 0 (Integer.int_of_string s);
|
23261
|
195 |
fun signint p x = if p then x else Integer.neg x
|
22964
|
196 |
|
|
197 |
val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr
|
23261
|
198 |
val s = Integer.int (size d2)
|
22964
|
199 |
|
23297
|
200 |
val q = signint p (str2int d1 * exp10 s + str2int d2)
|
|
201 |
val r = signint ep (str2int e) - s
|
22964
|
202 |
in
|
23261
|
203 |
approx_dec_by_bin (Integer.int n) (q,r)
|
22964
|
204 |
end
|
|
205 |
|
|
206 |
fun mk_float (a, b) = @{term "float"} $
|
|
207 |
HOLogic.mk_prod (pairself (HOLogic.mk_number HOLogic.intT) (a, b));
|
|
208 |
|
|
209 |
fun dest_float (Const ("Float.float", _) $ (Const ("Pair", _) $ a $ b)) =
|
|
210 |
pairself (snd o HOLogic.dest_number) (a, b)
|
23297
|
211 |
| dest_float t = ((snd o HOLogic.dest_number) t, 0);
|
22964
|
212 |
|
|
213 |
fun approx_float prec f value =
|
|
214 |
let
|
|
215 |
val interval = approx_decstr_by_bin prec value
|
|
216 |
val (flower, fupper) = f interval
|
|
217 |
in
|
|
218 |
(mk_float flower, mk_float fupper)
|
|
219 |
end;
|
|
220 |
|
|
221 |
end;
|