author  hoelzl 
Wed, 11 Mar 2009 10:58:18 +0100  
changeset 30439  57c68b3af2ea 
parent 30413  c41afa5607be 
child 30443  873fa77be5f0 
permissions  rwrr 
30439  1 
(* Title: HOL/Decision_Procs/Approximation.thy 
30122  2 
Author: Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009 
3 
*) 

4 

29805  5 
header {* Prove unequations about real numbers by computation *} 
30122  6 

29805  7 
theory Approximation 
29823
0ab754d13ccd
session Reflecion renamed to Decision_Procs, moved Dense_Linear_Order there
haftmann
parents:
29805
diff
changeset

8 
imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat 
29805  9 
begin 
10 

11 
section "Horner Scheme" 

12 

13 
subsection {* Define auxiliary helper @{text horner} function *} 

14 

15 
fun horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where 

16 
"horner F G 0 i k x = 0"  

17 
"horner F G (Suc n) i k x = 1 / real k  x * horner F G n (F i) (G i k) x" 

18 

19 
lemma horner_schema': fixes x :: real and a :: "nat \<Rightarrow> real" 

20 
shows "a 0  x * (\<Sum> i=0..<n. (1)^i * a (Suc i) * x^i) = (\<Sum> i=0..<Suc n. (1)^i * a i * x^i)" 

21 
proof  

22 
have shift_pow: "\<And>i.  (x * ((1)^i * a (Suc i) * x ^ i)) = (1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto 

23 
show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc] 

24 
setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (1)^n *a n * x^n"] by auto 

25 
qed 

26 

27 
lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat" 

28 
assumes f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)" 

29 
shows "horner F G n ((F^j') s) (f j') x = (\<Sum> j = 0..< n. 1^j * (1 / real (f (j' + j))) * x^j)" 

30 
proof (induct n arbitrary: i k j') 

31 
case (Suc n) 

32 

33 
show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc] 

34 
using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto 

35 
qed auto 

36 

37 
lemma horner_bounds': 

38 
assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)" 

39 
and lb_0: "\<And> i k x. lb 0 i k x = 0" 

40 
and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k)  x * (ub n (F i) (G i k) x)" 

41 
and ub_0: "\<And> i k x. ub 0 i k x = 0" 

42 
and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k)  x * (lb n (F i) (G i k) x)" 

43 
shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> horner F G n ((F^j') s) (f j') (Ifloat x) \<and> 

44 
horner F G n ((F^j') s) (f j') (Ifloat x) \<le> Ifloat (ub n ((F^j') s) (f j') x)" 

45 
(is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'") 

46 
proof (induct n arbitrary: j') 

47 
case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto 

48 
next 

49 
case (Suc n) 

50 
have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps Ifloat_sub diff_def 

51 
proof (rule add_mono) 

52 
show "Ifloat (lapprox_rat prec 1 (int (f j'))) \<le> 1 / real (f j')" using lapprox_rat[of prec 1 "int (f j')"] by auto 

53 
from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> Ifloat x` 

54 
show " Ifloat (x * ub n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) x) \<le>  (Ifloat x * horner F G n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) (Ifloat x))" 

55 
unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono) 

56 
qed 

57 
moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps Ifloat_sub diff_def 

58 
proof (rule add_mono) 

59 
show "1 / real (f j') \<le> Ifloat (rapprox_rat prec 1 (int (f j')))" using rapprox_rat[of 1 "int (f j')" prec] by auto 

60 
from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> Ifloat x` 

61 
show " (Ifloat x * horner F G n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) (Ifloat x)) \<le> 

62 
 Ifloat (x * lb n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) x)" 

63 
unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono) 

64 
qed 

65 
ultimately show ?case by blast 

66 
qed 

67 

68 
subsection "Theorems for floating point functions implementing the horner scheme" 

69 

70 
text {* 

71 

72 
Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are 

73 
all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}. 

74 

75 
*} 

76 

77 
lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" 

78 
assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)" 

79 
and lb_0: "\<And> i k x. lb 0 i k x = 0" 

80 
and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k)  x * (ub n (F i) (G i k) x)" 

81 
and ub_0: "\<And> i k x. ub 0 i k x = 0" 

82 
and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k)  x * (lb n (F i) (G i k) x)" 

83 
shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> (\<Sum>j=0..<n. 1^j * (1 / real (f (j' + j))) * (Ifloat x)^j)" (is "?lb") and 

84 
"(\<Sum>j=0..<n. 1^j * (1 / real (f (j' + j))) * (Ifloat x)^j) \<le> Ifloat (ub n ((F^j') s) (f j') x)" (is "?ub") 

85 
proof  

86 
have "?lb \<and> ?ub" 

87 
using horner_bounds'[where lb=lb, OF `0 \<le> Ifloat x` f_Suc lb_0 lb_Suc ub_0 ub_Suc] 

88 
unfolding horner_schema[where f=f, OF f_Suc] . 

89 
thus "?lb" and "?ub" by auto 

90 
qed 

91 

92 
lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" 

93 
assumes "Ifloat x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)" 

94 
and lb_0: "\<And> i k x. lb 0 i k x = 0" 

95 
and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) + x * (ub n (F i) (G i k) x)" 

96 
and ub_0: "\<And> i k x. ub 0 i k x = 0" 

97 
and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) + x * (lb n (F i) (G i k) x)" 

98 
shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j)" (is "?lb") and 

99 
"(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) \<le> Ifloat (ub n ((F^j') s) (f j') x)" (is "?ub") 

100 
proof  

101 
{ fix x y z :: float have "x  y * z = x +  y * z" 

102 
by (cases x, cases y, cases z, simp add: plus_float.simps minus_float.simps uminus_float.simps times_float.simps algebra_simps) 

103 
} note diff_mult_minus = this 

104 

105 
{ fix x :: float have " ( x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this 

106 

107 
have move_minus: "Ifloat (x) = 1 * Ifloat x" by auto 

108 

109 
have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) = 

110 
(\<Sum>j = 0..<n. 1 ^ j * (1 / real (f (j' + j))) * Ifloat ( x) ^ j)" 

111 
proof (rule setsum_cong, simp) 

112 
fix j assume "j \<in> {0 ..< n}" 

113 
show "1 / real (f (j' + j)) * Ifloat x ^ j = 1 ^ j * (1 / real (f (j' + j))) * Ifloat ( x) ^ j" 

114 
unfolding move_minus power_mult_distrib real_mult_assoc[symmetric] 

115 
unfolding real_mult_commute unfolding real_mult_assoc[of "1^j", symmetric] power_mult_distrib[symmetric] 

116 
by auto 

117 
qed 

118 

119 
have "0 \<le> Ifloat (x)" using assms by auto 

120 
from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec 

121 
and lb="\<lambda> n i k x. lb n i k (x)" and ub="\<lambda> n i k x. ub n i k (x)", unfolded lb_Suc ub_Suc diff_mult_minus, 

122 
OF this f_Suc lb_0 refl ub_0 refl] 

123 
show "?lb" and "?ub" unfolding minus_minus sum_eq 

124 
by auto 

125 
qed 

126 

127 
subsection {* Selectors for next even or odd number *} 

128 

129 
text {* 

130 

131 
The horner scheme computes alternating series. To get the upper and lower bounds we need to 

132 
guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}. 

133 

134 
*} 

135 

136 
definition get_odd :: "nat \<Rightarrow> nat" where 

137 
"get_odd n = (if odd n then n else (Suc n))" 

138 

139 
definition get_even :: "nat \<Rightarrow> nat" where 

140 
"get_even n = (if even n then n else (Suc n))" 

141 

142 
lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto) 

143 
lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto) 

144 
lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)" 

145 
proof (cases "odd n") 

146 
case True hence "0 < n" by (rule odd_pos) 

147 
from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto 

148 
thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast 

149 
next 

150 
case False hence "odd (Suc n)" by auto 

151 
thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast 

152 
qed 

153 

154 
lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] . 

155 
lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto 

156 

157 
section "Power function" 

158 

159 
definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where 

160 
"float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n) 

161 
else if u < 0 then (u ^ n, l ^ n) 

162 
else (0, (max (l) u) ^ n))" 

163 

164 
lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {Ifloat l .. Ifloat u}" 

165 
shows "x^n \<in> {Ifloat l1..Ifloat u1}" 

166 
proof (cases "even n") 

167 
case True 

168 
show ?thesis 

169 
proof (cases "0 < l") 

170 
case True hence "odd n \<or> 0 < l" and "0 \<le> Ifloat l" unfolding less_float_def by auto 

171 
have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto 

172 
have "Ifloat l^n \<le> x^n" and "x^n \<le> Ifloat u^n " using `0 \<le> Ifloat l` and assms unfolding atLeastAtMost_iff using power_mono[of "Ifloat l" x] power_mono[of x "Ifloat u"] by auto 

173 
thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto 

174 
next 

175 
case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto 

176 
show ?thesis 

177 
proof (cases "u < 0") 

178 
case True hence "0 \<le>  Ifloat u" and " Ifloat u \<le>  x" and "0 \<le>  x" and "x \<le>  Ifloat l" using assms unfolding less_float_def by auto 

179 
hence "Ifloat u^n \<le> x^n" and "x^n \<le> Ifloat l^n" using power_mono[of "x" "Ifloat l" n] power_mono[of "Ifloat u" "x" n] 

180 
unfolding power_minus_even[OF `even n`] by auto 

181 
moreover have u1: "u1 = l ^ n" and l1: "l1 = u ^ n" using assms unfolding float_power_bnds_def if_not_P[OF P] if_P[OF True] by auto 

182 
ultimately show ?thesis using float_power by auto 

183 
next 

184 
case False 

185 
have "\<bar>x\<bar> \<le> Ifloat (max (l) u)" 

186 
proof (cases "l \<le> u") 

187 
case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto 

188 
next 

189 
case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto 

190 
qed 

191 
hence x_abs: "\<bar>x\<bar> \<le> \<bar>Ifloat (max (l) u)\<bar>" by auto 

192 
have u1: "u1 = (max (l) u) ^ n" and l1: "l1 = 0" using assms unfolding float_power_bnds_def if_not_P[OF P] if_not_P[OF False] by auto 

193 
show ?thesis unfolding atLeastAtMost_iff l1 u1 float_power using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto 

194 
qed 

195 
qed 

196 
next 

197 
case False hence "odd n \<or> 0 < l" by auto 

198 
have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto 

199 
have "Ifloat l^n \<le> x^n" and "x^n \<le> Ifloat u^n " using assms unfolding atLeastAtMost_iff using power_mono_odd[OF False] by auto 

200 
thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto 

201 
qed 

202 

203 
lemma bnds_power: "\<forall> x l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {Ifloat l .. Ifloat u} \<longrightarrow> Ifloat l1 \<le> x^n \<and> x^n \<le> Ifloat u1" 

204 
using float_power_bnds by auto 

205 

206 
section "Square root" 

207 

208 
text {* 

209 

210 
The square root computation is implemented as newton iteration. As first first step we use the 

211 
nearest power of two greater than the square root. 

212 

213 
*} 

214 

215 
fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where 

216 
"sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)"  

217 
"sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x 

218 
in Float 1 1 * (y + float_divr prec x y))" 

219 

220 
definition ub_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where 

221 
"ub_sqrt prec x = (if 0 < x then Some (sqrt_iteration prec prec x) else if x < 0 then None else Some 0)" 

222 

223 
definition lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where 

224 
"lb_sqrt prec x = (if 0 < x then Some (float_divl prec x (sqrt_iteration prec prec x)) else if x < 0 then None else Some 0)" 

225 

226 
lemma sqrt_ub_pos_pos_1: 

227 
assumes "sqrt x < b" and "0 < b" and "0 < x" 

228 
shows "sqrt x < (b + x / b)/2" 

229 
proof  

230 
from assms have "0 < (b  sqrt x) ^ 2 " by simp 

231 
also have "\<dots> = b ^ 2  2 * b * sqrt x + (sqrt x) ^ 2" by algebra 

232 
also have "\<dots> = b ^ 2  2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2) 

233 
finally have "0 < b ^ 2  2 * b * sqrt x + x" by assumption 

234 
hence "0 < b / 2  sqrt x + x / (2 * b)" using assms 

235 
by (simp add: field_simps power2_eq_square) 

236 
thus ?thesis by (simp add: field_simps) 

237 
qed 

238 

239 
lemma sqrt_iteration_bound: assumes "0 < Ifloat x" 

240 
shows "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec n x)" 

241 
proof (induct n) 

242 
case 0 

243 
show ?case 

244 
proof (cases x) 

245 
case (Float m e) 

246 
hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto 

247 
hence "0 < sqrt (real m)" by auto 

248 

249 
have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto 

250 

251 
have "Ifloat x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))" 

252 
unfolding pow2_add pow2_int Float Ifloat.simps by auto 

253 
also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))" 

254 
proof (rule mult_strict_right_mono, auto) 

255 
show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2] 

256 
unfolding real_of_int_less_iff[of m, symmetric] by auto 

257 
qed 

258 
finally have "sqrt (Ifloat x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto 

259 
also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)" 

260 
proof  

261 
let ?E = "e + bitlen m" 

262 
have E_mod_pow: "pow2 (?E mod 2) < 4" 

263 
proof (cases "?E mod 2 = 1") 

264 
case True thus ?thesis by auto 

265 
next 

266 
case False 

267 
have "0 \<le> ?E mod 2" by auto 

268 
have "?E mod 2 < 2" by auto 

269 
from this[THEN zless_imp_add1_zle] 

270 
have "?E mod 2 \<le> 0" using False by auto 

271 
from xt1(5)[OF `0 \<le> ?E mod 2` this] 

272 
show ?thesis by auto 

273 
qed 

274 
hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto 

275 
hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto 

276 

277 
have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto 

278 
have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))" 

279 
unfolding E_eq unfolding pow2_add .. 

280 
also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))" 

281 
unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto 

282 
also have "\<dots> < pow2 (?E div 2) * 2" 

283 
by (rule mult_strict_left_mono, auto intro: E_mod_pow) 

284 
also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto 

285 
finally show ?thesis by auto 

286 
qed 

287 
finally show ?thesis 

288 
unfolding Float sqrt_iteration.simps Ifloat.simps by auto 

289 
qed 

290 
next 

291 
case (Suc n) 

292 
let ?b = "sqrt_iteration prec n x" 

293 
have "0 < sqrt (Ifloat x)" using `0 < Ifloat x` by auto 

294 
also have "\<dots> < Ifloat ?b" using Suc . 

295 
finally have "sqrt (Ifloat x) < (Ifloat ?b + Ifloat x / Ifloat ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < Ifloat x`] by auto 

296 
also have "\<dots> \<le> (Ifloat ?b + Ifloat (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr) 

297 
also have "\<dots> = Ifloat (Float 1 1) * (Ifloat ?b + Ifloat (float_divr prec x ?b))" by auto 

298 
finally show ?case unfolding sqrt_iteration.simps Let_def Ifloat_mult Ifloat_add right_distrib . 

299 
qed 

300 

301 
lemma sqrt_iteration_lower_bound: assumes "0 < Ifloat x" 

302 
shows "0 < Ifloat (sqrt_iteration prec n x)" (is "0 < ?sqrt") 

303 
proof  

304 
have "0 < sqrt (Ifloat x)" using assms by auto 

305 
also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] . 

306 
finally show ?thesis . 

307 
qed 

308 

309 
lemma lb_sqrt_lower_bound: assumes "0 \<le> Ifloat x" 

310 
shows "0 \<le> Ifloat (the (lb_sqrt prec x))" 

311 
proof (cases "0 < x") 

312 
case True hence "0 < Ifloat x" and "0 \<le> x" using `0 \<le> Ifloat x` unfolding less_float_def le_float_def by auto 

313 
hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto 

314 
hence "0 \<le> Ifloat (float_divl prec x (sqrt_iteration prec prec x))" using float_divl_lower_bound[OF `0 \<le> x`] unfolding le_float_def by auto 

315 
thus ?thesis unfolding lb_sqrt_def using True by auto 

316 
next 

317 
case False with `0 \<le> Ifloat x` have "Ifloat x = 0" unfolding less_float_def by auto 

318 
thus ?thesis unfolding lb_sqrt_def less_float_def by auto 

319 
qed 

320 

321 
lemma lb_sqrt_upper_bound: assumes "0 \<le> Ifloat x" 

322 
shows "Ifloat (the (lb_sqrt prec x)) \<le> sqrt (Ifloat x)" 

323 
proof (cases "0 < x") 

324 
case True hence "0 < Ifloat x" and "0 \<le> Ifloat x" unfolding less_float_def by auto 

325 
hence sqrt_gt0: "0 < sqrt (Ifloat x)" by auto 

326 
hence sqrt_ub: "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto 

327 

328 
have "Ifloat (float_divl prec x (sqrt_iteration prec prec x)) \<le> Ifloat x / Ifloat (sqrt_iteration prec prec x)" by (rule float_divl) 

329 
also have "\<dots> < Ifloat x / sqrt (Ifloat x)" 

330 
by (rule divide_strict_left_mono[OF sqrt_ub `0 < Ifloat x` mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]]) 

331 
also have "\<dots> = sqrt (Ifloat x)" unfolding inverse_eq_iff_eq[of _ "sqrt (Ifloat x)", symmetric] sqrt_divide_self_eq[OF `0 \<le> Ifloat x`, symmetric] by auto 

332 
finally show ?thesis unfolding lb_sqrt_def if_P[OF `0 < x`] by auto 

333 
next 

334 
case False with `0 \<le> Ifloat x` 

335 
have "\<not> x < 0" unfolding less_float_def le_float_def by auto 

336 
show ?thesis unfolding lb_sqrt_def if_not_P[OF False] if_not_P[OF `\<not> x < 0`] using assms by auto 

337 
qed 

338 

339 
lemma lb_sqrt: assumes "Some y = lb_sqrt prec x" 

340 
shows "Ifloat y \<le> sqrt (Ifloat x)" and "0 \<le> Ifloat x" 

341 
proof  

342 
show "0 \<le> Ifloat x" 

343 
proof (rule ccontr) 

344 
assume "\<not> 0 \<le> Ifloat x" 

345 
hence "lb_sqrt prec x = None" unfolding lb_sqrt_def less_float_def by auto 

346 
thus False using assms by auto 

347 
qed 

348 
from lb_sqrt_upper_bound[OF this, of prec] 

349 
show "Ifloat y \<le> sqrt (Ifloat x)" unfolding assms[symmetric] by auto 

350 
qed 

351 

352 
lemma ub_sqrt_lower_bound: assumes "0 \<le> Ifloat x" 

353 
shows "sqrt (Ifloat x) \<le> Ifloat (the (ub_sqrt prec x))" 

354 
proof (cases "0 < x") 

355 
case True hence "0 < Ifloat x" unfolding less_float_def by auto 

356 
hence "0 < sqrt (Ifloat x)" by auto 

357 
hence "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto 

358 
thus ?thesis unfolding ub_sqrt_def if_P[OF `0 < x`] by auto 

359 
next 

360 
case False with `0 \<le> Ifloat x` 

361 
have "Ifloat x = 0" unfolding less_float_def le_float_def by auto 

362 
thus ?thesis unfolding ub_sqrt_def less_float_def le_float_def by auto 

363 
qed 

364 

365 
lemma ub_sqrt: assumes "Some y = ub_sqrt prec x" 

366 
shows "sqrt (Ifloat x) \<le> Ifloat y" and "0 \<le> Ifloat x" 

367 
proof  

368 
show "0 \<le> Ifloat x" 

369 
proof (rule ccontr) 

370 
assume "\<not> 0 \<le> Ifloat x" 

371 
hence "ub_sqrt prec x = None" unfolding ub_sqrt_def less_float_def by auto 

372 
thus False using assms by auto 

373 
qed 

374 
from ub_sqrt_lower_bound[OF this, of prec] 

375 
show "sqrt (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto 

376 
qed 

377 

378 
lemma bnds_sqrt: "\<forall> x lx ux. (Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u" 

379 
proof (rule allI, rule allI, rule allI, rule impI) 

380 
fix x lx ux 

381 
assume "(Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}" 

382 
hence l: "Some l = lb_sqrt prec lx " and u: "Some u = ub_sqrt prec ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto 

383 

384 
have "Ifloat lx \<le> x" and "x \<le> Ifloat ux" using x by auto 

385 

386 
from lb_sqrt(1)[OF l] real_sqrt_le_mono[OF `Ifloat lx \<le> x`] 

387 
have "Ifloat l \<le> sqrt x" by (rule order_trans) 

388 
moreover 

389 
from real_sqrt_le_mono[OF `x \<le> Ifloat ux`] ub_sqrt(1)[OF u] 

390 
have "sqrt x \<le> Ifloat u" by (rule order_trans) 

391 
ultimately show "Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u" .. 

392 
qed 

393 

394 
section "Arcus tangens and \<pi>" 

395 

396 
subsection "Compute arcus tangens series" 

397 

398 
text {* 

399 

400 
As first step we implement the computation of the arcus tangens series. This is only valid in the range 

401 
@{term "{1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens. 

402 

403 
*} 

404 

405 
fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" 

406 
and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where 

407 
"ub_arctan_horner prec 0 k x = 0" 

408 
 "ub_arctan_horner prec (Suc n) k x = 

409 
(rapprox_rat prec 1 (int k))  x * (lb_arctan_horner prec n (k + 2) x)" 

410 
 "lb_arctan_horner prec 0 k x = 0" 

411 
 "lb_arctan_horner prec (Suc n) k x = 

412 
(lapprox_rat prec 1 (int k))  x * (ub_arctan_horner prec n (k + 2) x)" 

413 

414 
lemma arctan_0_1_bounds': assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1" and "even n" 

415 
shows "arctan (Ifloat x) \<in> {Ifloat (x * lb_arctan_horner prec n 1 (x * x)) .. Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x * x))}" 

416 
proof  

417 
let "?c i" = "1^i * (1 / real (i * 2 + 1) * Ifloat x ^ (i * 2 + 1))" 

418 
let "?S n" = "\<Sum> i=0..<n. ?c i" 

419 

420 
have "0 \<le> Ifloat (x * x)" by auto 

421 
from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto 

422 

423 
have "arctan (Ifloat x) \<in> { ?S n .. ?S (Suc n) }" 

424 
proof (cases "Ifloat x = 0") 

425 
case False 

426 
hence "0 < Ifloat x" using `0 \<le> Ifloat x` by auto 

427 
hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * Ifloat x ^ (0 * 2 + 1)" by auto 

428 

429 
have "\<bar> Ifloat x \<bar> \<le> 1" using `0 \<le> Ifloat x` `Ifloat x \<le> 1` by auto 

430 
from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this] monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded `2 * m = n`] 

431 
show ?thesis unfolding arctan_series[OF `\<bar> Ifloat x \<bar> \<le> 1`] Suc_plus1 . 

432 
qed auto 

433 
note arctan_bounds = this[unfolded atLeastAtMost_iff] 

434 

435 
have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto 

436 

437 
note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0 

438 
and lb="\<lambda>n i k x. lb_arctan_horner prec n k x" 

439 
and ub="\<lambda>n i k x. ub_arctan_horner prec n k x", 

440 
OF `0 \<le> Ifloat (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps] 

441 

442 
{ have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n" 

443 
using bounds(1) `0 \<le> Ifloat x` 

444 
unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric] 

445 
unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"] 

446 
by (auto intro!: mult_left_mono) 

447 
also have "\<dots> \<le> arctan (Ifloat x)" using arctan_bounds .. 

448 
finally have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (Ifloat x)" . } 

449 
moreover 

450 
{ have "arctan (Ifloat x) \<le> ?S (Suc n)" using arctan_bounds .. 

451 
also have "\<dots> \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))" 

452 
using bounds(2)[of "Suc n"] `0 \<le> Ifloat x` 

453 
unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric] 

454 
unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"] 

455 
by (auto intro!: mult_left_mono) 

456 
finally have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . } 

457 
ultimately show ?thesis by auto 

458 
qed 

459 

460 
lemma arctan_0_1_bounds: assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1" 

461 
shows "arctan (Ifloat x) \<in> {Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}" 

462 
proof (cases "even n") 

463 
case True 

464 
obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto 

465 
hence "even n'" unfolding even_nat_Suc by auto 

466 
have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))" 

467 
unfolding `Suc n' = get_odd n`[symmetric] using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto 

468 
moreover 

469 
have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)" 

470 
unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n`] by auto 

471 
ultimately show ?thesis by auto 

472 
next 

473 
case False hence "0 < n" by (rule odd_pos) 

474 
from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" .. 

475 
from False[unfolded this even_nat_Suc] 

476 
have "even n'" and "even (Suc (Suc n'))" by auto 

477 
have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` . 

478 

479 
have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))" 

480 
unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto 

481 
moreover 

482 
have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)" 

483 
unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even (Suc (Suc n'))`] by auto 

484 
ultimately show ?thesis by auto 

485 
qed 

486 

487 
subsection "Compute \<pi>" 

488 

489 
definition ub_pi :: "nat \<Rightarrow> float" where 

490 
"ub_pi prec = (let A = rapprox_rat prec 1 5 ; 

491 
B = lapprox_rat prec 1 239 

492 
in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A))  

493 
B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))" 

494 

495 
definition lb_pi :: "nat \<Rightarrow> float" where 

496 
"lb_pi prec = (let A = lapprox_rat prec 1 5 ; 

497 
B = rapprox_rat prec 1 239 

498 
in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A))  

499 
B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))" 

500 

501 
lemma pi_boundaries: "pi \<in> {Ifloat (lb_pi n) .. Ifloat (ub_pi n)}" 

502 
proof  

503 
have machin_pi: "pi = 4 * (4 * arctan (1 / 5)  arctan (1 / 239))" unfolding machin[symmetric] by auto 

504 

505 
{ fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto 

506 
let ?k = "rapprox_rat prec 1 k" 

507 
have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto 

508 

509 
have "0 \<le> Ifloat ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`) 

510 
have "Ifloat ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`] 

511 
by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`) 

512 

513 
have "1 / real k \<le> Ifloat ?k" using rapprox_rat[where x=1 and y=k] by auto 

514 
hence "arctan (1 / real k) \<le> arctan (Ifloat ?k)" by (rule arctan_monotone') 

515 
also have "\<dots> \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" 

516 
using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto 

517 
finally have "arctan (1 / (real k)) \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" . 

518 
} note ub_arctan = this 

519 

520 
{ fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto 

521 
let ?k = "lapprox_rat prec 1 k" 

522 
have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto 

523 
have "1 / real k \<le> 1" using `1 < k` by auto 

524 

525 
have "\<And>n. 0 \<le> Ifloat ?k" using lapprox_rat_bottom[where x=1 and y=k, OF zero_le_one `0 < k`] by (auto simp add: `1 div k = 0`) 

526 
have "\<And>n. Ifloat ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`) 

527 

528 
have "Ifloat ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto 

529 

530 
have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (Ifloat ?k)" 

531 
using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto 

532 
also have "\<dots> \<le> arctan (1 / real k)" using `Ifloat ?k \<le> 1 / real k` by (rule arctan_monotone') 

533 
finally have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" . 

534 
} note lb_arctan = this 

535 

536 
have "pi \<le> Ifloat (ub_pi n)" 

537 
unfolding ub_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub unfolding Float_num 

538 
using lb_arctan[of 239] ub_arctan[of 5] 

539 
by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps) 

540 
moreover 

541 
have "Ifloat (lb_pi n) \<le> pi" 

542 
unfolding lb_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub Float_num 

543 
using lb_arctan[of 5] ub_arctan[of 239] 

544 
by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps) 

545 
ultimately show ?thesis by auto 

546 
qed 

547 

548 
subsection "Compute arcus tangens in the entire domain" 

549 

550 
function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where 

551 
"lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ; 

552 
lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) 

553 
in (if x < 0 then  ub_arctan prec (x) else 

554 
if x \<le> Float 1 1 then lb_horner x else 

555 
if x \<le> Float 1 1 then Float 1 1 * lb_horner (float_divl prec x (1 + the (ub_sqrt prec (1 + x * x)))) 

556 
else (let inv = float_divr prec 1 x 

557 
in if inv > 1 then 0 

558 
else lb_pi prec * Float 1 1  ub_horner inv)))" 

559 

560 
 "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ; 

561 
ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) 

562 
in (if x < 0 then  lb_arctan prec (x) else 

563 
if x \<le> Float 1 1 then ub_horner x else 

564 
if x \<le> Float 1 1 then let y = float_divr prec x (1 + the (lb_sqrt prec (1 + x * x))) 

565 
in if y > 1 then ub_pi prec * Float 1 1 

566 
else Float 1 1 * ub_horner y 

567 
else ub_pi prec * Float 1 1  lb_horner (float_divl prec 1 x)))" 

568 
by pat_completeness auto 

569 
termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def) 

570 

571 
declare ub_arctan_horner.simps[simp del] 

572 
declare lb_arctan_horner.simps[simp del] 

573 

574 
lemma lb_arctan_bound': assumes "0 \<le> Ifloat x" 

575 
shows "Ifloat (lb_arctan prec x) \<le> arctan (Ifloat x)" 

576 
proof  

577 
have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto 

578 
let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)" 

579 
and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)" 

580 

581 
show ?thesis 

582 
proof (cases "x \<le> Float 1 1") 

583 
case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto 

584 
show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True] 

585 
using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto 

586 
next 

587 
case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto 

588 
let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)" 

589 
let ?fR = "1 + the (ub_sqrt prec (1 + x * x))" 

590 
let ?DIV = "float_divl prec x ?fR" 

591 

592 
have sqr_ge0: "0 \<le> 1 + Ifloat x * Ifloat x" using sum_power2_ge_zero[of 1 "Ifloat x", unfolded numeral_2_eq_2] by auto 

593 
hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg) 

594 

595 
have "sqrt (Ifloat (1 + x * x)) \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0) 

596 
hence "?R \<le> Ifloat ?fR" by auto 

597 
hence "0 < ?fR" and "0 < Ifloat ?fR" unfolding less_float_def using `0 < ?R` by auto 

598 

599 
have monotone: "Ifloat (float_divl prec x ?fR) \<le> Ifloat x / ?R" 

600 
proof  

601 
have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl) 

602 
also have "\<dots> \<le> Ifloat x / ?R" by (rule divide_left_mono[OF `?R \<le> Ifloat ?fR` `0 \<le> Ifloat x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> Ifloat ?fR`] divisor_gt0]]) 

603 
finally show ?thesis . 

604 
qed 

605 

606 
show ?thesis 

607 
proof (cases "x \<le> Float 1 1") 

608 
case True 

609 

610 
have "Ifloat x \<le> sqrt (Ifloat (1 + x * x))" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto 

611 
also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0) 

612 
finally have "Ifloat x \<le> Ifloat ?fR" by auto 

613 
moreover have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl) 

614 
ultimately have "Ifloat ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < Ifloat ?fR`, symmetric] by auto 

615 

616 
have "0 \<le> Ifloat ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto 

617 

618 
have "Ifloat (Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (Ifloat (float_divl prec x ?fR))" unfolding Ifloat_mult[of "Float 1 1"] Float_num 

619 
using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto 

620 
also have "\<dots> \<le> 2 * arctan (Ifloat x / ?R)" 

621 
using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono) 

30273
ecd6f0ca62ea
declare power_Suc [simp]; remove redundant typespecific versions of power_Suc
huffman
parents:
30122
diff
changeset

622 
also have "2 * arctan (Ifloat x / ?R) = arctan (Ifloat x)" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 . 
29805  623 
finally show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_P[OF True] . 
624 
next 

625 
case False 

626 
hence "2 < Ifloat x" unfolding le_float_def Float_num by auto 

627 
hence "1 \<le> Ifloat x" by auto 

628 

629 
let "?invx" = "float_divr prec 1 x" 

630 
have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto 

631 

632 
show ?thesis 

633 
proof (cases "1 < ?invx") 

634 
case True 

635 
show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False] if_P[OF True] 

636 
using `0 \<le> arctan (Ifloat x)` by auto 

637 
next 

638 
case False 

639 
hence "Ifloat ?invx \<le> 1" unfolding less_float_def by auto 

640 
have "0 \<le> Ifloat ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> Ifloat x`) 

641 

642 
have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto 

643 

644 
have "arctan (1 / Ifloat x) \<le> arctan (Ifloat ?invx)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divr) 

645 
also have "\<dots> \<le> Ifloat (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto 

646 
finally have "pi / 2  Ifloat (?ub_horner ?invx) \<le> arctan (Ifloat x)" 

647 
using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`] 

648 
unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto 

649 
moreover 

650 
have "Ifloat (lb_pi prec * Float 1 1) \<le> pi / 2" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto 

651 
ultimately 

652 
show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False] 

653 
by auto 

654 
qed 

655 
qed 

656 
qed 

657 
qed 

658 

659 
lemma ub_arctan_bound': assumes "0 \<le> Ifloat x" 

660 
shows "arctan (Ifloat x) \<le> Ifloat (ub_arctan prec x)" 

661 
proof  

662 
have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto 

663 

664 
let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)" 

665 
and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)" 

666 

667 
show ?thesis 

668 
proof (cases "x \<le> Float 1 1") 

669 
case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto 

670 
show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True] 

671 
using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto 

672 
next 

673 
case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto 

674 
let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)" 

675 
let ?fR = "1 + the (lb_sqrt prec (1 + x * x))" 

676 
let ?DIV = "float_divr prec x ?fR" 

677 

678 
have sqr_ge0: "0 \<le> 1 + Ifloat x * Ifloat x" using sum_power2_ge_zero[of 1 "Ifloat x", unfolded numeral_2_eq_2] by auto 

679 
hence "0 \<le> Ifloat (1 + x*x)" by auto 

680 

681 
hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg) 

682 

683 
have "Ifloat (the (lb_sqrt prec (1 + x * x))) \<le> sqrt (Ifloat (1 + x * x))" by (rule lb_sqrt_upper_bound, auto simp add: sqr_ge0) 

684 
hence "Ifloat ?fR \<le> ?R" by auto 

685 
have "0 < Ifloat ?fR" unfolding Ifloat_add Ifloat_1 by (rule order_less_le_trans[OF zero_less_one], auto simp add: lb_sqrt_lower_bound[OF `0 \<le> Ifloat (1 + x*x)`]) 

686 

687 
have monotone: "Ifloat x / ?R \<le> Ifloat (float_divr prec x ?fR)" 

688 
proof  

689 
from divide_left_mono[OF `Ifloat ?fR \<le> ?R` `0 \<le> Ifloat x` mult_pos_pos[OF divisor_gt0 `0 < Ifloat ?fR`]] 

690 
have "Ifloat x / ?R \<le> Ifloat x / Ifloat ?fR" . 

691 
also have "\<dots> \<le> Ifloat ?DIV" by (rule float_divr) 

692 
finally show ?thesis . 

693 
qed 

694 

695 
show ?thesis 

696 
proof (cases "x \<le> Float 1 1") 

697 
case True 

698 
show ?thesis 

699 
proof (cases "?DIV > 1") 

700 
case True 

701 
have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 1)" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto 

702 
from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le] 

703 
show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_P[OF `x \<le> Float 1 1`] if_P[OF True] . 

704 
next 

705 
case False 

706 
hence "Ifloat ?DIV \<le> 1" unfolding less_float_def by auto 

707 

708 
have "0 \<le> Ifloat x / ?R" using `0 \<le> Ifloat x` `0 < ?R` unfolding real_0_le_divide_iff by auto 

709 
hence "0 \<le> Ifloat ?DIV" using monotone by (rule order_trans) 

710 

30273
ecd6f0ca62ea
declare power_Suc [simp]; remove redundant typespecific versions of power_Suc
huffman
parents:
30122
diff
changeset

711 
have "arctan (Ifloat x) = 2 * arctan (Ifloat x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 . 
29805  712 
also have "\<dots> \<le> 2 * arctan (Ifloat ?DIV)" 
713 
using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono) 

714 
also have "\<dots> \<le> Ifloat (Float 1 1 * ?ub_horner ?DIV)" unfolding Ifloat_mult[of "Float 1 1"] Float_num 

715 
using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto 

716 
finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_P[OF `x \<le> Float 1 1`] if_not_P[OF False] . 

717 
qed 

718 
next 

719 
case False 

720 
hence "2 < Ifloat x" unfolding le_float_def Float_num by auto 

721 
hence "1 \<le> Ifloat x" by auto 

722 
hence "0 < Ifloat x" by auto 

723 
hence "0 < x" unfolding less_float_def by auto 

724 

725 
let "?invx" = "float_divl prec 1 x" 

726 
have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto 

727 

728 
have "Ifloat ?invx \<le> 1" unfolding less_float_def by (rule order_trans[OF float_divl], auto simp add: `1 \<le> Ifloat x` divide_le_eq_1_pos[OF `0 < Ifloat x`]) 

729 
have "0 \<le> Ifloat ?invx" unfolding Ifloat_0[symmetric] by (rule float_divl_lower_bound[unfolded le_float_def], auto simp add: `0 < x`) 

730 

731 
have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto 

732 

733 
have "Ifloat (?lb_horner ?invx) \<le> arctan (Ifloat ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto 

734 
also have "\<dots> \<le> arctan (1 / Ifloat x)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divl) 

735 
finally have "arctan (Ifloat x) \<le> pi / 2  Ifloat (?lb_horner ?invx)" 

736 
using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`] 

737 
unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto 

738 
moreover 

739 
have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 1)" unfolding Ifloat_mult Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto 

740 
ultimately 

741 
show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False] 

742 
by auto 

743 
qed 

744 
qed 

745 
qed 

746 

747 
lemma arctan_boundaries: 

748 
"arctan (Ifloat x) \<in> {Ifloat (lb_arctan prec x) .. Ifloat (ub_arctan prec x)}" 

749 
proof (cases "0 \<le> x") 

750 
case True hence "0 \<le> Ifloat x" unfolding le_float_def by auto 

751 
show ?thesis using ub_arctan_bound'[OF `0 \<le> Ifloat x`] lb_arctan_bound'[OF `0 \<le> Ifloat x`] unfolding atLeastAtMost_iff by auto 

752 
next 

753 
let ?mx = "x" 

754 
case False hence "x < 0" and "0 \<le> Ifloat ?mx" unfolding le_float_def less_float_def by auto 

755 
hence bounds: "Ifloat (lb_arctan prec ?mx) \<le> arctan (Ifloat ?mx) \<and> arctan (Ifloat ?mx) \<le> Ifloat (ub_arctan prec ?mx)" 

756 
using ub_arctan_bound'[OF `0 \<le> Ifloat ?mx`] lb_arctan_bound'[OF `0 \<le> Ifloat ?mx`] by auto 

757 
show ?thesis unfolding Ifloat_minus arctan_minus lb_arctan.simps[where x=x] ub_arctan.simps[where x=x] Let_def if_P[OF `x < 0`] 

758 
unfolding atLeastAtMost_iff using bounds[unfolded Ifloat_minus arctan_minus] by auto 

759 
qed 

760 

761 
lemma bnds_arctan: "\<forall> x lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u" 

762 
proof (rule allI, rule allI, rule allI, rule impI) 

763 
fix x lx ux 

764 
assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}" 

765 
hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto 

766 

767 
{ from arctan_boundaries[of lx prec, unfolded l] 

768 
have "Ifloat l \<le> arctan (Ifloat lx)" by (auto simp del: lb_arctan.simps) 

769 
also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone') 

770 
finally have "Ifloat l \<le> arctan x" . 

771 
} moreover 

772 
{ have "arctan x \<le> arctan (Ifloat ux)" using x by (auto intro: arctan_monotone') 

773 
also have "\<dots> \<le> Ifloat u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps) 

774 
finally have "arctan x \<le> Ifloat u" . 

775 
} ultimately show "Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u" .. 

776 
qed 

777 

778 
section "Sinus and Cosinus" 

779 

780 
subsection "Compute the cosinus and sinus series" 

781 

782 
fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" 

783 
and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where 

784 
"ub_sin_cos_aux prec 0 i k x = 0" 

785 
 "ub_sin_cos_aux prec (Suc n) i k x = 

786 
(rapprox_rat prec 1 (int k))  x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)" 

787 
 "lb_sin_cos_aux prec 0 i k x = 0" 

788 
 "lb_sin_cos_aux prec (Suc n) i k x = 

789 
(lapprox_rat prec 1 (int k))  x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)" 

790 

791 
lemma cos_aux: 

792 
shows "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. 1^i * (1/real (fact (2 * i))) * (Ifloat x)^(2 * i))" (is "?lb") 

793 
and "(\<Sum> i=0..<n. 1^i * (1/real (fact (2 * i))) * (Ifloat x)^(2 * i)) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub") 

794 
proof  

795 
have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto 

796 
let "?f n" = "fact (2 * n)" 

797 

798 
{ fix n 

799 
have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto) 

800 
have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 1 * (((\<lambda>i. i + 2) ^ n) 1 + 1)" 

801 
unfolding F by auto } note f_eq = this 

802 

803 
from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, 

804 
OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps] 

805 
show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "Ifloat x"]) 

806 
qed 

807 

808 
lemma cos_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2" 

809 
shows "cos (Ifloat x) \<in> {Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. Ifloat (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}" 

810 
proof (cases "Ifloat x = 0") 

811 
case False hence "Ifloat x \<noteq> 0" by auto 

812 
hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto 

813 
have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0 

814 
using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto 

815 

816 
{ fix x n have "(\<Sum> i=0..<n. 1^i * (1/real (fact (2 * i))) * x^(2 * i)) 

817 
= (\<Sum> i = 0 ..< 2 * n. (if even(i) then (1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum") 

818 
proof  

819 
have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto 

820 
also have "\<dots> = 

821 
(\<Sum> j = 0 ..< n. 1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto 

822 
also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)" 

823 
unfolding sum_split_even_odd .. 

824 
also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)" 

825 
by (rule setsum_cong2) auto 

826 
finally show ?thesis by assumption 

827 
qed } note morph_to_if_power = this 

828 

829 

830 
{ fix n :: nat assume "0 < n" 

831 
hence "0 < 2 * n" by auto 

832 
obtain t where "0 < t" and "t < Ifloat x" and 

833 
cos_eq: "cos (Ifloat x) = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (1 ^ (i div 2))/(real (fact i)) else 0) * (Ifloat x) ^ i) 

834 
+ (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (Ifloat x)^(2*n)" 

835 
(is "_ = ?SUM + ?rest / ?fact * ?pow") 

836 
using Maclaurin_cos_expansion2[OF `0 < Ifloat x` `0 < 2 * n`] by auto 

837 

838 
have "cos t * 1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto 

839 
also have "\<dots> = cos (t + real n * pi)" using cos_add by auto 

840 
also have "\<dots> = ?rest" by auto 

841 
finally have "cos t * 1^n = ?rest" . 

842 
moreover 

843 
have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto 

844 
hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto 

845 
ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le>  ?rest " by auto 

846 

847 
have "0 < ?fact" by auto 

848 
have "0 < ?pow" using `0 < Ifloat x` by auto 

849 

850 
{ 

851 
assume "even n" 

852 
have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM" 

853 
unfolding morph_to_if_power[symmetric] using cos_aux by auto 

854 
also have "\<dots> \<le> cos (Ifloat x)" 

855 
proof  

856 
from even[OF `even n`] `0 < ?fact` `0 < ?pow` 

857 
have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le) 

858 
thus ?thesis unfolding cos_eq by auto 

859 
qed 

860 
finally have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (Ifloat x)" . 

861 
} note lb = this 

862 

863 
{ 

864 
assume "odd n" 

865 
have "cos (Ifloat x) \<le> ?SUM" 

866 
proof  

867 
from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`] 

868 
have "0 \<le> ( ?rest) / ?fact * ?pow" 

869 
by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le) 

870 
thus ?thesis unfolding cos_eq by auto 

871 
qed 

872 
also have "\<dots> \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" 

873 
unfolding morph_to_if_power[symmetric] using cos_aux by auto 

874 
finally have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" . 

875 
} note ub = this and lb 

876 
} note ub = this(1) and lb = this(2) 

877 

878 
have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] . 

879 
moreover have "Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (Ifloat x)" 

880 
proof (cases "0 < get_even n") 

881 
case True show ?thesis using lb[OF True get_even] . 

882 
next 

883 
case False 

884 
hence "get_even n = 0" by auto 

885 
have " (pi / 2) \<le> Ifloat x" by (rule order_trans[OF _ `0 < Ifloat x`[THEN less_imp_le]], auto) 

886 
with `Ifloat x \<le> pi / 2` 

887 
show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps Ifloat_minus Ifloat_0 using cos_ge_zero by auto 

888 
qed 

889 
ultimately show ?thesis by auto 

890 
next 

891 
case True 

892 
show ?thesis 

893 
proof (cases "n = 0") 

894 
case True 

895 
thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat x = 0` lapprox_rat[where x="1" and y=1] by auto 

896 
next 

897 
case False with not0_implies_Suc obtain m where "n = Suc m" by blast 

898 
thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `Ifloat x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto) 

899 
qed 

900 
qed 

901 

902 
lemma sin_aux: assumes "0 \<le> Ifloat x" 

903 
shows "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> (\<Sum> i=0..<n. 1^i * (1/real (fact (2 * i + 1))) * (Ifloat x)^(2 * i + 1))" (is "?lb") 

904 
and "(\<Sum> i=0..<n. 1^i * (1/real (fact (2 * i + 1))) * (Ifloat x)^(2 * i + 1)) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub") 

905 
proof  

906 
have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto 

907 
let "?f n" = "fact (2 * n + 1)" 

908 

909 
{ fix n 

910 
have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto) 

911 
have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 2 * (((\<lambda>i. i + 2) ^ n) 2 + 1)" 

912 
unfolding F by auto } note f_eq = this 

913 

914 
from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, 

915 
OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps] 

916 
show "?lb" and "?ub" using `0 \<le> Ifloat x` unfolding Ifloat_mult 

917 
unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric] 

918 
unfolding real_mult_commute 

919 
by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "Ifloat x"]) 

920 
qed 

921 

922 
lemma sin_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2" 

923 
shows "sin (Ifloat x) \<in> {Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. Ifloat (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}" 

924 
proof (cases "Ifloat x = 0") 

925 
case False hence "Ifloat x \<noteq> 0" by auto 

926 
hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto 

927 
have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0 

928 
using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto 

929 

930 
{ fix x n have "(\<Sum> j = 0 ..< n. 1 ^ (((2 * j + 1)  Suc 0) div 2) / (real (fact (2 * j + 1))) * x ^(2 * j + 1)) 

931 
= (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (1 ^ ((i  Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _") 

932 
proof  

933 
have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto 

934 
have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto 

935 
also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else 1 ^ ((i  Suc 0) div 2) / (real (fact i)) * x ^ i)" 

936 
unfolding sum_split_even_odd .. 

937 
also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 0 else 1 ^ ((i  Suc 0) div 2) / (real (fact i))) * x ^ i)" 

938 
by (rule setsum_cong2) auto 

939 
finally show ?thesis by assumption 

940 
qed } note setsum_morph = this 

941 

942 
{ fix n :: nat assume "0 < n" 

943 
hence "0 < 2 * n + 1" by auto 

944 
obtain t where "0 < t" and "t < Ifloat x" and 

945 
sin_eq: "sin (Ifloat x) = (\<Sum> i = 0 ..< 2 * n + 1. (if even(i) then 0 else (1 ^ ((i  Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i) 

946 
+ (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (Ifloat x)^(2*n + 1)" 

947 
(is "_ = ?SUM + ?rest / ?fact * ?pow") 

948 
using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < Ifloat x`] by auto 

949 

950 
have "?rest = cos t * 1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto 

951 
moreover 

952 
have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto 

953 
hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto 

954 
ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le>  ?rest " by auto 

955 

956 
have "0 < ?fact" by (rule real_of_nat_fact_gt_zero) 

957 
have "0 < ?pow" using `0 < Ifloat x` by (rule zero_less_power) 

958 

959 
{ 

960 
assume "even n" 

961 
have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> 

962 
(\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (1 ^ ((i  Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)" 

963 
using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto 

964 
also have "\<dots> \<le> ?SUM" by auto 

965 
also have "\<dots> \<le> sin (Ifloat x)" 

966 
proof  

967 
from even[OF `even n`] `0 < ?fact` `0 < ?pow` 

968 
have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le) 

969 
thus ?thesis unfolding sin_eq by auto 

970 
qed 

971 
finally have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (Ifloat x)" . 

972 
} note lb = this 

973 

974 
{ 

975 
assume "odd n" 

976 
have "sin (Ifloat x) \<le> ?SUM" 

977 
proof  

978 
from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`] 

979 
have "0 \<le> ( ?rest) / ?fact * ?pow" 

980 
by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le) 

981 
thus ?thesis unfolding sin_eq by auto 

982 
qed 

983 
also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (1 ^ ((i  Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)" 

984 
by auto 

985 
also have "\<dots> \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" 

986 
using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto 

987 
finally have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" . 

988 
} note ub = this and lb 

989 
} note ub = this(1) and lb = this(2) 

990 

991 
have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] . 

992 
moreover have "Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (Ifloat x)" 

993 
proof (cases "0 < get_even n") 

994 
case True show ?thesis using lb[OF True get_even] . 

995 
next 

996 
case False 

997 
hence "get_even n = 0" by auto 

998 
with `Ifloat x \<le> pi / 2` `0 \<le> Ifloat x` 

999 
show ?thesis unfolding `get_even n = 0` ub_sin_cos_aux.simps Ifloat_minus Ifloat_0 using sin_ge_zero by auto 

1000 
qed 

1001 
ultimately show ?thesis by auto 

1002 
next 

1003 
case True 

1004 
show ?thesis 

1005 
proof (cases "n = 0") 

1006 
case True 

1007 
thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat x = 0` lapprox_rat[where x="1" and y=1] by auto 

1008 
next 

1009 
case False with not0_implies_Suc obtain m where "n = Suc m" by blast 

1010 
thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `Ifloat x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto) 

1011 
qed 

1012 
qed 

1013 

1014 
subsection "Compute the cosinus in the entire domain" 

1015 

1016 
definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where 

1017 
"lb_cos prec x = (let 

1018 
horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ; 

1019 
half = \<lambda> x. if x < 0 then  1 else Float 1 1 * x * x  1 

1020 
in if x < Float 1 1 then horner x 

1021 
else if x < 1 then half (horner (x * Float 1 1)) 

1022 
else half (half (horner (x * Float 1 2))))" 

1023 

1024 
definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where 

1025 
"ub_cos prec x = (let 

1026 
horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ; 

1027 
half = \<lambda> x. Float 1 1 * x * x  1 

1028 
in if x < Float 1 1 then horner x 

1029 
else if x < 1 then half (horner (x * Float 1 1)) 

1030 
else half (half (horner (x * Float 1 2))))" 

1031 

1032 
definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where 

1033 
"bnds_cos prec lx ux = (let lpi = lb_pi prec 

1034 
in if lx < lpi \<or> ux > lpi then (Float 1 0, Float 1 0) 

1035 
else if ux \<le> 0 then (lb_cos prec (lx), ub_cos prec (ux)) 

1036 
else if 0 \<le> lx then (lb_cos prec ux, ub_cos prec lx) 

1037 
else (min (lb_cos prec (lx)) (lb_cos prec ux), Float 1 0))" 

1038 

1039 
lemma lb_cos: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi" 

1040 
shows "cos (Ifloat x) \<in> {Ifloat (lb_cos prec x) .. Ifloat (ub_cos prec x)}" (is "?cos x \<in> { Ifloat (?lb x) .. Ifloat (?ub x) }") 

1041 
proof  

1042 
{ fix x :: real 

1043 
have "cos x = cos (x / 2 + x / 2)" by auto 

1044 
also have "\<dots> = cos (x / 2) * cos (x / 2) + sin (x / 2) * sin (x / 2)  sin (x / 2) * sin (x / 2) + cos (x / 2) * cos (x / 2)  1" 

1045 
unfolding cos_add by auto 

1046 
also have "\<dots> = 2 * cos (x / 2) * cos (x / 2)  1" by algebra 

1047 
finally have "cos x = 2 * cos (x / 2) * cos (x / 2)  1" . 

1048 
} note x_half = this[symmetric] 

1049 

1050 
have "\<not> x < 0" using `0 \<le> Ifloat x` unfolding less_float_def by auto 

1051 
let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)" 

1052 
let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)" 

1053 
let "?ub_half x" = "Float 1 1 * x * x  1" 

1054 
let "?lb_half x" = "if x < 0 then  1 else Float 1 1 * x * x  1" 

1055 

1056 
show ?thesis 

1057 
proof (cases "x < Float 1 1") 

1058 
case True hence "Ifloat x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto 

1059 
show ?thesis unfolding lb_cos_def[where x=x] ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_P[OF `x < Float 1 1`] Let_def 

1060 
using cos_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`] . 

1061 
next 

1062 
case False 

1063 

1064 
{ fix y x :: float let ?x2 = "Ifloat (x * Float 1 1)" 

1065 
assume "Ifloat y \<le> cos ?x2" and "pi \<le> Ifloat x" and "Ifloat x \<le> pi" 

1066 
hence " (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto 

1067 
hence "0 \<le> cos ?x2" by (rule cos_ge_zero) 

1068 

1069 
have "Ifloat (?lb_half y) \<le> cos (Ifloat x)" 

1070 
proof (cases "y < 0") 

1071 
case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto 

1072 
next 

1073 
case False 

1074 
hence "0 \<le> Ifloat y" unfolding less_float_def by auto 

1075 
from mult_mono[OF `Ifloat y \<le> cos ?x2` `Ifloat y \<le> cos ?x2` `0 \<le> cos ?x2` this] 

1076 
have "Ifloat y * Ifloat y \<le> cos ?x2 * cos ?x2" . 

1077 
hence "2 * Ifloat y * Ifloat y \<le> 2 * cos ?x2 * cos ?x2" by auto 

1078 
hence "2 * Ifloat y * Ifloat y  1 \<le> 2 * cos (Ifloat x / 2) * cos (Ifloat x / 2)  1" unfolding Float_num Ifloat_mult by auto 

1079 
thus ?thesis unfolding if_not_P[OF False] x_half Float_num Ifloat_mult Ifloat_sub by auto 

1080 
qed 

1081 
} note lb_half = this 

1082 

1083 
{ fix y x :: float let ?x2 = "Ifloat (x * Float 1 1)" 

1084 
assume ub: "cos ?x2 \<le> Ifloat y" and " pi \<le> Ifloat x" and "Ifloat x \<le> pi" 

1085 
hence " (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto 

1086 
hence "0 \<le> cos ?x2" by (rule cos_ge_zero) 

1087 

1088 
have "cos (Ifloat x) \<le> Ifloat (?ub_half y)" 

1089 
proof  

1090 
have "0 \<le> Ifloat y" using `0 \<le> cos ?x2` ub by (rule order_trans) 

1091 
from mult_mono[OF ub ub this `0 \<le> cos ?x2`] 

1092 
have "cos ?x2 * cos ?x2 \<le> Ifloat y * Ifloat y" . 

1093 
hence "2 * cos ?x2 * cos ?x2 \<le> 2 * Ifloat y * Ifloat y" by auto 

1094 
hence "2 * cos (Ifloat x / 2) * cos (Ifloat x / 2)  1 \<le> 2 * Ifloat y * Ifloat y  1" unfolding Float_num Ifloat_mult by auto 

1095 
thus ?thesis unfolding x_half Ifloat_mult Float_num Ifloat_sub by auto 

1096 
qed 

1097 
} note ub_half = this 

1098 

1099 
let ?x2 = "x * Float 1 1" 

1100 
let ?x4 = "x * Float 1 1 * Float 1 1" 

1101 

1102 
have "pi \<le> Ifloat x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] `0 \<le> Ifloat x` by (rule order_trans) 

1103 

1104 
show ?thesis 

1105 
proof (cases "x < 1") 

1106 
case True hence "Ifloat x \<le> 1" unfolding less_float_def by auto 

1107 
have "0 \<le> Ifloat ?x2" and "Ifloat ?x2 \<le> pi / 2" using pi_ge_two `0 \<le> Ifloat x` unfolding Ifloat_mult Float_num using assms by auto 

1108 
from cos_boundaries[OF this] 

1109 
have lb: "Ifloat (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> Ifloat (?ub_horner ?x2)" by auto 

1110 

1111 
have "Ifloat (?lb x) \<le> ?cos x" 

1112 
proof  

1113 
from lb_half[OF lb `pi \<le> Ifloat x` `Ifloat x \<le> pi`] 

1114 
show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 1` `x < 1` by auto 

1115 
qed 

1116 
moreover have "?cos x \<le> Ifloat (?ub x)" 

1117 
proof  

1118 
from ub_half[OF ub `pi \<le> Ifloat x` `Ifloat x \<le> pi`] 

1119 
show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 1` `x < 1` by auto 

1120 
qed 

1121 
ultimately show ?thesis by auto 

1122 
next 

1123 
case False 

1124 
have "0 \<le> Ifloat ?x4" and "Ifloat ?x4 \<le> pi / 2" using pi_ge_two `0 \<le> Ifloat x` `Ifloat x \<le> pi` unfolding Ifloat_mult Float_num by auto 

1125 
from cos_boundaries[OF this] 

1126 
have lb: "Ifloat (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> Ifloat (?ub_horner ?x4)" by auto 

1127 

1128 
have eq_4: "?x2 * Float 1 1 = x * Float 1 2" by (cases x, auto simp add: times_float.simps) 

1129 

1130 
have "Ifloat (?lb x) \<le> ?cos x" 

1131 
proof  

1132 
have "pi \<le> Ifloat ?x2" and "Ifloat ?x2 \<le> pi" unfolding Ifloat_mult Float_num using pi_ge_two `0 \<le> Ifloat x` `Ifloat x \<le> pi` by auto 

1133 
from lb_half[OF lb_half[OF lb this] `pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4] 

1134 
show ?thesis unfolding lb_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 1`] if_not_P[OF `\<not> x < 1`] Let_def . 

1135 
qed 

1136 
moreover have "?cos x \<le> Ifloat (?ub x)" 

1137 
proof  

1138 
have "pi \<le> Ifloat ?x2" and "Ifloat ?x2 \<le> pi" unfolding Ifloat_mult Float_num using pi_ge_two `0 \<le> Ifloat x` `Ifloat x \<le> pi` by auto 

1139 
from ub_half[OF ub_half[OF ub this] `pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4] 

1140 
show ?thesis unfolding ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 1`] if_not_P[OF `\<not> x < 1`] Let_def . 

1141 
qed 

1142 
ultimately show ?thesis by auto 

1143 
qed 

1144 
qed 

1145 
qed 

1146 

1147 
lemma lb_cos_minus: assumes "pi \<le> Ifloat x" and "Ifloat x \<le> 0" 

1148 
shows "cos (Ifloat (x)) \<in> {Ifloat (lb_cos prec (x)) .. Ifloat (ub_cos prec (x))}" 

1149 
proof  

1150 
have "0 \<le> Ifloat (x)" and "Ifloat (x) \<le> pi" using `pi \<le> Ifloat x` `Ifloat x \<le> 0` by auto 

1151 
from lb_cos[OF this] show ?thesis . 

1152 
qed 

1153 

1154 
lemma bnds_cos: "\<forall> x lx ux. (l, u) = bnds_cos prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> cos x \<and> cos x \<le> Ifloat u" 

1155 
proof (rule allI, rule allI, rule allI, rule impI) 

1156 
fix x lx ux 

1157 
assume "(l, u) = bnds_cos prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}" 

1158 
hence bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto 

1159 

1160 
let ?lpi = "lb_pi prec" 

1161 
have [intro!]: "Ifloat lx \<le> Ifloat ux" using x by auto 

1162 
hence "lx \<le> ux" unfolding le_float_def . 

1163 

1164 
show "Ifloat l \<le> cos x \<and> cos x \<le> Ifloat u" 

1165 
proof (cases "lx < ?lpi \<or> ux > ?lpi") 

1166 
case True 

1167 
show ?thesis using bnds unfolding bnds_cos_def if_P[OF True] Let_def using cos_le_one cos_ge_minus_one by auto 

1168 
next 

1169 
case False note not_out = this 

1170 
hence lpi_lx: " Ifloat ?lpi \<le> Ifloat lx" and lpi_ux: "Ifloat ux \<le> Ifloat ?lpi" unfolding le_float_def less_float_def by auto 

1171 

1172 
from pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1, THEN le_imp_neg_le] lpi_lx 

1173 
have " pi \<le> Ifloat lx" by (rule order_trans) 

1174 
hence " pi \<le> x" and " pi \<le> Ifloat ux" and "x \<le> Ifloat ux" using x by auto 

1175 

1176 
from lpi_ux pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1] 

1177 
have "Ifloat ux \<le> pi" by (rule order_trans) 

1178 
hence "x \<le> pi" and "Ifloat lx \<le> pi" and "Ifloat lx \<le> x" using x by auto 

1179 

1180 
note lb_cos_minus_bottom = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct1] 

1181 
note lb_cos_minus_top = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct2] 

1182 
note lb_cos_bottom = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct1] 

1183 
note lb_cos_top = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct2] 

1184 

1185 
show ?thesis 

1186 
proof (cases "ux \<le> 0") 

1187 
case True hence "Ifloat ux \<le> 0" unfolding le_float_def by auto 

1188 
hence "x \<le> 0" and "Ifloat lx \<le> 0" using x by auto 

1189 

1190 
{ have "Ifloat (lb_cos prec (lx)) \<le> cos (Ifloat (lx))" using lb_cos_minus_bottom[OF `pi \<le> Ifloat lx` `Ifloat lx \<le> 0`] . 

1191 
also have "\<dots> \<le> cos x" unfolding Ifloat_minus cos_minus using cos_monotone_minus_pi_0'[OF ` pi \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> 0`] . 

1192 
finally have "Ifloat (lb_cos prec (lx)) \<le> cos x" . } 

1193 
moreover 

1194 
{ have "cos x \<le> cos (Ifloat (ux))" unfolding Ifloat_minus cos_minus using cos_monotone_minus_pi_0'[OF ` pi \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> 0`] . 

1195 
also have "\<dots> \<le> Ifloat (ub_cos prec (ux))" using lb_cos_minus_top[OF `pi \<le> Ifloat ux` `Ifloat ux \<le> 0`] . 

1196 
finally have "cos x \<le> Ifloat (ub_cos prec (ux))" . } 

1197 
ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_P[OF True] by auto 

1198 
next 

1199 
case False note not_ux = this 

1200 

1201 
show ?thesis 

1202 
proof (cases "0 \<le> lx") 

1203 
case True hence "0 \<le> Ifloat lx" unfolding le_float_def by auto 

1204 
hence "0 \<le> x" and "0 \<le> Ifloat ux" using x by auto 

1205 

1206 
{ have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] . 

1207 
also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] . 

1208 
finally have "Ifloat (lb_cos prec ux) \<le> cos x" . } 

1209 
moreover 

1210 
{ have "cos x \<le> cos (Ifloat lx)" using cos_monotone_0_pi'[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> pi`] . 

1211 
also have "\<dots> \<le> Ifloat (ub_cos prec lx)" using lb_cos_top[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> pi`] . 

1212 
finally have "cos x \<le> Ifloat (ub_cos prec lx)" . } 

1213 
ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_not_P[OF not_ux] if_P[OF True] by auto 

1214 
next 

1215 
case False with not_ux 

1216 
have "Ifloat lx \<le> 0" and "0 \<le> Ifloat ux" unfolding le_float_def by auto 

1217 

1218 
have "Ifloat (min (lb_cos prec (lx)) (lb_cos prec ux)) \<le> cos x" 

1219 
proof (cases "x \<le> 0") 

1220 
case True 

1221 
have "Ifloat (lb_cos prec (lx)) \<le> cos (Ifloat (lx))" using lb_cos_minus_bottom[OF `pi \<le> Ifloat lx` `Ifloat lx \<le> 0`] . 

1222 
also have "\<dots> \<le> cos x" unfolding Ifloat_minus cos_minus using cos_monotone_minus_pi_0'[OF ` pi \<le> Ifloat lx` `Ifloat lx \<le> x` `x \<le> 0`] . 

1223 
finally show ?thesis unfolding Ifloat_min by auto 

1224 
next 

1225 
case False hence "0 \<le> x" by auto 

1226 
have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] . 

1227 
also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] . 

1228 
finally show ?thesis unfolding Ifloat_min by auto 

1229 
qed 

1230 
moreover have "cos x \<le> Ifloat (Float 1 0)" by auto 

1231 
ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_not_P[OF not_ux] if_not_P[OF False] by auto 

1232 
qed 

1233 
qed 

1234 
qed 

1235 
qed 

1236 

1237 
subsection "Compute the sinus in the entire domain" 

1238 

1239 
function lb_sin :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_sin :: "nat \<Rightarrow> float \<Rightarrow> float" where 

1240 
"lb_sin prec x = (let sqr_diff = \<lambda> x. if x > 1 then 0 else 1  x * x 

1241 
in if x < 0 then  ub_sin prec ( x) 

1242 
else if x \<le> Float 1 1 then x * lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 2 1 (x * x) 

1243 
else the (lb_sqrt prec (sqr_diff (ub_cos prec x))))"  

1244 

1245 
"ub_sin prec x = (let sqr_diff = \<lambda> x. if x < 0 then 1 else 1  x * x 

1246 
in if x < 0 then  lb_sin prec ( x) 

1247 
else if x \<le> Float 1 1 then x * ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 2 1 (x * x) 

1248 
else the (ub_sqrt prec (sqr_diff (lb_cos prec x))))" 

1249 
by pat_completeness auto 

1250 
termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def) 

1251 

1252 
definition bnds_sin :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where 

1253 
"bnds_sin prec lx ux = (let 

1254 
lpi = lb_pi prec ; 

1255 
half_pi = lpi * Float 1 1 

1256 
in if lx \<le>  half_pi \<or> half_pi \<le> ux then (Float 1 0, Float 1 0) 

1257 
else (lb_sin prec lx, ub_sin prec ux))" 

1258 

1259 
lemma lb_sin: assumes " (pi / 2) \<le> Ifloat x" and "Ifloat x \<le> pi / 2" 

1260 
shows "sin (Ifloat x) \<in> { Ifloat (lb_sin prec x) .. Ifloat (ub_sin prec x) }" (is "?sin x \<in> { ?lb x .. ?ub x}") 

1261 
proof  

1262 
{ fix x :: float assume "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2" 

1263 
hence "\<not> (x < 0)" and " (pi / 2) \<le> Ifloat x" unfolding less_float_def using pi_ge_two by auto 

1264 

1265 
have "Ifloat x \<le> pi" using `Ifloat x \<le> pi / 2` using pi_ge_two by auto 

1266 

1267 
have "?sin x \<in> { ?lb x .. ?ub x}" 

1268 
proof (cases "x \<le> Float 1 1") 

1269 
case True from sin_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`] 

1270 
show ?thesis unfolding lb_sin.simps[of prec x] ub_sin.simps[of prec x] if_not_P[OF `\<not> (x < 0)`] if_P[OF True] Let_def . 

1271 
next 

1272 
case False 

1273 
have "0 \<le> cos (Ifloat x)" using cos_ge_zero[OF _ `Ifloat x \<le> pi /2`] `0 \<le> Ifloat x` pi_ge_two by auto 

1274 
have "0 \<le> sin (Ifloat x)" using `0 \<le> Ifloat x` and `Ifloat x \<le> pi / 2` using sin_ge_zero by auto 

1275 

1276 
have "?sin x \<le> ?ub x" 

1277 
proof (cases "lb_cos prec x < 0") 

1278 
case True 

1279 
have "?sin x \<le> 1" using sin_le_one . 

1280 
also have "\<dots> \<le> Ifloat (the (ub_sqrt prec 1))" using ub_sqrt_lower_bound[where prec=prec and x=1] unfolding Ifloat_1 by auto 

1281 
finally show ?thesis unfolding ub_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_P[OF True] Let_def . 

1282 
next 

1283 
case False hence "0 \<le> Ifloat (lb_cos prec x)" unfolding less_float_def by auto 

1284 

1285 
have "sin (Ifloat x) = sqrt (1  cos (Ifloat x) ^ 2)" unfolding sin_squared_eq[symmetric] real_sqrt_abs using `0 \<le> sin (Ifloat x)` by auto 

1286 
also have "\<dots> \<le> sqrt (Ifloat (1  lb_cos prec x * lb_cos prec x))" 

1287 
proof (rule real_sqrt_le_mono) 

30273
ecd6f0ca62ea
declare power_Suc [simp]; remove redundant typespecific versions of power_Suc
huffman
parents:
30122
diff
changeset

1288 
have "Ifloat (lb_cos prec x * lb_cos prec x) \<le> cos (Ifloat x) ^ 2" unfolding numeral_2_eq_2 power_Suc2 power_0 Ifloat_mult 
29805  1289 
using `0 \<le> Ifloat (lb_cos prec x)` lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] `0 \<le> cos (Ifloat x)` by(auto intro!: mult_mono) 
1290 
thus "1  cos (Ifloat x) ^ 2 \<le> Ifloat (1  lb_cos prec x * lb_cos prec x)" unfolding Ifloat_sub Ifloat_1 by auto 

1291 
qed 

1292 
also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1  lb_cos prec x * lb_cos prec x)))" 

1293 
proof (rule ub_sqrt_lower_bound) 

1294 
have "Ifloat (lb_cos prec x) \<le> cos (Ifloat x)" using lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] by auto 

1295 
from mult_mono[OF order_trans[OF this cos_le_one] order_trans[OF this cos_le_one]] 

1296 
have "Ifloat (lb_cos prec x) * Ifloat (lb_cos prec x) \<le> 1" using `0 \<le> Ifloat (lb_cos prec x)` by auto 

1297 
thus "0 \<le> Ifloat (1  lb_cos prec x * lb_cos prec x)" by auto 

1298 
qed 

1299 
finally show ?thesis unfolding ub_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False] Let_def . 

1300 
qed 

1301 
moreover 

1302 
have "?lb x \<le> ?sin x" 

1303 
proof (cases "1 < ub_cos prec x") 

1304 
case True 

1305 
show ?thesis unfolding lb_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_P[OF True] Let_def 

1306 
by (rule order_trans[OF _ sin_ge_zero[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`]]) 

1307 
(auto simp add: lb_sqrt_upper_bound[where prec=prec and x=0, unfolded Ifloat_0 real_sqrt_zero]) 

1308 
next 

1309 
case False hence "Ifloat (ub_cos prec x) \<le> 1" unfolding less_float_def by auto 

1310 
have "0 \<le> Ifloat (ub_cos prec x)" using order_trans[OF `0 \<le> cos (Ifloat x)`] lb_cos `0 \<le> Ifloat x` `Ifloat x \<le> pi` by auto 

1311 

1312 
have "Ifloat (the (lb_sqrt prec (1  ub_cos prec x * ub_cos prec x))) \<le> sqrt (Ifloat (1  ub_cos prec x * ub_cos prec x))" 

1313 
proof (rule lb_sqrt_upper_bound) 

1314 
from mult_mono[OF `Ifloat (ub_cos prec x) \<le> 1` `Ifloat (ub_cos prec x) \<le> 1`] `0 \<le> Ifloat (ub_cos prec x)` 

1315 
have "Ifloat (ub_cos prec x) * Ifloat (ub_cos prec x) \<le> 1" by auto 

1316 
thus "0 \<le> Ifloat (1  ub_cos prec x * ub_cos prec x)" by auto 

1317 
qed 

1318 
also have "\<dots> \<le> sqrt (1  cos (Ifloat x) ^ 2)" 

1319 
proof (rule real_sqrt_le_mono) 

30273
ecd6f0ca62ea
declare power_Suc [simp]; remove redundant typespecific versions of power_Suc
huffman
parents:
30122
diff
changeset

1320 
have "cos (Ifloat x) ^ 2 \<le> Ifloat (ub_cos prec x * ub_cos prec x)" unfolding numeral_2_eq_2 power_Suc2 power_0 Ifloat_mult 
29805  1321 
using `0 \<le> Ifloat (ub_cos prec x)` lb_cos[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`] `0 \<le> cos (Ifloat x)` by(auto intro!: mult_mono) 
1322 
thus "Ifloat (1  ub_cos prec x * ub_cos prec x) \<le> 1  cos (Ifloat x) ^ 2" unfolding Ifloat_sub Ifloat_1 by auto 

1323 
qed 

1324 
also have "\<dots> = sin (Ifloat x)" unfolding sin_squared_eq[symmetric] real_sqrt_abs using `0 \<le> sin (Ifloat x)` by auto 

1325 
finally show ?thesis unfolding lb_sin.simps if_not_P[OF `\<not> (x < 0)`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False] Let_def . 

1326 
qed 

1327 
ultimately show ?thesis by auto 

1328 
qed 

1329 
} note for_pos = this 

1330 

1331 
show ?thesis 

1332 
proof (cases "x < 0") 

1333 
case True 

1334 
hence "0 \<le> Ifloat (x)" and "Ifloat ( x) \<le> pi / 2" using `(pi/2) \<le> Ifloat x` unfolding less_float_def by auto 

1335 
from for_pos[OF this] 

1336 
show ?thesis unfolding Ifloat_minus sin_minus lb_sin.simps[of prec x] ub_sin.simps[of prec x] if_P[OF True] Let_def atLeastAtMost_iff by auto 

1337 
next 

1338 
case False hence "0 \<le> Ifloat x" unfolding less_float_def by auto 

1339 
from for_pos[OF this `Ifloat x \<le> pi /2`] 

1340 
show ?thesis . 

1341 
qed 

1342 
qed 

1343 

1344 
lemma bnds_sin: "\<forall> x lx ux. (l, u) = bnds_sin prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> sin x \<and> sin x \<le> Ifloat u" 

1345 
proof (rule allI, rule allI, rule allI, rule impI) 

1346 
fix x lx ux 

1347 
assume "(l, u) = bnds_sin prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}" 

1348 
hence bnds: "(l, u) = bnds_sin prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto 

1349 
show "Ifloat l \<le> sin x \<and> sin x \<le> Ifloat u" 

1350 
proof (cases "lx \<le>  (lb_pi prec * Float 1 1) \<or> lb_pi prec * Float 1 1 \<le> ux") 

1351 
case True show ?thesis using bnds unfolding bnds_sin_def if_P[OF True] Let_def by auto 

1352 
next 

1353 
case False 

1354 
hence " lb_pi prec * Float 1 1 \<le> lx" and "ux \<le> lb_pi prec * Float 1 1" unfolding le_float_def by auto 

1355 
moreover have "Ifloat (lb_pi prec * Float 1 1) \<le> pi / 2" unfolding Ifloat_mult using pi_boundaries by auto 

1356 
ultimately have " (pi / 2) \<le> Ifloat lx" and "Ifloat ux \<le> pi / 2" and "Ifloat lx \<le> Ifloat ux" unfolding le_float_def using x by auto 

1357 
hence " (pi / 2) \<le> Ifloat ux" and "Ifloat lx \<le> pi / 2" by auto 

1358 

1359 
have " (pi / 2) \<le> x""x \<le> pi / 2" using `Ifloat ux \<le> pi / 2` ` (pi /2) \<le> Ifloat lx` x by auto 

1360 

1361 
{ have "Ifloat (lb_sin prec lx) \<le> sin (Ifloat lx)" using lb_sin[OF ` (pi / 2) \<le> Ifloat lx` `Ifloat lx \<le> pi / 2`] unfolding atLeastAtMost_iff by auto 

1362 
also have "\<dots> \<le> sin x" using sin_monotone_2pi' ` (pi / 2) \<le> Ifloat lx` x `x \<le> pi / 2` by auto 

1363 
finally have "Ifloat (lb_sin prec lx) \<le> sin x" . } 

1364 
moreover 

1365 
{ have "sin x \<le> sin (Ifloat ux)" using sin_monotone_2pi' ` (pi / 2) \<le> x` x `Ifloat ux \<le> pi / 2` by auto 

1366 
also have "\<dots> \<le> Ifloat (ub_sin prec ux)" using lb_sin[OF ` (pi / 2) \<le> Ifloat ux` `Ifloat ux \<le> pi / 2`] unfolding atLeastAtMost_iff by auto 

1367 
finally have "sin x \<le> Ifloat (ub_sin prec ux)" . } 

1368 
ultimately 

1369 
show ?thesis using bnds unfolding bnds_sin_def if_not_P[OF False] Let_def by auto 

1370 
qed 

1371 
qed 

1372 

1373 
section "Exponential function" 

1374 

1375 
subsection "Compute the series of the exponential function" 

1376 

1377 
fun ub_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and lb_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where 

1378 
"ub_exp_horner prec 0 i k x = 0"  

1379 
"ub_exp_horner prec (Suc n) i k x = rapprox_rat prec 1 (int k) + x * lb_exp_horner prec n (i + 1) (k * i) x"  

1380 
"lb_exp_horner prec 0 i k x = 0"  

1381 
"lb_exp_horner prec (Suc n) i k x = lapprox_rat prec 1 (int k) + x * ub_exp_horner prec n (i + 1) (k * i) x" 

1382 

1383 
lemma bnds_exp_horner: assumes "Ifloat x \<le> 0" 

1384 
shows "exp (Ifloat x) \<in> { Ifloat (lb_exp_horner prec (get_even n) 1 1 x) .. Ifloat (ub_exp_horner prec (get_odd n) 1 1 x) }" 

1385 
proof  

1386 
{ fix n 

1387 
have F: "\<And> m. ((\<lambda>i. i + 1) ^ n) m = n + m" by (induct n, auto) 

1388 
have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^ n) 1" unfolding F by auto } note f_eq = this 

1389 

1390 
note bounds = horner_bounds_nonpos[where f="fact" and lb="lb_exp_horner prec" and ub="ub_exp_horner prec" and j'=0 and s=1, 

1391 
OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps] 

1392 

1393 
{ have "Ifloat (lb_exp_horner prec (get_even n) 1 1 x) \<le> (\<Sum>j = 0..<get_even n. 1 / real (fact j) * Ifloat x ^ j)" 

1394 
using bounds(1) by auto 

1395 
also have "\<dots> \<le> exp (Ifloat x)" 

1396 
proof  

1397 
obtain t where "\<bar>t\<bar> \<le> \<bar>Ifloat x\<bar>" and "exp (Ifloat x) = (\<Sum>m = 0..<get_even n. (Ifloat x) ^ m / real (fact m)) + exp t / real (fact (get_even n)) * (Ifloat x) ^ (get_even n)" 

1398 
using Maclaurin_exp_le by blast 

1399 
moreover have "0 \<le> exp t / real (fact (get_even n)) * (Ifloat x) ^ (get_even n)" 

1400 
by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero) 

1401 
ultimately show ?thesis 

1402 
using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg) 

1403 
qed 

1404 
finally have "Ifloat (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (Ifloat x)" . 

1405 
} moreover 

1406 
{ 

1407 
have x_less_zero: "Ifloat x ^ get_odd n \<le> 0" 

1408 
proof (cases "Ifloat x = 0") 

1409 
case True 

1410 
have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto 

1411 
thus ?thesis unfolding True power_0_left by auto 

1412 
next 

1413 
case False hence "Ifloat x < 0" using `Ifloat x \<le> 0` by auto 

1414 
show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `Ifloat x < 0`) 

1415 
qed 

1416 

1417 
obtain t where "\<bar>t\<bar> \<le> \<bar>Ifloat x\<bar>" and "exp (Ifloat x) = (\<Sum>m = 0..<get_odd n. (Ifloat x) ^ m / real (fact m)) + exp t / real (fact (get_odd n)) * (Ifloat x) ^ (get_odd n)" 

1418 
using Maclaurin_exp_le by blast 

1419 
moreover have "exp t / real (fact (get_odd n)) * (Ifloat x) ^ (get_odd n) \<le> 0" 

1420 
by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero) 

1421 
ultimately have "exp (Ifloat x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * Ifloat x ^ j)" 

1422 
using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg) 

1423 
also have "\<dots> \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)" 

1424 
using bounds(2) by auto 

1425 
finally have "exp (Ifloat x) \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)" . 

1426 
} ultimately show ?thesis by auto 

1427 
qed 

1428 

1429 
subsection "Compute the exponential function on the entire domain" 

1430 

1431 
function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where 

1432 
"lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (x)) 

1433 
else let 

1434 
horner = (\<lambda> x. let y = lb_exp_horner prec (get_even (prec + 2)) 1 1 x in if y \<le> 0 then Float 1 2 else y) 

1435 
in if x <  1 then (case floor_fl x of (Float m e) \<Rightarrow> (horner (float_divl prec x ( Float m e))) ^ (nat (m) * 2 ^ nat e)) 

1436 
else horner x)"  

1437 
"ub_exp prec x = (if 0 < x then float_divr prec 1 (lb_exp prec (x)) 

a5da150bd0ab
Add approximation method
