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