author | haftmann |
Wed, 03 Dec 2008 15:58:44 +0100 | |
changeset 28952 | 15a4b2cf8c34 |
parent 25300 | src/HOL/Real/float_arith.ML@bc3a1c964704 |
child 29804 | e15b74577368 |
permissions | -rw-r--r-- |
24584 | 1 |
(* Title: HOL/Real/float_arith.ML |
2 |
ID: $Id$ |
|
22964 | 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 |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
12 |
val approx_dec_by_bin: int -> 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; |
|
23569 | 88 |
fun exp5 x = Integer.pow x 5; |
89 |
fun exp10 x = Integer.pow x 10; |
|
25300 | 90 |
fun exp2 x = Integer.pow x 2; |
22964 | 91 |
|
92 |
fun find_most_significant q r = |
|
93 |
let |
|
94 |
fun int2real i = |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
95 |
case (Real.fromString o string_of_int) i of |
22964 | 96 |
SOME r => r |
97 |
| NONE => raise (Floating_point "int2real") |
|
98 |
fun subtract (q, r) (q', r') = |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
99 |
if r <= r' then |
23297 | 100 |
(q - q' * exp10 (r' - r), r) |
22964 | 101 |
else |
23297 | 102 |
(q * exp10 (r - r') - q', r') |
22964 | 103 |
fun bin2dec d = |
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
104 |
if 0 <= d then |
25300 | 105 |
(exp2 d, 0) |
22964 | 106 |
else |
23297 | 107 |
(exp5 (~ d), d) |
22964 | 108 |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
109 |
val L = Real.floor (int2real (IntInf.log2 q) + int2real r * ln2_10) |
23297 | 110 |
val L1 = L + 1 |
22964 | 111 |
|
112 |
val (q1, r1) = subtract (q, r) (bin2dec L1) |
|
113 |
in |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
114 |
if 0 <= q1 then |
22964 | 115 |
let |
23297 | 116 |
val (q2, r2) = subtract (q, r) (bin2dec (L1 + 1)) |
22964 | 117 |
in |
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
118 |
if 0 <= q2 then |
22964 | 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 |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
127 |
if 0 <= q0 then |
22964 | 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 |
|
25300 | 137 |
| addseq acc d' (d::ds) = addseq (acc + exp2 (d - d')) d' ds |
22964 | 138 |
|
23297 | 139 |
fun seq2bin [] = (0, 0) |
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
140 |
| seq2bin (d::ds) = (addseq 0 d ds + 1, d) |
22964 | 141 |
|
142 |
fun approx d_seq d0 precision (q,r) = |
|
23297 | 143 |
if q = 0 then |
22964 | 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 |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
151 |
if precision < d0 - d then |
22964 | 152 |
let |
23297 | 153 |
val d' = d0 - precision |
22964 | 154 |
val x1 = seq2bin (d_seq) |
25300 | 155 |
val x2 = (fst x1 * exp2 (snd x1 - d') + 1, d') (* = seq2bin (d'::d_seq) *) |
22964 | 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) = |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
164 |
if q = 0 then |
23297 | 165 |
((0, 0), (0, 0)) |
22964 | 166 |
else |
167 |
let |
|
168 |
val (d, (q', r')) = find_most_significant q r |
|
169 |
in |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
170 |
if precision <= 0 then |
22964 | 171 |
let |
172 |
val x1 = seq2bin [d] |
|
173 |
in |
|
23297 | 174 |
if q' = 0 then |
22964 | 175 |
(x1, x1) |
176 |
else |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
177 |
(x1, seq2bin [d + 1]) |
22964 | 178 |
end |
179 |
else |
|
180 |
approx [d] d precision (q', r') |
|
181 |
end |
|
182 |
in |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
183 |
if 0 <= q then |
22964 | 184 |
approx_start n (q,r) |
185 |
else |
|
186 |
let |
|
23297 | 187 |
val ((a1,b1), (a2, b2)) = approx_start n (~ q, r) |
22964 | 188 |
in |
23297 | 189 |
((~ a2, b2), (~ a1, b1)) |
22964 | 190 |
end |
191 |
end |
|
192 |
||
193 |
fun approx_decstr_by_bin n decstr = |
|
194 |
let |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
195 |
fun str2int s = the_default 0 (Int.fromString s) |
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
196 |
fun signint p x = if p then x else ~ x |
22964 | 197 |
|
198 |
val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr |
|
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
199 |
val s = size d2 |
22964 | 200 |
|
23297 | 201 |
val q = signint p (str2int d1 * exp10 s + str2int d2) |
202 |
val r = signint ep (str2int e) - s |
|
22964 | 203 |
in |
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset
|
204 |
approx_dec_by_bin n (q,r) |
22964 | 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) |
|
23297 | 212 |
| dest_float t = ((snd o HOLogic.dest_number) t, 0); |
22964 | 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; |