src/HOL/Decision_Procs/Approximation.thy
author hoelzl
Wed, 11 Mar 2009 10:58:18 +0100
changeset 30439 57c68b3af2ea
parent 30413 c41afa5607be
child 30443 873fa77be5f0
permissions -rw-r--r--
Updated paths in Decision_Procs comments and NEWS
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
30439
57c68b3af2ea Updated paths in Decision_Procs comments and NEWS
hoelzl
parents: 30413
diff changeset
     1
(*  Title:      HOL/Decision_Procs/Approximation.thy
30122
1c912a9d8200 standard headers;
wenzelm
parents: 29823
diff changeset
     2
    Author:     Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009
1c912a9d8200 standard headers;
wenzelm
parents: 29823
diff changeset
     3
*)
1c912a9d8200 standard headers;
wenzelm
parents: 29823
diff changeset
     4
29805
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
     5
header {* Prove unequations about real numbers by computation *}
30122
1c912a9d8200 standard headers;
wenzelm
parents: 29823
diff changeset
     6
29805
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
     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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
     9
begin
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    10
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    11
section "Horner Scheme"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    12
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    13
subsection {* Define auxiliary helper @{text horner} function *}
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    14
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    15
fun horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    16
"horner F G 0 i k x       = 0" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    17
"horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    18
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    19
lemma horner_schema': fixes x :: real  and a :: "nat \<Rightarrow> real"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    21
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    23
  show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    24
    setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n  *a n * x^n"] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    25
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    26
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    27
lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    28
  assumes f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    30
proof (induct n arbitrary: i k j')
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    31
  case (Suc n)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    32
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    33
  show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    34
    using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    35
qed auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    36
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    37
lemma horner_bounds':
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    38
  assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    39
  and lb_0: "\<And> i k x. lb 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    41
  and ub_0: "\<And> i k x. ub 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    43
  shows "Ifloat (lb n ((F^j') s) (f j') x) \<le> horner F G n ((F^j') s) (f j') (Ifloat x) \<and> 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    44
         horner F G n ((F^j') s) (f j') (Ifloat x) \<le> Ifloat (ub n ((F^j') s) (f j') x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    45
  (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    46
proof (induct n arbitrary: j')
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    47
  case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    48
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    49
  case (Suc n)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    50
  have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps Ifloat_sub diff_def
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    51
  proof (rule add_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    53
    from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    55
      unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    56
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    57
  moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps Ifloat_sub diff_def
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    58
  proof (rule add_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    60
    from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    61
    show "- (Ifloat x * horner F G n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) (Ifloat x)) \<le> 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    62
          - Ifloat (x * lb n (F ((F ^ j') s)) (G ((F ^ j') s) (f j')) x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    63
      unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    64
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    65
  ultimately show ?case by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    66
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    67
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    68
subsection "Theorems for floating point functions implementing the horner scheme"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    69
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    70
text {*
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    71
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    72
Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    73
all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    74
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    75
*}
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    76
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    77
lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    78
  assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    79
  and lb_0: "\<And> i k x. lb 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    81
  and ub_0: "\<And> i k x. ub 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    85
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    86
  have "?lb  \<and> ?ub" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    87
    using horner_bounds'[where lb=lb, OF `0 \<le> Ifloat x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    88
    unfolding horner_schema[where f=f, OF f_Suc] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    89
  thus "?lb" and "?ub" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    90
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    91
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    92
lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    93
  assumes "Ifloat x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F^n) s) (f n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    94
  and lb_0: "\<And> i k x. lb 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    96
  and ub_0: "\<And> i k x. ub 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
    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")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   100
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   101
  { fix x y z :: float have "x - y * z = x + - y * z"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   102
      by (cases x, cases y, cases z, simp add: plus_float.simps minus_float.simps uminus_float.simps times_float.simps algebra_simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   103
  } note diff_mult_minus = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   104
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   105
  { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   106
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   107
  have move_minus: "Ifloat (-x) = -1 * Ifloat x" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   108
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   109
  have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x)^j) = 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   110
    (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   111
  proof (rule setsum_cong, simp)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   112
    fix j assume "j \<in> {0 ..< n}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   113
    show "1 / real (f (j' + j)) * Ifloat x ^ j = -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   114
      unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   115
      unfolding real_mult_commute unfolding real_mult_assoc[of "-1^j", symmetric] power_mult_distrib[symmetric]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   116
      by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   117
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   118
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   119
  have "0 \<le> Ifloat (-x)" using assms by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   120
  from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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,
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   122
    OF this f_Suc lb_0 refl ub_0 refl]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   123
  show "?lb" and "?ub" unfolding minus_minus sum_eq
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   124
    by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   125
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   126
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   127
subsection {* Selectors for next even or odd number *}
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   128
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   129
text {*
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   130
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   131
The horner scheme computes alternating series. To get the upper and lower bounds we need to
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   132
guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   133
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   134
*}
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   135
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   136
definition get_odd :: "nat \<Rightarrow> nat" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   137
  "get_odd n = (if odd n then n else (Suc n))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   138
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   139
definition get_even :: "nat \<Rightarrow> nat" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   140
  "get_even n = (if even n then n else (Suc n))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   141
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   142
lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   143
lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   144
lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   145
proof (cases "odd n")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   146
  case True hence "0 < n" by (rule odd_pos)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   147
  from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   148
  thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   149
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   150
  case False hence "odd (Suc n)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   151
  thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   152
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   153
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   154
lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   155
lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   156
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   157
section "Power function"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   158
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   159
definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   160
"float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   161
                      else if u < 0         then (u ^ n, l ^ n)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   162
                                            else (0, (max (-l) u) ^ n))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   163
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   164
lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {Ifloat l .. Ifloat u}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   165
  shows "x^n \<in> {Ifloat l1..Ifloat u1}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   166
proof (cases "even n")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   167
  case True 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   168
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   169
  proof (cases "0 < l")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   170
    case True hence "odd n \<or> 0 < l" and "0 \<le> Ifloat l" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   173
    thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   174
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   175
    case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   176
    show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   177
    proof (cases "u < 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   180
	unfolding power_minus_even[OF `even n`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   182
      ultimately show ?thesis using float_power by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   183
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   184
      case False 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   185
      have "\<bar>x\<bar> \<le> Ifloat (max (-l) u)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   186
      proof (cases "-l \<le> u")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   187
	case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   188
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   189
	case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   190
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   191
      hence x_abs: "\<bar>x\<bar> \<le> \<bar>Ifloat (max (-l) u)\<bar>" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   194
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   195
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   196
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   197
  case False hence "odd n \<or> 0 < l" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   200
  thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   201
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   202
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   204
  using float_power_bnds by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   205
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   206
section "Square root"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   207
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   208
text {*
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   209
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   210
The square root computation is implemented as newton iteration. As first first step we use the
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   211
nearest power of two greater than the square root.
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   212
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   213
*}
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   214
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   215
fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   216
"sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   217
"sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   218
                                  in Float 1 -1 * (y + float_divr prec x y))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   219
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   220
definition ub_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   221
"ub_sqrt prec x = (if 0 < x then Some (sqrt_iteration prec prec x) else if x < 0 then None else Some 0)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   222
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   223
definition lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   225
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   226
lemma sqrt_ub_pos_pos_1:
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   227
  assumes "sqrt x < b" and "0 < b" and "0 < x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   228
  shows "sqrt x < (b + x / b)/2"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   229
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   230
  from assms have "0 < (b - sqrt x) ^ 2 " by simp
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   231
  also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   232
  also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   233
  finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   234
  hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   235
    by (simp add: field_simps power2_eq_square)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   236
  thus ?thesis by (simp add: field_simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   237
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   238
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   239
lemma sqrt_iteration_bound: assumes "0 < Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   240
  shows "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec n x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   241
proof (induct n)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   242
  case 0
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   243
  show ?case
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   244
  proof (cases x)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   245
    case (Float m e)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   246
    hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   247
    hence "0 < sqrt (real m)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   248
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   249
    have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   250
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   251
    have "Ifloat x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   252
      unfolding pow2_add pow2_int Float Ifloat.simps by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   253
    also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   254
    proof (rule mult_strict_right_mono, auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   255
      show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2] 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   256
	unfolding real_of_int_less_iff[of m, symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   257
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   258
    finally have "sqrt (Ifloat x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   259
    also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   260
    proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   261
      let ?E = "e + bitlen m"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   262
      have E_mod_pow: "pow2 (?E mod 2) < 4"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   263
      proof (cases "?E mod 2 = 1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   264
	case True thus ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   265
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   266
	case False 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   267
	have "0 \<le> ?E mod 2" by auto 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   268
	have "?E mod 2 < 2" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   269
	from this[THEN zless_imp_add1_zle]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   270
	have "?E mod 2 \<le> 0" using False by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   271
	from xt1(5)[OF `0 \<le> ?E mod 2` this]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   272
	show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   273
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   274
      hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   275
      hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   276
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   277
      have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   278
      have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   279
	unfolding E_eq unfolding pow2_add ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   280
      also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   281
	unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   282
      also have "\<dots> < pow2 (?E div 2) * 2" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   283
	by (rule mult_strict_left_mono, auto intro: E_mod_pow)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   284
      also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   285
      finally show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   286
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   287
    finally show ?thesis 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   288
      unfolding Float sqrt_iteration.simps Ifloat.simps by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   289
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   290
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   291
  case (Suc n)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   292
  let ?b = "sqrt_iteration prec n x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   293
  have "0 < sqrt (Ifloat x)" using `0 < Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   294
  also have "\<dots> < Ifloat ?b" using Suc .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   296
  also have "\<dots> \<le> (Ifloat ?b + Ifloat (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   297
  also have "\<dots> = Ifloat (Float 1 -1) * (Ifloat ?b + Ifloat (float_divr prec x ?b))" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   298
  finally show ?case unfolding sqrt_iteration.simps Let_def Ifloat_mult Ifloat_add right_distrib .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   299
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   300
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   301
lemma sqrt_iteration_lower_bound: assumes "0 < Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   302
  shows "0 < Ifloat (sqrt_iteration prec n x)" (is "0 < ?sqrt")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   303
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   304
  have "0 < sqrt (Ifloat x)" using assms by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   305
  also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   306
  finally show ?thesis .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   307
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   308
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   309
lemma lb_sqrt_lower_bound: assumes "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   310
  shows "0 \<le> Ifloat (the (lb_sqrt prec x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   311
proof (cases "0 < x")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   313
  hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   315
  thus ?thesis unfolding lb_sqrt_def using True by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   316
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   317
  case False with `0 \<le> Ifloat x` have "Ifloat x = 0" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   318
  thus ?thesis unfolding lb_sqrt_def less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   319
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   320
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   321
lemma lb_sqrt_upper_bound: assumes "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   322
  shows "Ifloat (the (lb_sqrt prec x)) \<le> sqrt (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   323
proof (cases "0 < x")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   324
  case True hence "0 < Ifloat x" and "0 \<le> Ifloat x" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   325
  hence sqrt_gt0: "0 < sqrt (Ifloat x)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   326
  hence sqrt_ub: "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   327
  
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   329
  also have "\<dots> < Ifloat x / sqrt (Ifloat x)" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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]])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   332
  finally show ?thesis unfolding lb_sqrt_def if_P[OF `0 < x`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   333
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   334
  case False with `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   335
  have "\<not> x < 0" unfolding less_float_def le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   336
  show ?thesis unfolding lb_sqrt_def if_not_P[OF False] if_not_P[OF `\<not> x < 0`] using assms by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   337
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   338
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   339
lemma lb_sqrt: assumes "Some y = lb_sqrt prec x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   340
  shows "Ifloat y \<le> sqrt (Ifloat x)" and "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   341
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   342
  show "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   343
  proof (rule ccontr)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   344
    assume "\<not> 0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   345
    hence "lb_sqrt prec x = None" unfolding lb_sqrt_def less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   346
    thus False using assms by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   347
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   348
  from lb_sqrt_upper_bound[OF this, of prec]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   349
  show "Ifloat y \<le> sqrt (Ifloat x)" unfolding assms[symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   350
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   351
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   352
lemma ub_sqrt_lower_bound: assumes "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   353
  shows "sqrt (Ifloat x) \<le> Ifloat (the (ub_sqrt prec x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   354
proof (cases "0 < x")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   355
  case True hence "0 < Ifloat x" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   356
  hence "0 < sqrt (Ifloat x)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   357
  hence "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   358
  thus ?thesis unfolding ub_sqrt_def if_P[OF `0 < x`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   359
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   360
  case False with `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   361
  have "Ifloat x = 0" unfolding less_float_def le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   362
  thus ?thesis unfolding ub_sqrt_def less_float_def le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   363
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   364
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   365
lemma ub_sqrt: assumes "Some y = ub_sqrt prec x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   366
  shows "sqrt (Ifloat x) \<le> Ifloat y" and "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   367
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   368
  show "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   369
  proof (rule ccontr)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   370
    assume "\<not> 0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   371
    hence "ub_sqrt prec x = None" unfolding ub_sqrt_def less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   372
    thus False using assms by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   373
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   374
  from ub_sqrt_lower_bound[OF this, of prec]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   375
  show "sqrt (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   376
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   377
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   379
proof (rule allI, rule allI, rule allI, rule impI)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   380
  fix x lx ux
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   381
  assume "(Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   383
  
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   384
  have "Ifloat lx \<le> x" and "x \<le> Ifloat ux" using x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   385
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   386
  from lb_sqrt(1)[OF l] real_sqrt_le_mono[OF `Ifloat lx \<le> x`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   387
  have "Ifloat l \<le> sqrt x" by (rule order_trans)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   388
  moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   389
  from real_sqrt_le_mono[OF `x \<le> Ifloat ux`] ub_sqrt(1)[OF u]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   390
  have "sqrt x \<le> Ifloat u" by (rule order_trans)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   391
  ultimately show "Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u" ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   392
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   393
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   394
section "Arcus tangens and \<pi>"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   395
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   396
subsection "Compute arcus tangens series"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   397
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   398
text {*
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   399
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   400
As first step we implement the computation of the arcus tangens series. This is only valid in the range
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   401
@{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   402
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   403
*}
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   404
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   405
fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   406
and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   407
  "ub_arctan_horner prec 0 k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   408
| "ub_arctan_horner prec (Suc n) k x = 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   409
    (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   410
| "lb_arctan_horner prec 0 k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   411
| "lb_arctan_horner prec (Suc n) k x = 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   412
    (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   413
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   414
lemma arctan_0_1_bounds': assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1" and "even n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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))}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   416
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   417
  let "?c i" = "-1^i * (1 / real (i * 2 + 1) * Ifloat x ^ (i * 2 + 1))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   418
  let "?S n" = "\<Sum> i=0..<n. ?c i"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   419
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   420
  have "0 \<le> Ifloat (x * x)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   421
  from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   422
  
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   423
  have "arctan (Ifloat x) \<in> { ?S n .. ?S (Suc n) }"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   424
  proof (cases "Ifloat x = 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   425
    case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   426
    hence "0 < Ifloat x" using `0 \<le> Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   427
    hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * Ifloat x ^ (0 * 2 + 1)" by auto 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   428
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   429
    have "\<bar> Ifloat x \<bar> \<le> 1"  using `0 \<le> Ifloat x` `Ifloat x \<le> 1` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   431
    show ?thesis unfolding arctan_series[OF `\<bar> Ifloat x \<bar> \<le> 1`] Suc_plus1  .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   432
  qed auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   433
  note arctan_bounds = this[unfolded atLeastAtMost_iff]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   434
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   435
  have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   436
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   437
  note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   438
    and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   439
    and ub="\<lambda>n i k x. ub_arctan_horner prec n k x", 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   440
    OF `0 \<le> Ifloat (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   441
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   442
  { have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   443
      using bounds(1) `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   444
      unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   445
      unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   446
      by (auto intro!: mult_left_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   447
    also have "\<dots> \<le> arctan (Ifloat x)" using arctan_bounds ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   448
    finally have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (Ifloat x)" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   449
  moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   450
  { have "arctan (Ifloat x) \<le> ?S (Suc n)" using arctan_bounds ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   451
    also have "\<dots> \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   452
      using bounds(2)[of "Suc n"] `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   453
      unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   454
      unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   455
      by (auto intro!: mult_left_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   456
    finally have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   457
  ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   458
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   459
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   460
lemma arctan_0_1_bounds: assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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))}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   462
proof (cases "even n")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   463
  case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   464
  obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   465
  hence "even n'" unfolding even_nat_Suc by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   466
  have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   468
  moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   469
  have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   471
  ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   472
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   473
  case False hence "0 < n" by (rule odd_pos)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   474
  from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   475
  from False[unfolded this even_nat_Suc]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   476
  have "even n'" and "even (Suc (Suc n'))" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   477
  have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   478
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   479
  have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   481
  moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   482
  have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   484
  ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   485
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   486
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   487
subsection "Compute \<pi>"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   488
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   489
definition ub_pi :: "nat \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   490
  "ub_pi prec = (let A = rapprox_rat prec 1 5 ; 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   491
                     B = lapprox_rat prec 1 239
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   492
                 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) - 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   493
                                                  B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   494
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   495
definition lb_pi :: "nat \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   496
  "lb_pi prec = (let A = lapprox_rat prec 1 5 ; 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   497
                     B = rapprox_rat prec 1 239
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   498
                 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) - 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   499
                                                  B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   500
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   501
lemma pi_boundaries: "pi \<in> {Ifloat (lb_pi n) .. Ifloat (ub_pi n)}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   502
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   503
  have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   504
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   505
  { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   506
    let ?k = "rapprox_rat prec 1 k"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   507
    have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   508
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   509
    have "0 \<le> Ifloat ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   510
    have "Ifloat ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   511
      by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   512
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   513
    have "1 / real k \<le> Ifloat ?k" using rapprox_rat[where x=1 and y=k] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   514
    hence "arctan (1 / real k) \<le> arctan (Ifloat ?k)" by (rule arctan_monotone')
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   515
    also have "\<dots> \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   516
      using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   517
    finally have "arctan (1 / (real k)) \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   518
  } note ub_arctan = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   519
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   520
  { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   521
    let ?k = "lapprox_rat prec 1 k"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   522
    have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   523
    have "1 / real k \<le> 1" using `1 < k` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   524
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   526
    have "\<And>n. Ifloat ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   527
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   528
    have "Ifloat ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   529
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   530
    have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (Ifloat ?k)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   531
      using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   532
    also have "\<dots> \<le> arctan (1 / real k)" using `Ifloat ?k \<le> 1 / real k` by (rule arctan_monotone')
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   533
    finally have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   534
  } note lb_arctan = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   535
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   536
  have "pi \<le> Ifloat (ub_pi n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   537
    unfolding ub_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub unfolding Float_num
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   538
    using lb_arctan[of 239] ub_arctan[of 5]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   539
    by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   540
  moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   541
  have "Ifloat (lb_pi n) \<le> pi"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   542
    unfolding lb_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub Float_num
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   543
    using lb_arctan[of 5] ub_arctan[of 239]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   544
    by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   545
  ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   546
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   547
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   548
subsection "Compute arcus tangens in the entire domain"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   549
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   550
function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   551
  "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   552
                           lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   553
    in (if x < 0          then - ub_arctan prec (-x) else
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   554
        if x \<le> Float 1 -1 then lb_horner x else
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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))))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   556
                          else (let inv = float_divr prec 1 x 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   557
                                in if inv > 1 then 0 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   558
                                              else lb_pi prec * Float 1 -1 - ub_horner inv)))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   559
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   560
| "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   561
                           ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   562
    in (if x < 0          then - lb_arctan prec (-x) else
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   563
        if x \<le> Float 1 -1 then ub_horner x else
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   564
        if x \<le> Float 1 1  then let y = float_divr prec x (1 + the (lb_sqrt prec (1 + x * x)))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   565
                               in if y > 1 then ub_pi prec * Float 1 -1 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   566
                                           else Float 1 1 * ub_horner y 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   567
                          else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   568
by pat_completeness auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   570
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   571
declare ub_arctan_horner.simps[simp del]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   572
declare lb_arctan_horner.simps[simp del]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   573
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   574
lemma lb_arctan_bound': assumes "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   575
  shows "Ifloat (lb_arctan prec x) \<le> arctan (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   576
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   577
  have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   578
  let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   579
    and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   580
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   581
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   582
  proof (cases "x \<le> Float 1 -1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   583
    case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   584
    show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   585
      using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   586
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   587
    case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   588
    let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   589
    let ?fR = "1 + the (ub_sqrt prec (1 + x * x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   590
    let ?DIV = "float_divl prec x ?fR"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   591
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   593
    hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   594
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   596
    hence "?R \<le> Ifloat ?fR" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   597
    hence "0 < ?fR" and "0 < Ifloat ?fR" unfolding less_float_def using `0 < ?R` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   598
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   599
    have monotone: "Ifloat (float_divl prec x ?fR) \<le> Ifloat x / ?R"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   600
    proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   601
      have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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]])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   603
      finally show ?thesis .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   604
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   605
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   606
    show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   607
    proof (cases "x \<le> Float 1 1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   608
      case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   609
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   611
      also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   612
      finally have "Ifloat x \<le> Ifloat ?fR" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   613
      moreover have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   614
      ultimately have "Ifloat ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < Ifloat ?fR`, symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   615
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   616
      have "0 \<le> Ifloat ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   617
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   619
	using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   620
      also have "\<dots> \<le> 2 * arctan (Ifloat x / ?R)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   621
	using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
30273
ecd6f0ca62ea declare power_Suc [simp]; remove redundant type-specific 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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   624
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   625
      case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   626
      hence "2 < Ifloat x" unfolding le_float_def Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   627
      hence "1 \<le> Ifloat x" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   628
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   629
      let "?invx" = "float_divr prec 1 x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   631
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   632
      show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   633
      proof (cases "1 < ?invx")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   634
	case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   636
	  using `0 \<le> arctan (Ifloat x)` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   637
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   638
	case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   639
	hence "Ifloat ?invx \<le> 1" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   640
	have "0 \<le> Ifloat ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> Ifloat x`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   641
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   642
	have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   643
	
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   644
	have "arctan (1 / Ifloat x) \<le> arctan (Ifloat ?invx)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divr)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   646
	finally have "pi / 2 - Ifloat (?ub_horner ?invx) \<le> arctan (Ifloat x)" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   647
	  using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`] 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   648
	  unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   649
	moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   651
	ultimately
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   653
	  by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   654
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   655
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   656
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   657
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   658
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   659
lemma ub_arctan_bound': assumes "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   660
  shows "arctan (Ifloat x) \<le> Ifloat (ub_arctan prec x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   661
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   662
  have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   663
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   664
  let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   665
    and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   666
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   667
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   668
  proof (cases "x \<le> Float 1 -1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   669
    case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   670
    show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   671
      using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   672
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   673
    case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   674
    let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   675
    let ?fR = "1 + the (lb_sqrt prec (1 + x * x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   676
    let ?DIV = "float_divr prec x ?fR"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   677
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   679
    hence "0 \<le> Ifloat (1 + x*x)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   680
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   681
    hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   682
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   684
    hence "Ifloat ?fR \<le> ?R" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)`])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   686
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   687
    have monotone: "Ifloat x / ?R \<le> Ifloat (float_divr prec x ?fR)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   688
    proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   689
      from divide_left_mono[OF `Ifloat ?fR \<le> ?R` `0 \<le> Ifloat x` mult_pos_pos[OF divisor_gt0 `0 < Ifloat ?fR`]]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   690
      have "Ifloat x / ?R \<le> Ifloat x / Ifloat ?fR" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   691
      also have "\<dots> \<le> Ifloat ?DIV" by (rule float_divr)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   692
      finally show ?thesis .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   693
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   694
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   695
    show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   696
    proof (cases "x \<le> Float 1 1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   697
      case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   698
      show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   699
      proof (cases "?DIV > 1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   700
	case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   702
	from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   704
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   705
	case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   706
	hence "Ifloat ?DIV \<le> 1" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   707
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   708
	have "0 \<le> Ifloat x / ?R" using `0 \<le> Ifloat x` `0 < ?R` unfolding real_0_le_divide_iff by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   709
	hence "0 \<le> Ifloat ?DIV" using monotone by (rule order_trans)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   710
30273
ecd6f0ca62ea declare power_Suc [simp]; remove redundant type-specific 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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   712
	also have "\<dots> \<le> 2 * arctan (Ifloat ?DIV)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   713
	  using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   714
	also have "\<dots> \<le> Ifloat (Float 1 1 * ?ub_horner ?DIV)" unfolding Ifloat_mult[of "Float 1 1"] Float_num
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   715
	  using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   717
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   718
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   719
      case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   720
      hence "2 < Ifloat x" unfolding le_float_def Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   721
      hence "1 \<le> Ifloat x" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   722
      hence "0 < Ifloat x" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   723
      hence "0 < x" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   724
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   725
      let "?invx" = "float_divl prec 1 x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   727
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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`])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   730
	
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   731
      have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   732
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   734
      also have "\<dots> \<le> arctan (1 / Ifloat x)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divl)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   735
      finally have "arctan (Ifloat x) \<le> pi / 2 - Ifloat (?lb_horner ?invx)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   736
	using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`] 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   737
	unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   738
      moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   740
      ultimately
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   742
	by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   743
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   744
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   745
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   746
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   747
lemma arctan_boundaries:
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   748
  "arctan (Ifloat x) \<in> {Ifloat (lb_arctan prec x) .. Ifloat (ub_arctan prec x)}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   749
proof (cases "0 \<le> x")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   750
  case True hence "0 \<le> Ifloat x" unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   752
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   753
  let ?mx = "-x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   754
  case False hence "x < 0" and "0 \<le> Ifloat ?mx" unfolding le_float_def less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   755
  hence bounds: "Ifloat (lb_arctan prec ?mx) \<le> arctan (Ifloat ?mx) \<and> arctan (Ifloat ?mx) \<le> Ifloat (ub_arctan prec ?mx)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   756
    using ub_arctan_bound'[OF `0 \<le> Ifloat ?mx`] lb_arctan_bound'[OF `0 \<le> Ifloat ?mx`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   758
    unfolding atLeastAtMost_iff using bounds[unfolded Ifloat_minus arctan_minus] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   759
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   760
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   762
proof (rule allI, rule allI, rule allI, rule impI)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   763
  fix x lx ux
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   764
  assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   766
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   767
  { from arctan_boundaries[of lx prec, unfolded l]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   768
    have "Ifloat l \<le> arctan (Ifloat lx)" by (auto simp del: lb_arctan.simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   769
    also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   770
    finally have "Ifloat l \<le> arctan x" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   771
  } moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   772
  { have "arctan x \<le> arctan (Ifloat ux)" using x by (auto intro: arctan_monotone')
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   773
    also have "\<dots> \<le> Ifloat u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   774
    finally have "arctan x \<le> Ifloat u" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   775
  } ultimately show "Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u" ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   776
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   777
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   778
section "Sinus and Cosinus"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   779
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   780
subsection "Compute the cosinus and sinus series"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   781
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   782
fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   783
and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   784
  "ub_sin_cos_aux prec 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   785
| "ub_sin_cos_aux prec (Suc n) i k x = 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   786
    (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   787
| "lb_sin_cos_aux prec 0 i k x = 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   788
| "lb_sin_cos_aux prec (Suc n) i k x = 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   789
    (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   790
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   791
lemma cos_aux:
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   794
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   795
  have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   796
  let "?f n" = "fact (2 * n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   797
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   798
  { fix n 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   799
    have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   800
    have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 1 * (((\<lambda>i. i + 2) ^ n) 1 + 1)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   801
      unfolding F by auto } note f_eq = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   802
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   803
  from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   804
    OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   805
  show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "Ifloat x"])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   806
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   807
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   808
lemma cos_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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))}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   810
proof (cases "Ifloat x = 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   811
  case False hence "Ifloat x \<noteq> 0" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   812
  hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   813
  have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   814
    using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   815
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   816
  { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x^(2 * i))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   817
    = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   818
  proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   819
    have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   820
    also have "\<dots> = 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   821
      (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   822
    also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   823
      unfolding sum_split_even_odd ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   824
    also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   825
      by (rule setsum_cong2) auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   826
    finally show ?thesis by assumption
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   827
  qed } note morph_to_if_power = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   828
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   829
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   830
  { fix n :: nat assume "0 < n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   831
    hence "0 < 2 * n" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   832
    obtain t where "0 < t" and "t < Ifloat x" and
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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) 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   834
      + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (Ifloat x)^(2*n)" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   835
      (is "_ = ?SUM + ?rest / ?fact * ?pow")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   836
      using Maclaurin_cos_expansion2[OF `0 < Ifloat x` `0 < 2 * n`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   837
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   838
    have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   839
    also have "\<dots> = cos (t + real n * pi)"  using cos_add by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   840
    also have "\<dots> = ?rest" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   841
    finally have "cos t * -1^n = ?rest" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   842
    moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   843
    have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   844
    hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   845
    ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   846
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   847
    have "0 < ?fact" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   848
    have "0 < ?pow" using `0 < Ifloat x` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   849
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   850
    {
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   851
      assume "even n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   852
      have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   853
	unfolding morph_to_if_power[symmetric] using cos_aux by auto 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   854
      also have "\<dots> \<le> cos (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   855
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   856
	from even[OF `even n`] `0 < ?fact` `0 < ?pow`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   857
	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   858
	thus ?thesis unfolding cos_eq by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   859
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   860
      finally have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (Ifloat x)" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   861
    } note lb = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   862
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   863
    {
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   864
      assume "odd n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   865
      have "cos (Ifloat x) \<le> ?SUM"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   866
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   867
	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   868
	have "0 \<le> (- ?rest) / ?fact * ?pow"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   869
	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   870
	thus ?thesis unfolding cos_eq by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   871
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   872
      also have "\<dots> \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   873
	unfolding morph_to_if_power[symmetric] using cos_aux by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   874
      finally have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   875
    } note ub = this and lb
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   876
  } note ub = this(1) and lb = this(2)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   877
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   879
  moreover have "Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (Ifloat x)" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   880
  proof (cases "0 < get_even n")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   881
    case True show ?thesis using lb[OF True get_even] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   882
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   883
    case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   884
    hence "get_even n = 0" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   885
    have "- (pi / 2) \<le> Ifloat x" by (rule order_trans[OF _ `0 < Ifloat x`[THEN less_imp_le]], auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   886
    with `Ifloat x \<le> pi / 2`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   887
    show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps Ifloat_minus Ifloat_0 using cos_ge_zero by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   888
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   889
  ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   890
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   891
  case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   892
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   893
  proof (cases "n = 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   894
    case True 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   896
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   897
    case False with not0_implies_Suc obtain m where "n = Suc m" by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   899
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   900
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   901
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   902
lemma sin_aux: assumes "0 \<le> Ifloat x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   905
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   906
  have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   907
  let "?f n" = "fact (2 * n + 1)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   908
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   909
  { fix n 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   910
    have F: "\<And>m. ((\<lambda>i. i + 2) ^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   911
    have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^ n) 2 * (((\<lambda>i. i + 2) ^ n) 2 + 1)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   912
      unfolding F by auto } note f_eq = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   913
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   914
  from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   915
    OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   916
  show "?lb" and "?ub" using `0 \<le> Ifloat x` unfolding Ifloat_mult
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   917
    unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   918
    unfolding real_mult_commute
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   919
    by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "Ifloat x"])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   920
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   921
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   922
lemma sin_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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))}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   924
proof (cases "Ifloat x = 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   925
  case False hence "Ifloat x \<noteq> 0" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   926
  hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   927
  have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   928
    using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   929
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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 = _")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   932
    proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   933
      have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   934
      have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   936
	unfolding sum_split_even_odd ..
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   938
	by (rule setsum_cong2) auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   939
      finally show ?thesis by assumption
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   940
    qed } note setsum_morph = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   941
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   942
  { fix n :: nat assume "0 < n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   943
    hence "0 < 2 * n + 1" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   944
    obtain t where "0 < t" and "t < Ifloat x" and
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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) 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   946
      + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (Ifloat x)^(2*n + 1)" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   947
      (is "_ = ?SUM + ?rest / ?fact * ?pow")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   948
      using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < Ifloat x`] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   949
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   950
    have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   951
    moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   952
    have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   953
    hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   954
    ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   955
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   956
    have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   957
    have "0 < ?pow" using `0 < Ifloat x` by (rule zero_less_power)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   958
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   959
    {
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   960
      assume "even n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   961
      have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   962
            (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   963
	using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   964
      also have "\<dots> \<le> ?SUM" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   965
      also have "\<dots> \<le> sin (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   966
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   967
	from even[OF `even n`] `0 < ?fact` `0 < ?pow`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   968
	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   969
	thus ?thesis unfolding sin_eq by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   970
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   971
      finally have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (Ifloat x)" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   972
    } note lb = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   973
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   974
    {
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   975
      assume "odd n"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   976
      have "sin (Ifloat x) \<le> ?SUM"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   977
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   978
	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   979
	have "0 \<le> (- ?rest) / ?fact * ?pow"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   980
	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   981
	thus ?thesis unfolding sin_eq by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   982
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   984
	 by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   985
      also have "\<dots> \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   986
	using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   987
      finally have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   988
    } note ub = this and lb
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   989
  } note ub = this(1) and lb = this(2)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   990
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   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] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   992
  moreover have "Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (Ifloat x)" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   993
  proof (cases "0 < get_even n")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   994
    case True show ?thesis using lb[OF True get_even] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   995
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   996
    case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   997
    hence "get_even n = 0" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   998
    with `Ifloat x \<le> pi / 2` `0 \<le> Ifloat x`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
   999
    show ?thesis unfolding `get_even n = 0` ub_sin_cos_aux.simps Ifloat_minus Ifloat_0 using sin_ge_zero by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1000
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1001
  ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1002
next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1003
  case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1004
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1005
  proof (cases "n = 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1006
    case True 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1008
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1009
    case False with not0_implies_Suc obtain m where "n = Suc m" by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1011
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1012
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1013
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1014
subsection "Compute the cosinus in the entire domain"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1015
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1016
definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1017
"lb_cos prec x = (let
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1018
    horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1019
    half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1020
  in if x < Float 1 -1 then horner x
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1021
else if x < 1          then half (horner (x * Float 1 -1))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1022
                       else half (half (horner (x * Float 1 -2))))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1023
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1024
definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1025
"ub_cos prec x = (let
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1026
    horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1027
    half = \<lambda> x. Float 1 1 * x * x - 1
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1028
  in if x < Float 1 -1 then horner x
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1029
else if x < 1          then half (horner (x * Float 1 -1))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1030
                       else half (half (horner (x * Float 1 -2))))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1031
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1032
definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1033
"bnds_cos prec lx ux = (let  lpi = lb_pi prec
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1034
  in   if lx < -lpi \<or> ux > lpi   then (Float -1 0, Float 1 0)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1035
  else if ux \<le> 0                 then (lb_cos prec (-lx), ub_cos prec (-ux))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1036
  else if 0 \<le> lx                 then (lb_cos prec ux, ub_cos prec lx)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1037
                                 else (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1038
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1039
lemma lb_cos: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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) }")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1041
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1042
  { fix x :: real
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1043
    have "cos x = cos (x / 2 + x / 2)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1045
      unfolding cos_add by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1046
    also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1047
    finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1048
  } note x_half = this[symmetric]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1049
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1050
  have "\<not> x < 0" using `0 \<le> Ifloat x` unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1051
  let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1052
  let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1053
  let "?ub_half x" = "Float 1 1 * x * x - 1"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1054
  let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1055
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1056
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1057
  proof (cases "x < Float 1 -1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1058
    case True hence "Ifloat x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1060
      using cos_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1061
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1062
    case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1063
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1064
    { fix y x :: float let ?x2 = "Ifloat (x * Float 1 -1)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1065
      assume "Ifloat y \<le> cos ?x2" and "-pi \<le> Ifloat x" and "Ifloat x \<le> pi"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1066
      hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1067
      hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1068
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1069
      have "Ifloat (?lb_half y) \<le> cos (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1070
      proof (cases "y < 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1071
	case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1072
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1073
	case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1074
	hence "0 \<le> Ifloat y" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1075
	from mult_mono[OF `Ifloat y \<le> cos ?x2` `Ifloat y \<le> cos ?x2` `0 \<le> cos ?x2` this]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1076
	have "Ifloat y * Ifloat y \<le> cos ?x2 * cos ?x2" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1077
	hence "2 * Ifloat y * Ifloat y \<le> 2 * cos ?x2 * cos ?x2" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1079
	thus ?thesis unfolding if_not_P[OF False] x_half Float_num Ifloat_mult Ifloat_sub by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1080
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1081
    } note lb_half = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1082
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1083
    { fix y x :: float let ?x2 = "Ifloat (x * Float 1 -1)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1084
      assume ub: "cos ?x2 \<le> Ifloat y" and "- pi \<le> Ifloat x" and "Ifloat x \<le> pi"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1085
      hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding Ifloat_mult Float_num by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1086
      hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1087
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1088
      have "cos (Ifloat x) \<le> Ifloat (?ub_half y)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1089
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1090
	have "0 \<le> Ifloat y" using `0 \<le> cos ?x2` ub by (rule order_trans)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1091
	from mult_mono[OF ub ub this `0 \<le> cos ?x2`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1092
	have "cos ?x2 * cos ?x2 \<le> Ifloat y * Ifloat y" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1093
	hence "2 * cos ?x2 * cos ?x2 \<le> 2 * Ifloat y * Ifloat y" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1095
	thus ?thesis unfolding x_half Ifloat_mult Float_num Ifloat_sub by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1096
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1097
    } note ub_half = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1098
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1099
    let ?x2 = "x * Float 1 -1"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1100
    let ?x4 = "x * Float 1 -1 * Float 1 -1"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1101
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1103
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1104
    show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1105
    proof (cases "x < 1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1106
      case True hence "Ifloat x \<le> 1" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1108
      from cos_boundaries[OF this]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1109
      have lb: "Ifloat (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> Ifloat (?ub_horner ?x2)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1110
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1111
      have "Ifloat (?lb x) \<le> ?cos x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1112
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1113
	from lb_half[OF lb `-pi \<le> Ifloat x` `Ifloat x \<le> pi`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1115
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1116
      moreover have "?cos x \<le> Ifloat (?ub x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1117
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1118
	from ub_half[OF ub `-pi \<le> Ifloat x` `Ifloat x \<le> pi`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1120
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1121
      ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1122
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1123
      case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1125
      from cos_boundaries[OF this]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1126
      have lb: "Ifloat (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> Ifloat (?ub_horner ?x4)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1127
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1128
      have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1129
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1130
      have "Ifloat (?lb x) \<le> ?cos x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1131
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1133
	from lb_half[OF lb_half[OF lb this] `-pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1135
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1136
      moreover have "?cos x \<le> Ifloat (?ub x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1137
      proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1139
	from ub_half[OF ub_half[OF ub this] `-pi \<le> Ifloat x` `Ifloat x \<le> pi`, unfolded eq_4]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1141
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1142
      ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1143
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1144
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1145
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1146
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1147
lemma lb_cos_minus: assumes "-pi \<le> Ifloat x" and "Ifloat x \<le> 0" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1148
  shows "cos (Ifloat (-x)) \<in> {Ifloat (lb_cos prec (-x)) .. Ifloat (ub_cos prec (-x))}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1149
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1150
  have "0 \<le> Ifloat (-x)" and "Ifloat (-x) \<le> pi" using `-pi \<le> Ifloat x` `Ifloat x \<le> 0` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1151
  from lb_cos[OF this] show ?thesis .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1152
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1153
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1155
proof (rule allI, rule allI, rule allI, rule impI)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1156
  fix x lx ux
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1157
  assume "(l, u) = bnds_cos prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1158
  hence bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1159
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1160
  let ?lpi = "lb_pi prec"  
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1161
  have [intro!]: "Ifloat lx \<le> Ifloat ux" using x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1162
  hence "lx \<le> ux" unfolding le_float_def .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1163
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1164
  show "Ifloat l \<le> cos x \<and> cos x \<le> Ifloat u"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1165
  proof (cases "lx < -?lpi \<or> ux > ?lpi")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1166
    case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1168
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1169
    case False note not_out = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1171
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1172
    from pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1, THEN le_imp_neg_le] lpi_lx
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1173
    have "- pi \<le> Ifloat lx" by (rule order_trans)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1174
    hence "- pi \<le> x" and "- pi \<le> Ifloat ux" and "x \<le> Ifloat ux" using x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1175
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1176
    from lpi_ux pi_boundaries[unfolded atLeastAtMost_iff, THEN conjunct1]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1177
    have "Ifloat ux \<le> pi" by (rule order_trans)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1178
    hence "x \<le> pi" and "Ifloat lx \<le> pi" and "Ifloat lx \<le> x" using x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1179
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1180
    note lb_cos_minus_bottom = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct1]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1181
    note lb_cos_minus_top = lb_cos_minus[unfolded atLeastAtMost_iff, THEN conjunct2]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1182
    note lb_cos_bottom = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct1]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1183
    note lb_cos_top = lb_cos[unfolded atLeastAtMost_iff, THEN conjunct2]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1184
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1185
    show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1186
    proof (cases "ux \<le> 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1187
      case True hence "Ifloat ux \<le> 0" unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1188
      hence "x \<le> 0" and "Ifloat lx \<le> 0" using x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1189
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1192
	finally have "Ifloat (lb_cos prec (-lx)) \<le> cos x" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1193
      moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1195
	also have "\<dots> \<le> Ifloat (ub_cos prec (-ux))" using lb_cos_minus_top[OF `-pi \<le> Ifloat ux` `Ifloat ux \<le> 0`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1196
	finally have "cos x \<le> Ifloat (ub_cos prec (-ux))" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1197
      ultimately show ?thesis using bnds unfolding bnds_cos_def Let_def if_not_P[OF not_out] if_P[OF True] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1198
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1199
      case False note not_ux = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1200
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1201
      show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1202
      proof (cases "0 \<le> lx")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1203
	case True hence "0 \<le> Ifloat lx" unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1204
	hence "0 \<le> x" and "0 \<le> Ifloat ux" using x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1205
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1206
	{ have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1207
	  also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1208
	  finally have "Ifloat (lb_cos prec ux) \<le> cos x" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1209
	moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1211
	  also have "\<dots> \<le> Ifloat (ub_cos prec lx)" using lb_cos_top[OF `0 \<le> Ifloat lx` `Ifloat lx \<le> pi`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1212
	  finally have "cos x \<le> Ifloat (ub_cos prec lx)" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1214
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1215
	case False with not_ux
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1216
	have "Ifloat lx \<le> 0" and "0 \<le> Ifloat ux" unfolding le_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1217
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1218
	have "Ifloat (min (lb_cos prec (-lx)) (lb_cos prec ux)) \<le> cos x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1219
	proof (cases "x \<le> 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1220
	  case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1223
	  finally show ?thesis unfolding Ifloat_min by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1224
	next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1225
	  case False hence "0 \<le> x" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1226
	  have "Ifloat (lb_cos prec ux) \<le> cos (Ifloat ux)" using lb_cos_bottom[OF `0 \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1227
	  also have "\<dots> \<le> cos x" using cos_monotone_0_pi'[OF `0 \<le> x` `x \<le> Ifloat ux` `Ifloat ux \<le> pi`] .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1228
	  finally show ?thesis unfolding Ifloat_min by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1229
	qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1230
	moreover have "cos x \<le> Ifloat (Float 1 0)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1232
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1233
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1234
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1235
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1236
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1237
subsection "Compute the sinus in the entire domain"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1238
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1239
function lb_sin :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_sin :: "nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1240
"lb_sin prec x = (let sqr_diff = \<lambda> x. if x > 1 then 0 else 1 - x * x 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1241
  in if x < 0           then - ub_sin prec (- x)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1243
                        else the (lb_sqrt prec (sqr_diff (ub_cos prec x))))" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1244
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1245
"ub_sin prec x = (let sqr_diff = \<lambda> x. if x < 0 then 1 else 1 - x * x
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1246
  in if x < 0           then - lb_sin prec (- x)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1248
                        else the (ub_sqrt prec (sqr_diff (lb_cos prec x))))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1249
by pat_completeness auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1251
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1252
definition bnds_sin :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1253
"bnds_sin prec lx ux = (let 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1254
    lpi = lb_pi prec ;
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1255
    half_pi = lpi * Float 1 -1
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1256
  in if lx \<le> - half_pi \<or> half_pi \<le> ux then (Float -1 0, Float 1 0)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1257
                                       else (lb_sin prec lx, ub_sin prec ux))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1258
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1259
lemma lb_sin: assumes "- (pi / 2) \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1260
  shows "sin (Ifloat x) \<in> { Ifloat (lb_sin prec x) .. Ifloat (ub_sin prec x) }" (is "?sin x \<in> { ?lb x .. ?ub x}")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1261
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1262
  { fix x :: float assume "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1263
    hence "\<not> (x < 0)" and "- (pi / 2) \<le> Ifloat x" unfolding less_float_def using pi_ge_two by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1264
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1265
    have "Ifloat x \<le> pi" using `Ifloat x \<le> pi / 2` using pi_ge_two by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1266
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1267
    have "?sin x \<in> { ?lb x .. ?ub x}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1268
    proof (cases "x \<le> Float 1 -1")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1269
      case True from sin_boundaries[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi / 2`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1271
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1272
      case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1274
      have "0 \<le> sin (Ifloat x)" using `0 \<le> Ifloat x` and `Ifloat x \<le> pi / 2` using sin_ge_zero by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1275
      
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1276
      have "?sin x \<le> ?ub x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1277
      proof (cases "lb_cos prec x < 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1278
	case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1279
	have "?sin x \<le> 1" using sin_le_one .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1282
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1283
	case False hence "0 \<le> Ifloat (lb_cos prec x)" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1284
	
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1286
	also have "\<dots> \<le> sqrt (Ifloat (1 - lb_cos prec x * lb_cos prec x))" 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1287
	proof (rule real_sqrt_le_mono)
30273
ecd6f0ca62ea declare power_Suc [simp]; remove redundant type-specific 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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1291
	qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1292
	also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 - lb_cos prec x * lb_cos prec x)))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1293
	proof (rule ub_sqrt_lower_bound)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1295
	  from mult_mono[OF order_trans[OF this cos_le_one] order_trans[OF this cos_le_one]]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1296
	  have "Ifloat (lb_cos prec x) * Ifloat (lb_cos prec x) \<le> 1" using `0 \<le> Ifloat (lb_cos prec x)` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1297
	  thus "0 \<le> Ifloat (1 - lb_cos prec x * lb_cos prec x)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1298
	qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1300
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1301
      moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1302
      have "?lb x \<le> ?sin x"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1303
      proof (cases "1 < ub_cos prec x")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1304
	case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1306
	  by (rule order_trans[OF _ sin_ge_zero[OF `0 \<le> Ifloat x` `Ifloat x \<le> pi`]]) 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1307
        (auto simp add: lb_sqrt_upper_bound[where prec=prec and x=0, unfolded Ifloat_0 real_sqrt_zero])
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1308
      next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1309
	case False hence "Ifloat (ub_cos prec x) \<le> 1" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1311
	
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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))"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1313
	proof (rule lb_sqrt_upper_bound)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)`
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1315
	  have "Ifloat (ub_cos prec x) * Ifloat (ub_cos prec x) \<le> 1" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1316
	  thus "0 \<le> Ifloat (1 - ub_cos prec x * ub_cos prec x)" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1317
	qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1318
	also have "\<dots> \<le> sqrt (1 - cos (Ifloat x) ^ 2)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1319
	proof (rule real_sqrt_le_mono)
30273
ecd6f0ca62ea declare power_Suc [simp]; remove redundant type-specific 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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1323
	qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1324
	also have "\<dots> = sin (Ifloat x)" unfolding sin_squared_eq[symmetric] real_sqrt_abs using `0 \<le> sin (Ifloat x)` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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 .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1326
      qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1327
      ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1328
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1329
  } note for_pos = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1330
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1331
  show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1332
  proof (cases "x < 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1333
    case True 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1334
    hence "0 \<le> Ifloat (-x)" and "Ifloat (- x) \<le> pi / 2" using `-(pi/2) \<le> Ifloat x` unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1335
    from for_pos[OF this]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1337
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1338
    case False hence "0 \<le> Ifloat x" unfolding less_float_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1339
    from for_pos[OF this `Ifloat x \<le> pi /2`]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1340
    show ?thesis .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1341
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1342
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1343
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1345
proof (rule allI, rule allI, rule allI, rule impI)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1346
  fix x lx ux
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1347
  assume "(l, u) = bnds_sin prec lx ux \<and> x \<in> {Ifloat lx .. Ifloat ux}"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1348
  hence bnds: "(l, u) = bnds_sin prec lx ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1349
  show "Ifloat l \<le> sin x \<and> sin x \<le> Ifloat u"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1350
  proof (cases "lx \<le> - (lb_pi prec * Float 1 -1) \<or> lb_pi prec * Float 1 -1 \<le> ux")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1351
    case True show ?thesis using bnds unfolding bnds_sin_def if_P[OF True] Let_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1352
  next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1353
    case False
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1355
    moreover have "Ifloat (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding Ifloat_mult using pi_boundaries by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1357
    hence "- (pi / 2) \<le> Ifloat ux" and "Ifloat lx \<le> pi / 2" by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1358
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1359
    have "- (pi / 2) \<le> x""x \<le> pi / 2" using `Ifloat ux \<le> pi / 2` `- (pi /2) \<le> Ifloat lx` x by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1360
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1362
      also have "\<dots> \<le> sin x" using sin_monotone_2pi' `- (pi / 2) \<le> Ifloat lx` x `x \<le> pi / 2` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1363
      finally have "Ifloat (lb_sin prec lx) \<le> sin x" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1364
    moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1365
    { have "sin x \<le> sin (Ifloat ux)" using sin_monotone_2pi' `- (pi / 2) \<le> x` x `Ifloat ux \<le> pi / 2` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1367
      finally have "sin x \<le> Ifloat (ub_sin prec ux)" . }
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1368
    ultimately
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1369
    show ?thesis using bnds unfolding bnds_sin_def if_not_P[OF False] Let_def by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1370
  qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1371
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1372
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1373
section "Exponential function"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1374
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1375
subsection "Compute the series of the exponential function"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1376
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1378
"ub_exp_horner prec 0 i k x       = 0" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1380
"lb_exp_horner prec 0 i k x       = 0" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1382
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1383
lemma bnds_exp_horner: assumes "Ifloat x \<le> 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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) }"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1385
proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1386
  { fix n
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1387
    have F: "\<And> m. ((\<lambda>i. i + 1) ^ n) m = n + m" by (induct n, auto)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1388
    have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^ n) 1" unfolding F by auto } note f_eq = this
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1389
    
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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,
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1391
    OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1392
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1394
      using bounds(1) by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1395
    also have "\<dots> \<le> exp (Ifloat x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1396
    proof -
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1398
	using Maclaurin_exp_le by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1399
      moreover have "0 \<le> exp t / real (fact (get_even n)) * (Ifloat x) ^ (get_even n)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1400
	by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1401
      ultimately show ?thesis
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1402
	using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1403
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1404
    finally have "Ifloat (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (Ifloat x)" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1405
  } moreover
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1406
  { 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1407
    have x_less_zero: "Ifloat x ^ get_odd n \<le> 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1408
    proof (cases "Ifloat x = 0")
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1409
      case True
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1410
      have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1411
      thus ?thesis unfolding True power_0_left by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1412
    next
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1413
      case False hence "Ifloat x < 0" using `Ifloat x \<le> 0` by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1414
      show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `Ifloat x < 0`)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1415
    qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1416
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1418
      using Maclaurin_exp_le by blast
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1419
    moreover have "exp t / real (fact (get_odd n)) * (Ifloat x) ^ (get_odd n) \<le> 0"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1420
      by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1421
    ultimately have "exp (Ifloat x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * Ifloat x ^ j)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1422
      using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1423
    also have "\<dots> \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1424
      using bounds(2) by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1425
    finally have "exp (Ifloat x) \<le> Ifloat (ub_exp_horner prec (get_odd n) 1 1 x)" .
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1426
  } ultimately show ?thesis by auto
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1427
qed
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1428
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1429
subsection "Compute the exponential function on the entire domain"
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1430
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1431
function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1432
"lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1433
             else let 
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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)
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  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))
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1436
                           else horner x)" |
a5da150bd0ab Add approximation method
hoelzl
parents:
diff changeset
  1437
"ub_exp prec x = (if 0 < x    then float_divr prec 1 (lb_exp prec (-x))
a5da150bd0ab Add approximation method