src/HOL/Number_Theory/Fib.thy
author haftmann
Sat Nov 11 18:41:08 2017 +0000 (19 months ago)
changeset 67051 e7e54a0b9197
parent 65393 079a6f850c02
child 69597 ff784d5a5bfb
permissions -rw-r--r--
dedicated definition for coprimality
     1 (*  Title:      HOL/Number_Theory/Fib.thy
     2     Author:     Lawrence C. Paulson
     3     Author:     Jeremy Avigad
     4     Author:     Manuel Eberl
     5 *)
     6 
     7 section \<open>The fibonacci function\<close>
     8 
     9 theory Fib
    10   imports Complex_Main
    11 begin
    12 
    13 
    14 subsection \<open>Fibonacci numbers\<close>
    15 
    16 fun fib :: "nat \<Rightarrow> nat"
    17   where
    18     fib0: "fib 0 = 0"
    19   | fib1: "fib (Suc 0) = 1"
    20   | fib2: "fib (Suc (Suc n)) = fib (Suc n) + fib n"
    21 
    22 
    23 subsection \<open>Basic Properties\<close>
    24 
    25 lemma fib_1 [simp]: "fib 1 = 1"
    26   by (metis One_nat_def fib1)
    27 
    28 lemma fib_2 [simp]: "fib 2 = 1"
    29   using fib.simps(3) [of 0] by (simp add: numeral_2_eq_2)
    30 
    31 lemma fib_plus_2: "fib (n + 2) = fib (n + 1) + fib n"
    32   by (metis Suc_eq_plus1 add_2_eq_Suc' fib.simps(3))
    33 
    34 lemma fib_add: "fib (Suc (n + k)) = fib (Suc k) * fib (Suc n) + fib k * fib n"
    35   by (induct n rule: fib.induct) (auto simp add: field_simps)
    36 
    37 lemma fib_neq_0_nat: "n > 0 \<Longrightarrow> fib n > 0"
    38   by (induct n rule: fib.induct) auto
    39 
    40 
    41 subsection \<open>More efficient code\<close>
    42 
    43 text \<open>
    44   The naive approach is very inefficient since the branching recursion leads to many
    45   values of @{term fib} being computed multiple times. We can avoid this by ``remembering''
    46   the last two values in the sequence, yielding a tail-recursive version.
    47   This is far from optimal (it takes roughly $O(n\cdot M(n))$ time where $M(n)$ is the
    48   time required to multiply two $n$-bit integers), but much better than the naive version,
    49   which is exponential.
    50 \<close>
    51 
    52 fun gen_fib :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat"
    53   where
    54     "gen_fib a b 0 = a"
    55   | "gen_fib a b (Suc 0) = b"
    56   | "gen_fib a b (Suc (Suc n)) = gen_fib b (a + b) (Suc n)"
    57 
    58 lemma gen_fib_recurrence: "gen_fib a b (Suc (Suc n)) = gen_fib a b n + gen_fib a b (Suc n)"
    59   by (induct a b n rule: gen_fib.induct) simp_all
    60 
    61 lemma gen_fib_fib: "gen_fib (fib n) (fib (Suc n)) m = fib (n + m)"
    62   by (induct m rule: fib.induct) (simp_all del: gen_fib.simps(3) add: gen_fib_recurrence)
    63 
    64 lemma fib_conv_gen_fib: "fib n = gen_fib 0 1 n"
    65   using gen_fib_fib[of 0 n] by simp
    66 
    67 declare fib_conv_gen_fib [code]
    68 
    69 
    70 subsection \<open>A Few Elementary Results\<close>
    71 
    72 text \<open>
    73   \<^medskip> Concrete Mathematics, page 278: Cassini's identity.  The proof is
    74   much easier using integers, not natural numbers!
    75 \<close>
    76 
    77 lemma fib_Cassini_int: "int (fib (Suc (Suc n)) * fib n) - int((fib (Suc n))\<^sup>2) = - ((-1)^n)"
    78   by (induct n rule: fib.induct) (auto simp add: field_simps power2_eq_square power_add)
    79 
    80 lemma fib_Cassini_nat:
    81   "fib (Suc (Suc n)) * fib n =
    82      (if even n then (fib (Suc n))\<^sup>2 - 1 else (fib (Suc n))\<^sup>2 + 1)"
    83   using fib_Cassini_int [of n] by (auto simp del: of_nat_mult of_nat_power)
    84 
    85 
    86 subsection \<open>Law 6.111 of Concrete Mathematics\<close>
    87 
    88 lemma coprime_fib_Suc_nat: "coprime (fib n) (fib (Suc n))"
    89   apply (induct n rule: fib.induct)
    90     apply (simp_all add: coprime_iff_gcd_eq_1 algebra_simps)
    91   apply (simp add: add.assoc [symmetric])
    92   done
    93 
    94 lemma gcd_fib_add:
    95   "gcd (fib m) (fib (n + m)) = gcd (fib m) (fib n)"
    96 proof (cases m)
    97   case 0
    98   then show ?thesis
    99     by simp
   100 next
   101   case (Suc q)
   102   from coprime_fib_Suc_nat [of q]
   103   have "coprime (fib (Suc q)) (fib q)"
   104     by (simp add: ac_simps)
   105   have "gcd (fib q) (fib (Suc q)) = Suc 0"
   106     using coprime_fib_Suc_nat [of q] by simp
   107   then have *: "gcd (fib n * fib q) (fib n * fib (Suc q)) = fib n"
   108     by (simp add: gcd_mult_distrib_nat [symmetric])
   109   moreover have "gcd (fib (Suc q)) (fib n * fib q + fib (Suc n) * fib (Suc q)) =
   110     gcd (fib (Suc q)) (fib n * fib q)"
   111     using gcd_add_mult [of "fib (Suc q)"] by (simp add: ac_simps)
   112   moreover have "gcd (fib (Suc q)) (fib n * fib (Suc q)) = fib (Suc q)"
   113     by simp
   114   ultimately show ?thesis
   115     using Suc \<open>coprime (fib (Suc q)) (fib q)\<close>
   116     by (auto simp add: fib_add algebra_simps gcd_mult_right_right_cancel)
   117 qed
   118 
   119 lemma gcd_fib_diff: "m \<le> n \<Longrightarrow> gcd (fib m) (fib (n - m)) = gcd (fib m) (fib n)"
   120   by (simp add: gcd_fib_add [symmetric, of _ "n-m"])
   121 
   122 lemma gcd_fib_mod: "0 < m \<Longrightarrow> gcd (fib m) (fib (n mod m)) = gcd (fib m) (fib n)"
   123 proof (induct n rule: less_induct)
   124   case (less n)
   125   show "gcd (fib m) (fib (n mod m)) = gcd (fib m) (fib n)"
   126   proof (cases "m < n")
   127     case True
   128     then have "m \<le> n" by auto
   129     with \<open>0 < m\<close> have "0 < n" by auto
   130     with \<open>0 < m\<close> \<open>m < n\<close> have *: "n - m < n" by auto
   131     have "gcd (fib m) (fib (n mod m)) = gcd (fib m) (fib ((n - m) mod m))"
   132       by (simp add: mod_if [of n]) (use \<open>m < n\<close> in auto)
   133     also have "\<dots> = gcd (fib m)  (fib (n - m))"
   134       by (simp add: less.hyps * \<open>0 < m\<close>)
   135     also have "\<dots> = gcd (fib m) (fib n)"
   136       by (simp add: gcd_fib_diff \<open>m \<le> n\<close>)
   137     finally show "gcd (fib m) (fib (n mod m)) = gcd (fib m) (fib n)" .
   138   next
   139     case False
   140     then show "gcd (fib m) (fib (n mod m)) = gcd (fib m) (fib n)"
   141       by (cases "m = n") auto
   142   qed
   143 qed
   144 
   145 lemma fib_gcd: "fib (gcd m n) = gcd (fib m) (fib n)"  \<comment> \<open>Law 6.111\<close>
   146   by (induct m n rule: gcd_nat_induct) (simp_all add: gcd_non_0_nat gcd.commute gcd_fib_mod)
   147 
   148 theorem fib_mult_eq_sum_nat: "fib (Suc n) * fib n = (\<Sum>k \<in> {..n}. fib k * fib k)"
   149   by (induct n rule: nat.induct) (auto simp add:  field_simps)
   150 
   151 
   152 subsection \<open>Closed form\<close>
   153 
   154 lemma fib_closed_form:
   155   fixes \<phi> \<psi> :: real
   156   defines "\<phi> \<equiv> (1 + sqrt 5) / 2"
   157     and "\<psi> \<equiv> (1 - sqrt 5) / 2"
   158   shows "of_nat (fib n) = (\<phi> ^ n - \<psi> ^ n) / sqrt 5"
   159 proof (induct n rule: fib.induct)
   160   fix n :: nat
   161   assume IH1: "of_nat (fib n) = (\<phi> ^ n - \<psi> ^ n) / sqrt 5"
   162   assume IH2: "of_nat (fib (Suc n)) = (\<phi> ^ Suc n - \<psi> ^ Suc n) / sqrt 5"
   163   have "of_nat (fib (Suc (Suc n))) = of_nat (fib (Suc n)) + of_nat (fib n)" by simp
   164   also have "\<dots> = (\<phi>^n * (\<phi> + 1) - \<psi>^n * (\<psi> + 1)) / sqrt 5"
   165     by (simp add: IH1 IH2 field_simps)
   166   also have "\<phi> + 1 = \<phi>\<^sup>2" by (simp add: \<phi>_def field_simps power2_eq_square)
   167   also have "\<psi> + 1 = \<psi>\<^sup>2" by (simp add: \<psi>_def field_simps power2_eq_square)
   168   also have "\<phi>^n * \<phi>\<^sup>2 - \<psi>^n * \<psi>\<^sup>2 = \<phi> ^ Suc (Suc n) - \<psi> ^ Suc (Suc n)"
   169     by (simp add: power2_eq_square)
   170   finally show "of_nat (fib (Suc (Suc n))) = (\<phi> ^ Suc (Suc n) - \<psi> ^ Suc (Suc n)) / sqrt 5" .
   171 qed (simp_all add: \<phi>_def \<psi>_def field_simps)
   172 
   173 lemma fib_closed_form':
   174   fixes \<phi> \<psi> :: real
   175   defines "\<phi> \<equiv> (1 + sqrt 5) / 2"
   176     and "\<psi> \<equiv> (1 - sqrt 5) / 2"
   177   assumes "n > 0"
   178   shows "fib n = round (\<phi> ^ n / sqrt 5)"
   179 proof (rule sym, rule round_unique')
   180   have "\<bar>\<phi> ^ n / sqrt 5 - of_int (int (fib n))\<bar> = \<bar>\<psi>\<bar> ^ n / sqrt 5"
   181     by (simp add: fib_closed_form[folded \<phi>_def \<psi>_def] field_simps power_abs)
   182   also {
   183     from assms have "\<bar>\<psi>\<bar>^n \<le> \<bar>\<psi>\<bar>^1"
   184       by (intro power_decreasing) (simp_all add: algebra_simps real_le_lsqrt)
   185     also have "\<dots> < sqrt 5 / 2" by (simp add: \<psi>_def field_simps)
   186     finally have "\<bar>\<psi>\<bar>^n / sqrt 5 < 1/2" by (simp add: field_simps)
   187   }
   188   finally show "\<bar>\<phi> ^ n / sqrt 5 - of_int (int (fib n))\<bar> < 1/2" .
   189 qed
   190 
   191 lemma fib_asymptotics:
   192   fixes \<phi> :: real
   193   defines "\<phi> \<equiv> (1 + sqrt 5) / 2"
   194   shows "(\<lambda>n. real (fib n) / (\<phi> ^ n / sqrt 5)) \<longlonglongrightarrow> 1"
   195 proof -
   196   define \<psi> :: real where "\<psi> \<equiv> (1 - sqrt 5) / 2"
   197   have "\<phi> > 1" by (simp add: \<phi>_def)
   198   then have *: "\<phi> \<noteq> 0" by auto
   199   have "(\<lambda>n. (\<psi> / \<phi>) ^ n) \<longlonglongrightarrow> 0"
   200     by (rule LIMSEQ_power_zero) (simp_all add: \<phi>_def \<psi>_def field_simps add_pos_pos)
   201   then have "(\<lambda>n. 1 - (\<psi> / \<phi>) ^ n) \<longlonglongrightarrow> 1 - 0"
   202     by (intro tendsto_diff tendsto_const)
   203   with * show ?thesis
   204     by (simp add: divide_simps fib_closed_form [folded \<phi>_def \<psi>_def])
   205 qed
   206 
   207 
   208 subsection \<open>Divide-and-Conquer recurrence\<close>
   209 
   210 text \<open>
   211   The following divide-and-conquer recurrence allows for a more efficient computation
   212   of Fibonacci numbers; however, it requires memoisation of values to be reasonably
   213   efficient, cutting the number of values to be computed to logarithmically many instead of
   214   linearly many. The vast majority of the computation time is then actually spent on the
   215   multiplication, since the output number is exponential in the input number.
   216 \<close>
   217 
   218 lemma fib_rec_odd:
   219   fixes \<phi> \<psi> :: real
   220   defines "\<phi> \<equiv> (1 + sqrt 5) / 2"
   221     and "\<psi> \<equiv> (1 - sqrt 5) / 2"
   222   shows "fib (Suc (2 * n)) = fib n^2 + fib (Suc n)^2"
   223 proof -
   224   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"
   225     by (simp add: fib_closed_form[folded \<phi>_def \<psi>_def] field_simps power2_eq_square)
   226   also
   227   let ?A = "\<phi>^(2 * n) + \<psi>^(2 * n) - 2*(\<phi> * \<psi>)^n + \<phi>^(2 * n + 2) + \<psi>^(2 * n + 2) - 2*(\<phi> * \<psi>)^(n + 1)"
   228   have "(\<phi> ^ n - \<psi> ^ n)\<^sup>2 + (\<phi> * \<phi> ^ n - \<psi> * \<psi> ^ n)\<^sup>2 = ?A"
   229     by (simp add: power2_eq_square algebra_simps power_mult power_mult_distrib)
   230   also have "\<phi> * \<psi> = -1"
   231     by (simp add: \<phi>_def \<psi>_def field_simps)
   232   then have "?A = \<phi>^(2 * n + 1) * (\<phi> + inverse \<phi>) + \<psi>^(2 * n + 1) * (\<psi> + inverse \<psi>)"
   233     by (auto simp: field_simps power2_eq_square)
   234   also have "1 + sqrt 5 > 0"
   235     by (auto intro: add_pos_pos)
   236   then have "\<phi> + inverse \<phi> = sqrt 5"
   237     by (simp add: \<phi>_def field_simps)
   238   also have "\<psi> + inverse \<psi> = -sqrt 5"
   239     by (simp add: \<psi>_def field_simps)
   240   also have "(\<phi> ^ (2 * n + 1) * sqrt 5 + \<psi> ^ (2 * n + 1) * - sqrt 5) / 5 =
   241     (\<phi> ^ (2 * n + 1) - \<psi> ^ (2 * n + 1)) * (sqrt 5 / 5)"
   242     by (simp add: field_simps)
   243   also have "sqrt 5 / 5 = inverse (sqrt 5)"
   244     by (simp add: field_simps)
   245   also have "(\<phi> ^ (2 * n + 1) - \<psi> ^ (2 * n + 1)) * \<dots> = of_nat (fib (Suc (2 * n)))"
   246     by (simp add: fib_closed_form[folded \<phi>_def \<psi>_def] divide_inverse)
   247   finally show ?thesis
   248     by (simp only: of_nat_eq_iff)
   249 qed
   250 
   251 lemma fib_rec_even: "fib (2 * n) = (fib (n - 1) + fib (n + 1)) * fib n"
   252 proof (induct n)
   253   case 0
   254   then show ?case by simp
   255 next
   256   case (Suc n)
   257   let ?rfib = "\<lambda>x. real (fib x)"
   258   have "2 * (Suc n) = Suc (Suc (2 * n))" by simp
   259   also have "real (fib \<dots>) = ?rfib n^2 + ?rfib (Suc n)^2 + (?rfib (n - 1) + ?rfib (n + 1)) * ?rfib n"
   260     by (simp add: fib_rec_odd Suc)
   261   also have "(?rfib (n - 1) + ?rfib (n + 1)) * ?rfib n = (2 * ?rfib (n + 1) - ?rfib n) * ?rfib n"
   262     by (cases n) simp_all
   263   also have "?rfib n^2 + ?rfib (Suc n)^2 + \<dots> = (?rfib (Suc n) + 2 * ?rfib n) * ?rfib (Suc n)"
   264     by (simp add: algebra_simps power2_eq_square)
   265   also have "\<dots> = real ((fib (Suc n - 1) + fib (Suc n + 1)) * fib (Suc n))" by simp
   266   finally show ?case by (simp only: of_nat_eq_iff)
   267 qed
   268 
   269 lemma fib_rec_even': "fib (2 * n) = (2 * fib (n - 1) + fib n) * fib n"
   270   by (subst fib_rec_even, cases n) simp_all
   271 
   272 lemma fib_rec:
   273   "fib n =
   274     (if n = 0 then 0 else if n = 1 then 1
   275      else if even n then let n' = n div 2; fn = fib n' in (2 * fib (n' - 1) + fn) * fn
   276      else let n' = n div 2 in fib n' ^ 2 + fib (Suc n') ^ 2)"
   277   by (auto elim: evenE oddE simp: fib_rec_odd fib_rec_even' Let_def)
   278 
   279 
   280 subsection \<open>Fibonacci and Binomial Coefficients\<close>
   281 
   282 lemma sum_drop_zero: "(\<Sum>k = 0..Suc n. if 0<k then (f (k - 1)) else 0) = (\<Sum>j = 0..n. f j)"
   283   by (induct n) auto
   284 
   285 lemma sum_choose_drop_zero:
   286   "(\<Sum>k = 0..Suc n. if k = 0 then 0 else (Suc n - k) choose (k - 1)) =
   287     (\<Sum>j = 0..n. (n-j) choose j)"
   288   by (rule trans [OF sum.cong sum_drop_zero]) auto
   289 
   290 lemma ne_diagonal_fib: "(\<Sum>k = 0..n. (n-k) choose k) = fib (Suc n)"
   291 proof (induct n rule: fib.induct)
   292   case 1
   293   show ?case by simp
   294 next
   295   case 2
   296   show ?case by simp
   297 next
   298   case (3 n)
   299   have "(\<Sum>k = 0..Suc n. Suc (Suc n) - k choose k) =
   300      (\<Sum>k = 0..Suc n. (Suc n - k choose k) + (if k = 0 then 0 else (Suc n - k choose (k - 1))))"
   301     by (rule sum.cong) (simp_all add: choose_reduce_nat)
   302   also have "\<dots> =
   303     (\<Sum>k = 0..Suc n. Suc n - k choose k) +
   304     (\<Sum>k = 0..Suc n. if k=0 then 0 else (Suc n - k choose (k - 1)))"
   305     by (simp add: sum.distrib)
   306   also have "\<dots> = (\<Sum>k = 0..Suc n. Suc n - k choose k) + (\<Sum>j = 0..n. n - j choose j)"
   307     by (metis sum_choose_drop_zero)
   308   finally show ?case using 3
   309     by simp
   310 qed
   311 
   312 end