|
1 (* Title: HOL/Real/float_arith.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 |
|
12 val approx_dec_by_bin: int -> Float.float -> Float.float * Float.float |
|
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; |
|
88 fun exp5 x = Integer.pow x 5; |
|
89 fun exp10 x = Integer.pow x 10; |
|
90 fun exp2 x = Integer.pow x 2; |
|
91 |
|
92 fun find_most_significant q r = |
|
93 let |
|
94 fun int2real i = |
|
95 case (Real.fromString o string_of_int) i of |
|
96 SOME r => r |
|
97 | NONE => raise (Floating_point "int2real") |
|
98 fun subtract (q, r) (q', r') = |
|
99 if r <= r' then |
|
100 (q - q' * exp10 (r' - r), r) |
|
101 else |
|
102 (q * exp10 (r - r') - q', r') |
|
103 fun bin2dec d = |
|
104 if 0 <= d then |
|
105 (exp2 d, 0) |
|
106 else |
|
107 (exp5 (~ d), d) |
|
108 |
|
109 val L = Real.floor (int2real (IntInf.log2 q) + int2real r * ln2_10) |
|
110 val L1 = L + 1 |
|
111 |
|
112 val (q1, r1) = subtract (q, r) (bin2dec L1) |
|
113 in |
|
114 if 0 <= q1 then |
|
115 let |
|
116 val (q2, r2) = subtract (q, r) (bin2dec (L1 + 1)) |
|
117 in |
|
118 if 0 <= q2 then |
|
119 raise (Floating_point "find_most_significant") |
|
120 else |
|
121 (L1, (q1, r1)) |
|
122 end |
|
123 else |
|
124 let |
|
125 val (q0, r0) = subtract (q, r) (bin2dec L) |
|
126 in |
|
127 if 0 <= q0 then |
|
128 (L, (q0, r0)) |
|
129 else |
|
130 raise (Floating_point "find_most_significant") |
|
131 end |
|
132 end |
|
133 |
|
134 fun approx_dec_by_bin n (q,r) = |
|
135 let |
|
136 fun addseq acc d' [] = acc |
|
137 | addseq acc d' (d::ds) = addseq (acc + exp2 (d - d')) d' ds |
|
138 |
|
139 fun seq2bin [] = (0, 0) |
|
140 | seq2bin (d::ds) = (addseq 0 d ds + 1, d) |
|
141 |
|
142 fun approx d_seq d0 precision (q,r) = |
|
143 if q = 0 then |
|
144 let val x = seq2bin d_seq in |
|
145 (x, x) |
|
146 end |
|
147 else |
|
148 let |
|
149 val (d, (q', r')) = find_most_significant q r |
|
150 in |
|
151 if precision < d0 - d then |
|
152 let |
|
153 val d' = d0 - precision |
|
154 val x1 = seq2bin (d_seq) |
|
155 val x2 = (fst x1 * exp2 (snd x1 - d') + 1, d') (* = seq2bin (d'::d_seq) *) |
|
156 in |
|
157 (x1, x2) |
|
158 end |
|
159 else |
|
160 approx (d::d_seq) d0 precision (q', r') |
|
161 end |
|
162 |
|
163 fun approx_start precision (q, r) = |
|
164 if q = 0 then |
|
165 ((0, 0), (0, 0)) |
|
166 else |
|
167 let |
|
168 val (d, (q', r')) = find_most_significant q r |
|
169 in |
|
170 if precision <= 0 then |
|
171 let |
|
172 val x1 = seq2bin [d] |
|
173 in |
|
174 if q' = 0 then |
|
175 (x1, x1) |
|
176 else |
|
177 (x1, seq2bin [d + 1]) |
|
178 end |
|
179 else |
|
180 approx [d] d precision (q', r') |
|
181 end |
|
182 in |
|
183 if 0 <= q then |
|
184 approx_start n (q,r) |
|
185 else |
|
186 let |
|
187 val ((a1,b1), (a2, b2)) = approx_start n (~ q, r) |
|
188 in |
|
189 ((~ a2, b2), (~ a1, b1)) |
|
190 end |
|
191 end |
|
192 |
|
193 fun approx_decstr_by_bin n decstr = |
|
194 let |
|
195 fun str2int s = the_default 0 (Int.fromString s) |
|
196 fun signint p x = if p then x else ~ x |
|
197 |
|
198 val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr |
|
199 val s = size d2 |
|
200 |
|
201 val q = signint p (str2int d1 * exp10 s + str2int d2) |
|
202 val r = signint ep (str2int e) - s |
|
203 in |
|
204 approx_dec_by_bin n (q,r) |
|
205 end |
|
206 |
|
207 fun mk_float (a, b) = @{term "float"} $ |
|
208 HOLogic.mk_prod (pairself (HOLogic.mk_number HOLogic.intT) (a, b)); |
|
209 |
|
210 fun dest_float (Const ("Float.float", _) $ (Const ("Pair", _) $ a $ b)) = |
|
211 pairself (snd o HOLogic.dest_number) (a, b) |
|
212 | dest_float t = ((snd o HOLogic.dest_number) t, 0); |
|
213 |
|
214 fun approx_float prec f value = |
|
215 let |
|
216 val interval = approx_decstr_by_bin prec value |
|
217 val (flower, fupper) = f interval |
|
218 in |
|
219 (mk_float flower, mk_float fupper) |
|
220 end; |
|
221 |
|
222 end; |