author  wenzelm 
Sun, 01 Mar 2009 23:36:12 +0100  
changeset 30190  479806475f3c 
parent 29804  e15b74577368 
child 37391  476270a6c2dc 
permissions  rwrr 
29804
e15b74577368
Added new Float theory and moved old Library/Float.thy to ComputeFloat
hoelzl
parents:
28952
diff
changeset

1 
(* Title: HOL/Tools/float_arith.ML 
22964  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 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

11 
val approx_dec_by_bin: int > Float.float > Float.float * Float.float 
22964  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; 

23569  87 
fun exp5 x = Integer.pow x 5; 
88 
fun exp10 x = Integer.pow x 10; 

25300  89 
fun exp2 x = Integer.pow x 2; 
22964  90 

91 
fun find_most_significant q r = 

92 
let 

93 
fun int2real i = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

94 
case (Real.fromString o string_of_int) i of 
22964  95 
SOME r => r 
96 
 NONE => raise (Floating_point "int2real") 

97 
fun subtract (q, r) (q', r') = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

98 
if r <= r' then 
23297  99 
(q  q' * exp10 (r'  r), r) 
22964  100 
else 
23297  101 
(q * exp10 (r  r')  q', r') 
22964  102 
fun bin2dec d = 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

103 
if 0 <= d then 
25300  104 
(exp2 d, 0) 
22964  105 
else 
23297  106 
(exp5 (~ d), d) 
22964  107 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

108 
val L = Real.floor (int2real (IntInf.log2 q) + int2real r * ln2_10) 
23297  109 
val L1 = L + 1 
22964  110 

111 
val (q1, r1) = subtract (q, r) (bin2dec L1) 

112 
in 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

113 
if 0 <= q1 then 
22964  114 
let 
23297  115 
val (q2, r2) = subtract (q, r) (bin2dec (L1 + 1)) 
22964  116 
in 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

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 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

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 

25300  136 
 addseq acc d' (d::ds) = addseq (acc + exp2 (d  d')) d' ds 
22964  137 

23297  138 
fun seq2bin [] = (0, 0) 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

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 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

150 
if precision < d0  d then 
22964  151 
let 
23297  152 
val d' = d0  precision 
22964  153 
val x1 = seq2bin (d_seq) 
25300  154 
val x2 = (fst x1 * exp2 (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) = 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

163 
if q = 0 then 
23297  164 
((0, 0), (0, 0)) 
22964  165 
else 
166 
let 

167 
val (d, (q', r')) = find_most_significant q r 

168 
in 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

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 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

176 
(x1, seq2bin [d + 1]) 
22964  177 
end 
178 
else 

179 
approx [d] d precision (q', r') 

180 
end 

181 
in 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

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 

24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

194 
fun str2int s = the_default 0 (Int.fromString s) 
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

195 
fun signint p x = if p then x else ~ x 
22964  196 

197 
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

198 
val s = 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 
24630
351a308ab58d
simplified type int (eliminated IntInf.int, integer);
wenzelm
parents:
24584
diff
changeset

203 
approx_dec_by_bin 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; 