src/HOL/Number_Theory/Fib.thy
changeset 64317 029e6247210e
parent 64267 b9a1486e79be
child 65393 079a6f850c02
     1.1 --- a/src/HOL/Number_Theory/Fib.thy	Thu Oct 20 11:05:16 2016 +0200
     1.2 +++ b/src/HOL/Number_Theory/Fib.thy	Thu Oct 20 13:53:36 2016 +0200
     1.3 @@ -1,12 +1,13 @@
     1.4  (*  Title:      HOL/Number_Theory/Fib.thy
     1.5      Author:     Lawrence C. Paulson
     1.6      Author:     Jeremy Avigad
     1.7 +    Author:     Manuel Eberl
     1.8  *)
     1.9  
    1.10  section \<open>The fibonacci function\<close>
    1.11  
    1.12  theory Fib
    1.13 -imports Main GCD Binomial
    1.14 +  imports Complex_Main
    1.15  begin
    1.16  
    1.17  
    1.18 @@ -38,6 +39,34 @@
    1.19    by (induct n rule: fib.induct) (auto simp add: )
    1.20  
    1.21  
    1.22 +subsection \<open>More efficient code\<close>
    1.23 +
    1.24 +text \<open>
    1.25 +  The naive approach is very inefficient since the branching recursion leads to many
    1.26 +  values of @{term fib} being computed multiple times. We can avoid this by ``remembering''
    1.27 +  the last two values in the sequence, yielding a tail-recursive version.
    1.28 +  This is far from optimal (it takes roughly $O(n\cdot M(n))$ time where $M(n)$ is the 
    1.29 +  time required to multiply two $n$-bit integers), but much better than the naive version,
    1.30 +  which is exponential.
    1.31 +\<close>
    1.32 +
    1.33 +fun gen_fib :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat" where
    1.34 +  "gen_fib a b 0 = a"
    1.35 +| "gen_fib a b (Suc 0) = b"
    1.36 +| "gen_fib a b (Suc (Suc n)) = gen_fib b (a + b) (Suc n)"
    1.37 +
    1.38 +lemma gen_fib_recurrence: "gen_fib a b (Suc (Suc n)) = gen_fib a b n + gen_fib a b (Suc n)"
    1.39 +  by (induction a b n rule: gen_fib.induct) simp_all
    1.40 +  
    1.41 +lemma gen_fib_fib: "gen_fib (fib n) (fib (Suc n)) m = fib (n + m)"
    1.42 +  by (induction m rule: fib.induct) (simp_all del: gen_fib.simps(3) add: gen_fib_recurrence)
    1.43 +
    1.44 +lemma fib_conv_gen_fib: "fib n = gen_fib 0 1 n"
    1.45 +  using gen_fib_fib[of 0 n] by simp
    1.46 +
    1.47 +declare fib_conv_gen_fib [code]
    1.48 +
    1.49 +
    1.50  subsection \<open>A Few Elementary Results\<close>
    1.51  
    1.52  text \<open>
    1.53 @@ -104,6 +133,114 @@
    1.54    by (induct n rule: nat.induct) (auto simp add:  field_simps)
    1.55  
    1.56  
    1.57 +subsection \<open>Closed form\<close>
    1.58 +
    1.59 +lemma fib_closed_form:
    1.60 +  defines "\<phi> \<equiv> (1 + sqrt 5) / (2::real)" and "\<psi> \<equiv> (1 - sqrt 5) / (2::real)"
    1.61 +  shows   "of_nat (fib n) = (\<phi> ^ n - \<psi> ^ n) / sqrt 5"
    1.62 +proof (induction n rule: fib.induct)
    1.63 +  fix n :: nat
    1.64 +  assume IH1: "of_nat (fib n) = (\<phi> ^ n - \<psi> ^ n) / sqrt 5"
    1.65 +  assume IH2: "of_nat (fib (Suc n)) = (\<phi> ^ Suc n - \<psi> ^ Suc n) / sqrt 5"
    1.66 +  have "of_nat (fib (Suc (Suc n))) = of_nat (fib (Suc n)) + of_nat (fib n)" by simp
    1.67 +  also have "... = (\<phi>^n*(\<phi> + 1) - \<psi>^n*(\<psi> + 1)) / sqrt 5"
    1.68 +    by (simp add: IH1 IH2 field_simps)
    1.69 +  also have "\<phi> + 1 = \<phi>\<^sup>2" by (simp add: \<phi>_def field_simps power2_eq_square)
    1.70 +  also have "\<psi> + 1 = \<psi>\<^sup>2" by (simp add: \<psi>_def field_simps power2_eq_square)
    1.71 +  also have "\<phi>^n * \<phi>\<^sup>2 - \<psi>^n * \<psi>\<^sup>2 = \<phi> ^ Suc (Suc n) - \<psi> ^ Suc (Suc n)"  
    1.72 +    by (simp add: power2_eq_square)
    1.73 +  finally show "of_nat (fib (Suc (Suc n))) = (\<phi> ^ Suc (Suc n) - \<psi> ^ Suc (Suc n)) / sqrt 5" .
    1.74 +qed (simp_all add: \<phi>_def \<psi>_def field_simps)
    1.75 +
    1.76 +lemma fib_closed_form':
    1.77 +  defines "\<phi> \<equiv> (1 + sqrt 5) / (2 :: real)" and "\<psi> \<equiv> (1 - sqrt 5) / (2 :: real)"
    1.78 +  assumes "n > 0"
    1.79 +  shows   "fib n = round (\<phi> ^ n / sqrt 5)"
    1.80 +proof (rule sym, rule round_unique')
    1.81 +  have "\<bar>\<phi> ^ n / sqrt 5 - of_int (int (fib n))\<bar> = \<bar>\<psi>\<bar> ^ n / sqrt 5"
    1.82 +    by (simp add: fib_closed_form[folded \<phi>_def \<psi>_def] field_simps power_abs)
    1.83 +  also {
    1.84 +    from assms have "\<bar>\<psi>\<bar>^n \<le> \<bar>\<psi>\<bar>^1"
    1.85 +      by (intro power_decreasing) (simp_all add: algebra_simps real_le_lsqrt)
    1.86 +    also have "... < sqrt 5 / 2" by (simp add: \<psi>_def field_simps)
    1.87 +    finally have "\<bar>\<psi>\<bar>^n / sqrt 5 < 1/2" by (simp add: field_simps)
    1.88 +  }
    1.89 +  finally show "\<bar>\<phi> ^ n / sqrt 5 - of_int (int (fib n))\<bar> < 1/2" .
    1.90 +qed
    1.91 +
    1.92 +lemma fib_asymptotics:
    1.93 +  defines "\<phi> \<equiv> (1 + sqrt 5) / (2 :: real)"
    1.94 +  shows   "(\<lambda>n. real (fib n) / (\<phi> ^ n / sqrt 5)) \<longlonglongrightarrow> 1"
    1.95 +proof -
    1.96 +  define \<psi> where "\<psi> \<equiv> (1 - sqrt 5) / (2 :: real)"
    1.97 +  have "\<phi> > 1" by (simp add: \<phi>_def)
    1.98 +  hence A: "\<phi> \<noteq> 0" by auto
    1.99 +  have "(\<lambda>n. (\<psi> / \<phi>) ^ n) \<longlonglongrightarrow> 0"
   1.100 +    by (rule LIMSEQ_power_zero) (simp_all add: \<phi>_def \<psi>_def field_simps add_pos_pos)
   1.101 +  hence "(\<lambda>n. 1 - (\<psi> / \<phi>) ^ n) \<longlonglongrightarrow> 1 - 0" by (intro tendsto_diff tendsto_const)
   1.102 +  with A show ?thesis
   1.103 +    by (simp add: divide_simps fib_closed_form [folded \<phi>_def \<psi>_def])
   1.104 +qed
   1.105 +
   1.106 +
   1.107 +subsection \<open>Divide-and-Conquer recurrence\<close>
   1.108 +
   1.109 +text \<open>
   1.110 +  The following divide-and-conquer recurrence allows for a more efficient computation 
   1.111 +  of Fibonacci numbers; however, it requires memoisation of values to be reasonably 
   1.112 +  efficient, cutting the number of values to be computed to logarithmically many instead of
   1.113 +  linearly many. The vast majority of the computation time is then actually spent on the 
   1.114 +  multiplication, since the output number is exponential in the input number.
   1.115 +\<close>
   1.116 +
   1.117 +lemma fib_rec_odd:
   1.118 +  defines "\<phi> \<equiv> (1 + sqrt 5) / (2 :: real)" and "\<psi> \<equiv> (1 - sqrt 5) / (2 :: real)"
   1.119 +  shows   "fib (Suc (2*n)) = fib n^2 + fib (Suc n)^2"
   1.120 +proof -
   1.121 +  have "of_nat (fib n^2 + fib (Suc n)^2) = ((\<phi> ^ n - \<psi> ^ n)\<^sup>2 + (\<phi> * \<phi> ^ n - \<psi> * \<psi> ^ n)\<^sup>2)/5"
   1.122 +    by (simp add: fib_closed_form[folded \<phi>_def \<psi>_def] field_simps power2_eq_square)
   1.123 +  also have "(\<phi> ^ n - \<psi> ^ n)\<^sup>2 + (\<phi> * \<phi> ^ n - \<psi> * \<psi> ^ n)\<^sup>2 = 
   1.124 +    \<phi>^(2*n) + \<psi>^(2*n) - 2*(\<phi>*\<psi>)^n + \<phi>^(2*n+2) + \<psi>^(2*n+2) - 2*(\<phi>*\<psi>)^(n+1)" (is "_ = ?A")
   1.125 +      by (simp add: power2_eq_square algebra_simps power_mult power_mult_distrib)
   1.126 +  also have "\<phi> * \<psi> = -1" by (simp add: \<phi>_def \<psi>_def field_simps)
   1.127 +  hence "?A = \<phi>^(2*n+1) * (\<phi> + inverse \<phi>) + \<psi>^(2*n+1) * (\<psi> + inverse \<psi>)" 
   1.128 +    by (auto simp: field_simps power2_eq_square)
   1.129 +  also have "1 + sqrt 5 > 0" by (auto intro: add_pos_pos)
   1.130 +  hence "\<phi> + inverse \<phi> = sqrt 5" by (simp add: \<phi>_def field_simps)
   1.131 +  also have "\<psi> + inverse \<psi> = -sqrt 5" by (simp add: \<psi>_def field_simps)
   1.132 +  also have "(\<phi> ^ (2*n+1) * sqrt 5 + \<psi> ^ (2*n+1)* - sqrt 5) / 5 =
   1.133 +               (\<phi> ^ (2*n+1) - \<psi> ^ (2*n+1)) * (sqrt 5 / 5)" by (simp add: field_simps)
   1.134 +  also have "sqrt 5 / 5 = inverse (sqrt 5)" by (simp add: field_simps)
   1.135 +  also have "(\<phi> ^ (2*n+1) - \<psi> ^ (2*n+1)) * ... = of_nat (fib (Suc (2*n)))"
   1.136 +    by (simp add: fib_closed_form[folded \<phi>_def \<psi>_def] divide_inverse)
   1.137 +  finally show ?thesis by (simp only: of_nat_eq_iff)
   1.138 +qed
   1.139 +
   1.140 +lemma fib_rec_even: "fib (2*n) = (fib (n - 1) + fib (n + 1)) * fib n"
   1.141 +proof (induction n)
   1.142 +  case (Suc n)
   1.143 +  let ?rfib = "\<lambda>x. real (fib x)"
   1.144 +  have "2 * (Suc n) = Suc (Suc (2*n))" by simp
   1.145 +  also have "real (fib ...) = ?rfib n^2 + ?rfib (Suc n)^2 + (?rfib (n - 1) + ?rfib (n + 1)) * ?rfib n" 
   1.146 +    by (simp add: fib_rec_odd Suc)
   1.147 +  also have "(?rfib (n - 1) + ?rfib (n + 1)) * ?rfib n = (2 * ?rfib (n + 1) - ?rfib n) * ?rfib n"
   1.148 +    by (cases n) simp_all
   1.149 +  also have "?rfib n^2 + ?rfib (Suc n)^2 + ... = (?rfib (Suc n) + 2 * ?rfib n) * ?rfib (Suc n)"
   1.150 +    by (simp add: algebra_simps power2_eq_square)
   1.151 +  also have "... = real ((fib (Suc n - 1) + fib (Suc n + 1)) * fib (Suc n))" by simp
   1.152 +  finally show ?case by (simp only: of_nat_eq_iff)
   1.153 +qed simp
   1.154 +
   1.155 +lemma fib_rec_even': "fib (2*n) = (2*fib (n - 1) + fib n) * fib n"
   1.156 +  by (subst fib_rec_even, cases n) simp_all
   1.157 +
   1.158 +lemma fib_rec:
   1.159 +  "fib n = (if n = 0 then 0 else if n = 1 then 1 else
   1.160 +            if even n then let n' = n div 2; fn = fib n' in (2 * fib (n' - 1) + fn) * fn
   1.161 +            else let n' = n div 2 in fib n' ^ 2 + fib (Suc n') ^ 2)"
   1.162 +  by (auto elim: evenE oddE simp: fib_rec_odd fib_rec_even' Let_def)
   1.163 +
   1.164 +
   1.165  subsection \<open>Fibonacci and Binomial Coefficients\<close>
   1.166  
   1.167  lemma sum_drop_zero: "(\<Sum>k = 0..Suc n. if 0<k then (f (k - 1)) else 0) = (\<Sum>j = 0..n. f j)"