src/HOL/Decision_Procs/Approximation.thy
changeset 54269 dcdfec41a325
parent 54230 b1d955791529
child 54489 03ff4d1e6784
equal deleted inserted replaced
54268:807532d15d16 54269:dcdfec41a325
   130   "get_even n = (if even n then n else (Suc n))"
   130   "get_even n = (if even n then n else (Suc n))"
   131 
   131 
   132 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
   132 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
   133 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
   133 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
   134 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
   134 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
   135 proof (cases "odd n")
   135   by (auto simp: get_odd_def odd_pos intro!: exI[of _ "n - 1"])
   136   case True hence "0 < n" by (rule odd_pos)
       
   137   from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto
       
   138   thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
       
   139 next
       
   140   case False hence "odd (Suc n)" by auto
       
   141   thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
       
   142 qed
       
   143 
   136 
   144 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
   137 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
   145 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
   138 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
   146 
   139 
   147 section "Power function"
   140 section "Power function"
   149 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
   142 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
   150 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
   143 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
   151                       else if u < 0         then (u ^ n, l ^ n)
   144                       else if u < 0         then (u ^ n, l ^ n)
   152                                             else (0, (max (-l) u) ^ n))"
   145                                             else (0, (max (-l) u) ^ n))"
   153 
   146 
   154 lemma float_power_bnds: fixes x :: real
   147 lemma float_power_bnds: "(l1, u1) = float_power_bnds n l u \<Longrightarrow> x \<in> {l .. u} \<Longrightarrow> (x::real) ^ n \<in> {l1..u1}"
   155   assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {l .. u}"
   148   by (auto simp: float_power_bnds_def max_def split: split_if_asm
   156   shows "x ^ n \<in> {l1..u1}"
   149            intro: power_mono_odd power_mono power_mono_even zero_le_even_power)
   157 proof (cases "even n")
       
   158   case True
       
   159   show ?thesis
       
   160   proof (cases "0 < l")
       
   161     case True hence "odd n \<or> 0 < l" and "0 \<le> real l" by auto
       
   162     have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms
       
   163       unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto
       
   164     have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using `0 \<le> real l` assms
       
   165       by (auto simp: power_mono)
       
   166     thus ?thesis using assms `0 < l` unfolding l1 u1 by auto
       
   167   next
       
   168     case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
       
   169     show ?thesis
       
   170     proof (cases "u < 0")
       
   171       case True hence "0 \<le> - real u" and "- real u \<le> - x" and "0 \<le> - x" and "-x \<le> - real l" using assms  by auto
       
   172       hence "real u ^ n \<le> x ^ n" and "x ^ n \<le> real l ^ n" using power_mono[of  "-x" "-real l" n] power_mono[of "-real u" "-x" n]
       
   173         unfolding power_minus_even[OF `even n`] by auto
       
   174       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
       
   175       ultimately show ?thesis by auto
       
   176     next
       
   177       case False
       
   178       have "\<bar>x\<bar> \<le> real (max (-l) u)"
       
   179       proof (cases "-l \<le> u")
       
   180         case True thus ?thesis unfolding max_def if_P[OF True] using assms by auto
       
   181       next
       
   182         case False thus ?thesis unfolding max_def if_not_P[OF False] using assms by auto
       
   183       qed
       
   184       hence x_abs: "\<bar>x\<bar> \<le> \<bar>real (max (-l) u)\<bar>" by auto
       
   185       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
       
   186       show ?thesis unfolding atLeastAtMost_iff l1 u1 using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto
       
   187     qed
       
   188   qed
       
   189 next
       
   190   case False hence "odd n \<or> 0 < l" by auto
       
   191   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
       
   192   have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using assms unfolding atLeastAtMost_iff using power_mono_odd[OF False] by auto
       
   193   thus ?thesis unfolding atLeastAtMost_iff l1 u1 less_float_def by auto
       
   194 qed
       
   195 
   150 
   196 lemma bnds_power: "\<forall> (x::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1"
   151 lemma bnds_power: "\<forall> (x::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1"
   197   using float_power_bnds by auto
   152   using float_power_bnds by auto
   198 
   153 
   199 section "Square root"
   154 section "Square root"
   242     by (simp add: field_simps power2_eq_square)
   197     by (simp add: field_simps power2_eq_square)
   243   thus ?thesis by (simp add: field_simps)
   198   thus ?thesis by (simp add: field_simps)
   244 qed
   199 qed
   245 
   200 
   246 lemma sqrt_iteration_bound: assumes "0 < real x"
   201 lemma sqrt_iteration_bound: assumes "0 < real x"
   247   shows "sqrt x < (sqrt_iteration prec n x)"
   202   shows "sqrt x < sqrt_iteration prec n x"
   248 proof (induct n)
   203 proof (induct n)
   249   case 0
   204   case 0
   250   show ?case
   205   show ?case
   251   proof (cases x)
   206   proof (cases x)
   252     case (Float m e)
   207     case (Float m e)
   354     hence "sqrt x \<le> ub_sqrt prec x"
   309     hence "sqrt x \<le> ub_sqrt prec x"
   355       unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }
   310       unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }
   356   note ub = this
   311   note ub = this
   357 
   312 
   358   show ?thesis
   313   show ?thesis
   359   proof (cases "0 < x")
   314     using lb[of "-x"] ub[of "-x"] lb[of x] ub[of x]
   360     case True with lb ub show ?thesis by auto
   315     by (auto simp add: lb_sqrt.simps ub_sqrt.simps real_sqrt_minus)
   361   next case False show ?thesis
       
   362   proof (cases "real x = 0")
       
   363     case True thus ?thesis
       
   364       by (auto simp add: lb_sqrt.simps ub_sqrt.simps)
       
   365   next
       
   366     case False with `\<not> 0 < x` have "x < 0" and "0 < -x"
       
   367       by auto
       
   368 
       
   369     with `\<not> 0 < x`
       
   370     show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`]
       
   371       by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps)
       
   372   qed qed
       
   373 qed
   316 qed
   374 
   317 
   375 lemma bnds_sqrt: "\<forall> (x::real) lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> sqrt x \<and> sqrt x \<le> u"
   318 lemma bnds_sqrt: "\<forall> (x::real) lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {lx .. ux} \<longrightarrow> l \<le> sqrt x \<and> sqrt x \<le> u"
   376 proof ((rule allI) +, rule impI, erule conjE, rule conjI)
   319 proof ((rule allI) +, rule impI, erule conjE, rule conjI)
   377   fix x :: real fix lx ux
   320   fix x :: real fix lx ux
   410 
   353 
   411 lemma arctan_0_1_bounds':
   354 lemma arctan_0_1_bounds':
   412   assumes "0 \<le> real x" "real x \<le> 1" and "even n"
   355   assumes "0 \<le> real x" "real x \<le> 1" and "even n"
   413   shows "arctan x \<in> {(x * lb_arctan_horner prec n 1 (x * x)) .. (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"
   356   shows "arctan x \<in> {(x * lb_arctan_horner prec n 1 (x * x)) .. (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"
   414 proof -
   357 proof -
   415   let "?c i" = "-1^i * (1 / (i * 2 + (1::nat)) * real x ^ (i * 2 + 1))"
   358   let ?c = "\<lambda>i. -1^i * (1 / (i * 2 + (1::nat)) * real x ^ (i * 2 + 1))"
   416   let "?S n" = "\<Sum> i=0..<n. ?c i"
   359   let ?S = "\<lambda>n. \<Sum> i=0..<n. ?c i"
   417 
   360 
   418   have "0 \<le> real (x * x)" by auto
   361   have "0 \<le> real (x * x)" by auto
   419   from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
   362   from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
   420 
   363 
   421   have "arctan x \<in> { ?S n .. ?S (Suc n) }"
   364   have "arctan x \<in> { ?S n .. ?S (Suc n) }"
   455   ultimately show ?thesis by auto
   398   ultimately show ?thesis by auto
   456 qed
   399 qed
   457 
   400 
   458 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"
   401 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"
   459   shows "arctan x \<in> {(x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"
   402   shows "arctan x \<in> {(x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"
   460 proof (cases "even n")
   403   using
   461   case True
   404     arctan_0_1_bounds'[OF assms, of n prec]
   462   obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
   405     arctan_0_1_bounds'[OF assms, of "n + 1" prec]
   463   hence "even n'" unfolding even_Suc by auto
   406     arctan_0_1_bounds'[OF assms, of "n - 1" prec]
   464   have "arctan x \<le> x * ub_arctan_horner prec (get_odd n) 1 (x * x)"
   407   by (auto simp: get_even_def get_odd_def odd_pos simp del: ub_arctan_horner.simps lb_arctan_horner.simps)
   465     unfolding `Suc n' = get_odd n`[symmetric] using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto
       
   466   moreover
       
   467   have "x * lb_arctan_horner prec (get_even n) 1 (x * x) \<le> arctan x"
       
   468     unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n`] by auto
       
   469   ultimately show ?thesis by auto
       
   470 next
       
   471   case False hence "0 < n" by (rule odd_pos)
       
   472   from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
       
   473   from False[unfolded this even_Suc]
       
   474   have "even n'" and "even (Suc (Suc n'))" by auto
       
   475   have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
       
   476 
       
   477   have "arctan x \<le> x * ub_arctan_horner prec (get_odd n) 1 (x * x)"
       
   478     unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto
       
   479   moreover
       
   480   have "(x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan x"
       
   481     unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even (Suc (Suc n'))`] by auto
       
   482   ultimately show ?thesis by auto
       
   483 qed
       
   484 
   408 
   485 subsection "Compute \<pi>"
   409 subsection "Compute \<pi>"
   486 
   410 
   487 definition ub_pi :: "nat \<Rightarrow> float" where
   411 definition ub_pi :: "nat \<Rightarrow> float" where
   488   "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
   412   "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
   528       using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
   452       using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
   529     also have "\<dots> \<le> arctan (1 / k)" using `?k \<le> 1 / k` by (rule arctan_monotone')
   453     also have "\<dots> \<le> arctan (1 / k)" using `?k \<le> 1 / k` by (rule arctan_monotone')
   530     finally have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan (1 / k)" .
   454     finally have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan (1 / k)" .
   531   } note lb_arctan = this
   455   } note lb_arctan = this
   532 
   456 
   533   have "pi \<le> ub_pi n"
   457   have "pi \<le> ub_pi n \<and> lb_pi n \<le> pi"
   534     unfolding ub_pi_def machin_pi Let_def unfolding Float_num
   458     unfolding lb_pi_def ub_pi_def machin_pi Let_def unfolding Float_num
   535     using lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2]
   459     using lb_arctan[of 5] ub_arctan[of 239] lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2]
   536     by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff)
   460     by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff)
   537   moreover
   461   then show ?thesis by auto
   538   have "lb_pi n \<le> pi"
       
   539     unfolding lb_pi_def machin_pi Let_def Float_num
       
   540     using lb_arctan[of 5] ub_arctan[of 239] powr_realpow[of 2 2]
       
   541     by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff)
       
   542   ultimately show ?thesis by auto
       
   543 qed
   462 qed
   544 
   463 
   545 subsection "Compute arcus tangens in the entire domain"
   464 subsection "Compute arcus tangens in the entire domain"
   546 
   465 
   547 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
   466 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
  1806   using powr_gt_zero[of 2 "exponent x"]
  1725   using powr_gt_zero[of 2 "exponent x"]
  1807   apply simp
  1726   apply simp
  1808   done
  1727   done
  1809 
  1728 
  1810 lemma Float_pos_eq_mantissa_pos:  "Float m e > 0 \<longleftrightarrow> m > 0"
  1729 lemma Float_pos_eq_mantissa_pos:  "Float m e > 0 \<longleftrightarrow> m > 0"
  1811   apply (auto simp add: zero_less_mult_iff zero_float_def powr_gt_zero[of 2 "exponent x"] dest: less_zeroE)
       
  1812   using powr_gt_zero[of 2 "e"]
  1730   using powr_gt_zero[of 2 "e"]
  1813   apply simp
  1731   by (auto simp add: zero_less_mult_iff zero_float_def simp del: powr_gt_zero dest: less_zeroE)
  1814   done
       
  1815 
  1732 
  1816 lemma Float_representation_aux:
  1733 lemma Float_representation_aux:
  1817   fixes m e
  1734   fixes m e
  1818   defines "x \<equiv> Float m e"
  1735   defines "x \<equiv> Float m e"
  1819   assumes "x > 0"
  1736   assumes "x > 0"