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