Theory Echelon_Form

(*  
    Title:      Echelon_Form.thy
    Author:     Jose Divasón <jose.divasonm at unirioja.es>
    Author:     Jesús Aransay <jesus-maria.aransay at unirioja.es>
*)

section‹Echelon Form›

theory Echelon_Form
imports 
  Rings2
  Gauss_Jordan.Determinants2
  Cayley_Hamilton_Compatible
begin


subsection‹Definition of Echelon Form›

text‹Echelon form up to column k (NOT INCLUDED).›

definition 
  echelon_form_upt_k :: "'a::{bezout_ring}^'cols::{mod_type}^'rows::{finite, ord}  nat  bool" 
  where 
  "echelon_form_upt_k A k = (
    (i. is_zero_row_upt_k i k A 
           ¬ (j. j>i  ¬ is_zero_row_upt_k j k A)) 
      
    (i j. i<j  ¬ (is_zero_row_upt_k i k A)  ¬ (is_zero_row_upt_k j k A) 
           ((LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0))))"
  
definition "echelon_form A = echelon_form_upt_k A (ncols A)"

text‹Some properties of matrices in echelon form.›

lemma echelon_form_upt_k_intro:
  assumes "(i. is_zero_row_upt_k i k A  ¬ (j. j>i  ¬ is_zero_row_upt_k j k A))"
  and "(i j. i<j  ¬ (is_zero_row_upt_k i k A)  ¬ (is_zero_row_upt_k j k A) 
   ((LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)))"
  shows "echelon_form_upt_k A k" using assms unfolding echelon_form_upt_k_def by fast

lemma echelon_form_upt_k_condition1:
  assumes "echelon_form_upt_k A k" "is_zero_row_upt_k i k A"
  shows "¬ (j. j>i  ¬ is_zero_row_upt_k j k A)"
  using assms unfolding echelon_form_upt_k_def by auto

lemma echelon_form_upt_k_condition1':
  assumes "echelon_form_upt_k A k" "is_zero_row_upt_k i k A" and "i<j"
  shows "is_zero_row_upt_k j k A"
  using assms unfolding echelon_form_upt_k_def by auto

lemma echelon_form_upt_k_condition2:
  assumes "echelon_form_upt_k A k" "i<j"
  and  "¬ (is_zero_row_upt_k i k A)" "¬ (is_zero_row_upt_k j k A)"
  shows "(LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)"
  using assms unfolding echelon_form_upt_k_def by auto

lemma echelon_form_upt_k_if_equal:
  assumes e: "echelon_form_upt_k A k"
  and eq: "a. b<from_nat k. A $ a $ b = B $ a $ b"
  and k: "k < ncols A"
  shows "echelon_form_upt_k B k"
  unfolding echelon_form_upt_k_def
proof (auto)
  fix i j assume zero_iB: "is_zero_row_upt_k i k B" and ij: "i < j" 
  have zero_iA: "is_zero_row_upt_k i k A" 
  proof (unfold is_zero_row_upt_k_def, clarify)
    fix ja::'b assume ja_k: "to_nat ja < k"
    have ja_k2: "ja < from_nat k" 
      by (metis (full_types) add_to_nat_def k from_nat_mono 
        ja_k monoid_add_class.add.right_neutral ncols_def to_nat_0)
    have "A $ i $ ja = B $ i $ ja" using eq ja_k2 by auto
    also have "... = 0" using zero_iB ja_k unfolding is_zero_row_upt_k_def by simp    
    finally show "A $ i $ ja = 0" .
  qed
  hence zero_jA: "is_zero_row_upt_k j k A" by (metis e echelon_form_upt_k_condition1 ij)
  show "is_zero_row_upt_k j k B"
  proof (unfold is_zero_row_upt_k_def, clarify)
    fix ja::'b assume ja_k: "to_nat ja < k"
    have ja_k2: "ja < from_nat k" 
      by (metis (full_types) add_to_nat_def k from_nat_mono 
        ja_k monoid_add_class.add.right_neutral ncols_def to_nat_0)
    have "B $ j $ ja = A $ j $ ja" using eq ja_k2 by auto
    also have "... = 0" using zero_jA ja_k unfolding is_zero_row_upt_k_def by simp    
    finally show "B $ j $ ja = 0" .
  qed
next
  fix i j
  assume ij: "i < j"
    and not_zero_iB: "¬ is_zero_row_upt_k i k B" 
    and not_zero_jB: "¬ is_zero_row_upt_k j k B"
  obtain a where Bia: "B $ i $ a  0" and ak: "a<from_nat k" 
    using not_zero_iB k unfolding is_zero_row_upt_k_def ncols_def
    by (metis add_to_nat_def from_nat_mono monoid_add_class.add.right_neutral to_nat_0)
  have Aia: "A $ i $ a  0" by (metis ak Bia eq)
  obtain b where Bjb: "B $ j $ b  0" and bk: "b<from_nat k" 
    using not_zero_jB k unfolding is_zero_row_upt_k_def ncols_def
    by (metis add_to_nat_def from_nat_mono monoid_add_class.add.right_neutral to_nat_0)
  have Ajb: "A $ j $ b  0" by (metis bk Bjb eq)
  have not_zero_iA: "¬ is_zero_row_upt_k i k A" 
    by (metis (full_types) Aia ak is_zero_row_upt_k_def to_nat_le)
  have not_zero_jA: "¬ is_zero_row_upt_k j k A" 
    by (metis (full_types) Ajb bk is_zero_row_upt_k_def to_nat_le)
  have "(LEAST n. A $ i $ n  0) = (LEAST n. B $ i $ n  0)"
  proof (rule Least_equality)
    have "(LEAST n. B $ i $ n  0)  a" by (rule Least_le, simp add: Bia)
    hence least_bi_less: "(LEAST n. B $ i $ n  0) < from_nat k" using ak by simp
    thus "A $ i $ (LEAST n. B $ i $ n  0)  0" 
      by (metis (mono_tags, lifting) LeastI  eq is_zero_row_upt_k_def not_zero_iB)
    fix y assume "A $ i $ y  0"
    thus "(LEAST n. B $ i $ n  0)  y"
      by (metis (mono_tags, lifting) Least_le  dual_order.strict_trans2 eq least_bi_less linear) 
  qed
  moreover have "(LEAST n. A $ j $ n  0) = (LEAST n. B $ j $ n  0)"
  proof (rule Least_equality)
    have "(LEAST n. B $ j $ n  0)  b" by (rule Least_le, simp add: Bjb)
    hence least_bi_less: "(LEAST n. B $ j $ n  0) < from_nat k" using bk by simp
    thus "A $ j $ (LEAST n. B $ j $ n  0)  0" 
      by (metis (mono_tags, lifting) LeastI eq is_zero_row_upt_k_def not_zero_jB)
    fix y assume "A $ j $ y  0"
    thus "(LEAST n. B $ j $ n  0)  y"
      by (metis (mono_tags, lifting) Least_le  dual_order.strict_trans2 eq least_bi_less linear) 
  qed
  moreover have "(LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)" 
    by (rule echelon_form_upt_k_condition2[OF e ij not_zero_iA not_zero_jA])  
  ultimately show "(LEAST n. B $ i $ n  0) < (LEAST n. B $ j $ n  0)" by auto
qed

lemma echelon_form_upt_k_0: "echelon_form_upt_k A 0"
  unfolding echelon_form_upt_k_def is_zero_row_upt_k_def by auto

lemma echelon_form_condition1:
  assumes r: "echelon_form A"
  shows "(i. is_zero_row i A  ¬ (j. j>i  ¬ is_zero_row j A))" 
  using r unfolding echelon_form_def 
  by (metis echelon_form_upt_k_condition1' is_zero_row_def)

lemma echelon_form_condition2:
  assumes r: "echelon_form A"
  shows "(i. i<j  ¬ (is_zero_row i A)  ¬ (is_zero_row j A) 
   ((LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)))"
  using r unfolding echelon_form_def 
  by (metis echelon_form_upt_k_condition2 is_zero_row_def)


lemma echelon_form_condition2_explicit:
  assumes rref_A: "echelon_form A"
  and i_le: "i < j"
  and "¬ is_zero_row i A" and "¬ is_zero_row j A"
  shows "(LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)"
  using echelon_form_condition2 assms by blast

lemma echelon_form_intro:
  assumes 1: "(i. is_zero_row i A  ¬ (j. j>i  ¬ is_zero_row j A))"
  and 2: "(i j. i<j  ¬ (is_zero_row i A)  ¬ (is_zero_row j A) 
   ((LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)))"
  shows "echelon_form A"  
proof (unfold echelon_form_def, rule echelon_form_upt_k_intro, auto) 
  fix i j assume "is_zero_row_upt_k i (ncols A) A" and "i < j"
  thus "is_zero_row_upt_k j (ncols A) A"
    using 1 is_zero_row_imp_is_zero_row_upt by (metis is_zero_row_def)
next
  fix i j
  assume "i < j" and "¬ is_zero_row_upt_k i (ncols A) A" and "¬ is_zero_row_upt_k j (ncols A) A"
  thus "(LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)"
    using 2 by (metis is_zero_row_imp_is_zero_row_upt)
qed

lemma echelon_form_implies_echelon_form_upt:
  fixes A::"'a::{bezout_ring}^'cols::{mod_type}^'rows::{mod_type}"
  assumes rref: "echelon_form A"
  shows "echelon_form_upt_k A k"
proof (rule echelon_form_upt_k_intro)
  show "i. is_zero_row_upt_k i k A  ¬ (j>i. ¬ is_zero_row_upt_k j k A)"
  proof (auto, rule ccontr)
    fix i j assume zero_i_k: "is_zero_row_upt_k i k A" and i_less_j: "i < j"
      and not_zero_j_k:"¬ is_zero_row_upt_k j k A"
    have not_zero_j: "¬ is_zero_row j A"
      using is_zero_row_imp_is_zero_row_upt not_zero_j_k by blast
    hence not_zero_i: "¬ is_zero_row i A"
      using echelon_form_condition1[OF rref] i_less_j by blast
    have Least_less: "(LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)"
      by (rule echelon_form_condition2_explicit[OF rref i_less_j not_zero_i not_zero_j])
    moreover have "(LEAST n. A $ j $ n  0) < (LEAST n. A $ i $ n  0)"
    proof (rule LeastI2_ex)   
      show "a. A $ i $ a  0"
        using not_zero_i unfolding is_zero_row_def is_zero_row_upt_k_def by fast
      fix x assume Aix_not_0: "A $ i $ x  0"
      have k_less_x: "k  to_nat x" 
        using zero_i_k Aix_not_0 unfolding is_zero_row_upt_k_def by force
      hence k_less_ncols: "k < ncols A"
        unfolding ncols_def using to_nat_less_card[of x] by simp
      obtain s where Ajs_not_zero: "A $ j $ s  0" and s_less_k: "to_nat s < k"
        using not_zero_j_k unfolding is_zero_row_upt_k_def by blast
      have "(LEAST n. A $ j $ n  0)  s" using Ajs_not_zero Least_le by fast
      also have "... = from_nat (to_nat s)" unfolding from_nat_to_nat_id ..
      also have "... < from_nat k"
        by (rule from_nat_mono[OF s_less_k k_less_ncols[unfolded ncols_def]])
      also have "...  x" using k_less_x leD not_le_imp_less to_nat_le by fast
      finally show "(LEAST n. A $ j $ n  0) < x" .
    qed
    ultimately show False by fastforce
  qed
  show "i j. i < j  ¬ is_zero_row_upt_k i k A  ¬ is_zero_row_upt_k j k A 
     (LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)"
    using echelon_form_condition2[OF rref] is_zero_row_imp_is_zero_row_upt by blast
qed

lemma upper_triangular_upt_k_def':
  assumes "i j. to_nat j  k  A $ i $ j  0  ji"
  shows "upper_triangular_upt_k A k"
  using assms 
  unfolding upper_triangular_upt_k_def
  by (metis leD linear) 


lemma echelon_form_imp_upper_triagular_upt:
  fixes A::"'a::{bezout_ring}^'n::{mod_type}^'n::{mod_type}"
  assumes "echelon_form A"
  shows "upper_triangular_upt_k A k"
proof (induct k)
  case 0
  show ?case unfolding upper_triangular_upt_k_def by simp
next
  case (Suc k)
  show ?case
    unfolding upper_triangular_upt_k_def 
  proof (clarify)
    fix i j::'n assume j_less_i: "j < i" and j_less_suc_k: "to_nat j < Suc k"
    show "A $ i $ j = 0"
    proof (cases "to_nat j < k")
      case True
      thus ?thesis 
        using Suc.hyps
        unfolding upper_triangular_upt_k_def using j_less_i True by auto
    next
      case False
      hence j_eq_k: "to_nat j = k" using j_less_suc_k by simp
      hence j_eq_k2: "from_nat k = j" by (metis from_nat_to_nat_id)
      have rref_suc: "echelon_form_upt_k A (Suc k)"
        by (metis assms echelon_form_implies_echelon_form_upt)   
      have zero_j_k: "is_zero_row_upt_k j k A" 
        unfolding is_zero_row_upt_k_def
        by (metis (opaque_lifting, mono_tags) Suc.hyps leD le_less_linear 
          j_eq_k to_nat_mono' upper_triangular_upt_k_def)
      hence zero_i_k: "is_zero_row_upt_k i k A" 
        by (metis (poly_guards_query) assms echelon_form_implies_echelon_form_upt 
          echelon_form_upt_k_condition1' j_less_i)
      show ?thesis
      proof (cases "A $ j $ j = 0")
        case True
        have "is_zero_row_upt_k j (Suc k) A"
          by (rule is_zero_row_upt_k_suc[OF zero_j_k], simp add: True j_eq_k2)
        hence "is_zero_row_upt_k i (Suc k) A" 
          by (metis echelon_form_upt_k_condition1' j_less_i rref_suc)
        thus ?thesis by (metis is_zero_row_upt_k_def j_eq_k lessI)
      next
        case False note Ajj_not_zero=False
        hence not_zero_j:"¬ is_zero_row_upt_k j (Suc k) A" 
          by (metis is_zero_row_upt_k_def j_eq_k lessI)
        show ?thesis
        proof (rule ccontr)
          assume Aij_not_zero: "A $ i $ j  0"
          hence not_zero_i: "¬ is_zero_row_upt_k i (Suc k) A" 
            by (metis is_zero_row_upt_k_def j_eq_k lessI)
          have Least_eq: "(LEAST n. A $ i $ n  0) = from_nat k"
          proof (rule Least_equality)
            show "A $ i $ from_nat k  0" using Aij_not_zero j_eq_k2 by simp
            show "y. A $ i $ y  0  from_nat k  y" 
              by (metis (full_types) is_zero_row_upt_k_def not_le_imp_less to_nat_le zero_i_k)
          qed
          moreover have Least_eq2: "(LEAST n. A $ j $ n  0) = from_nat k"
          proof (rule Least_equality)
            show "A $ j $ from_nat k  0" using Ajj_not_zero j_eq_k2 by simp
            show "y. A $ j $ y  0  from_nat k  y" 
              by (metis (full_types) is_zero_row_upt_k_def not_le_imp_less to_nat_le zero_j_k)
          qed
          ultimately show False 
            using echelon_form_upt_k_condition2[OF rref_suc j_less_i not_zero_j not_zero_i]
            by simp
        qed
      qed
    qed
  qed
qed

text‹A matrix in echelon form is upper triangular.›
lemma echelon_form_imp_upper_triagular:
  fixes A::"'a::{bezout_ring}^'n::{mod_type}^'n::{mod_type}"
  assumes "echelon_form A"
  shows "upper_triangular A"
  using echelon_form_imp_upper_triagular_upt[OF assms] 
  by (metis upper_triangular_upt_imp_upper_triangular)


lemma echelon_form_upt_k_interchange:
  fixes A::"'a::{bezout_ring}^'c::{mod_type}^'b::{mod_type}"
  assumes e: "echelon_form_upt_k A k"
  and zero_ikA: "is_zero_row_upt_k (from_nat i) k A"
  and Amk_not_0: "A $ m $ from_nat k  0"
  and i_le_m: "(from_nat i)m"
  and k: "k<ncols A"
  shows "echelon_form_upt_k (interchange_rows A (from_nat i) 
  (LEAST n. A $ n $ from_nat k  0  (from_nat i)  n)) k"
proof (rule echelon_form_upt_k_if_equal[OF e _ k], auto)
  fix a and b::'c
  assume b: "b < from_nat k"
  let ?least = "(LEAST n. A $ n $ from_nat k  0  (from_nat i)  n)"
  let ?interchange = "(interchange_rows A (from_nat i) ?least)"
  have "(from_nat i)?least" by (metis (mono_tags, lifting) Amk_not_0 LeastI_ex i_le_m)
  hence zero_leastkA: "is_zero_row_upt_k ?least k A"
    using echelon_form_upt_k_condition1[OF e zero_ikA] 
    by (metis (poly_guards_query) dual_order.strict_iff_order zero_ikA)
  show "A $ a $ b = ?interchange $ a $ b"
  proof (cases "a=from_nat i")
    case True
    hence "?interchange $ a $ b = A $ ?least $ b" unfolding interchange_rows_def by auto
    also have "... = 0" using zero_leastkA unfolding is_zero_row_upt_k_def 
      by (metis (mono_tags) b to_nat_le)
    finally have "?interchange $ a $ b = 0" .
    moreover have "A $ a $ b = 0" 
      by (metis True b is_zero_row_upt_k_def to_nat_le zero_ikA)
    ultimately show ?thesis by simp
  next
    case False note a_not_i=False
    show ?thesis
    proof (cases "a=?least")
      case True
      hence "?interchange $ a $ b = A $ (from_nat i) $ b" unfolding interchange_rows_def by auto
      also have "... = 0" using zero_ikA  unfolding is_zero_row_upt_k_def 
        by (metis (poly_guards_query) b to_nat_le)
      finally have "?interchange $ a $ b = 0" .
      moreover have "A $ a $ b = 0" by (metis True b is_zero_row_upt_k_def to_nat_le zero_leastkA)
      ultimately show ?thesis by simp
    next
      case False
      thus ?thesis using a_not_i unfolding interchange_rows_def by auto
    qed
  qed
qed


text‹There are similar theorems to the following ones in the Gauss-Jordan developments, but
for matrices in reduced row echelon form. It is possible to prove that reduced row echelon form 
implies echelon form. Then the theorems in the Gauss-Jordan development could be 
obtained with ease.›

lemma greatest_less_zero_row:
  fixes A::"'a::{bezout_ring}^'cols::{mod_type}^'rows::{finite, wellorder}"
  assumes r: "echelon_form_upt_k A k"
  and zero_i: "is_zero_row_upt_k i k A"
  and not_all_zero: "¬ (a. is_zero_row_upt_k a k A)"
  shows "(GREATEST m. ¬ is_zero_row_upt_k m k A) < i"
proof (rule ccontr)
  assume not_less_i: "¬ (GREATEST m. ¬ is_zero_row_upt_k m k A) < i"
  have i_less_greatest: "i < (GREATEST m. ¬ is_zero_row_upt_k m k A)"
    by (metis not_less_i neq_iff GreatestI not_all_zero zero_i)
  have "is_zero_row_upt_k (GREATEST m. ¬ is_zero_row_upt_k m k A) k A"
    using r zero_i i_less_greatest unfolding echelon_form_upt_k_def by blast
  thus False using GreatestI_ex not_all_zero by fast
qed

lemma greatest_ge_nonzero_row':
  fixes A::"'a::{bezout_ring}^'cols::{mod_type}^'rows::{mod_type}"
  assumes r: "echelon_form_upt_k A k"
  and i: "i  (GREATEST m. ¬ is_zero_row_upt_k m k A)"
  and not_all_zero: "¬ (a. is_zero_row_upt_k a k A)"
  shows "¬ is_zero_row_upt_k i k A"
  using greatest_less_zero_row[OF r] i not_all_zero by fastforce

lemma rref_imp_ef: 
  fixes A::"'a::{bezout_ring}^'cols::{mod_type}^'rows::{mod_type}"
  assumes rref: "reduced_row_echelon_form A"
  shows "echelon_form A"
proof (rule echelon_form_intro)
  show "i. is_zero_row i A  ¬ (j>i. ¬ is_zero_row j A)"
    by (simp add: rref rref_condition1)
  show "i j. i < j  ¬ is_zero_row i A  ¬ is_zero_row j A 
     (LEAST n. A $ i $ n  0) < (LEAST n. A $ j $ n  0)"
    by (simp add: rref_condition3_equiv rref)
qed

subsection‹Computing the echelon form of a matrix›

subsubsection‹Demonstration over principal ideal rings›
       
text‹Important remark:

We want to prove that there exist the echelon form of any matrix whose elements belong to a bezout 
domain. In addition, we want to compute the echelon form, so we will need computable gcd 
and bezout operations which is possible over euclidean domains. 
Our approach consists of demonstrating the correctness over bezout domains
and executing over euclidean domains.

To do that, we have studied several options:

\begin{enumerate}
  \item We could define a gcd in bezout rings (bezout_ring_gcd›) as follows:
  gcd_bezout_ring a b = (SOME d. d dvd a ∧ d dvd b ∧ (∀d'. d' dvd a ∧ d' dvd b ⟶ d' dvd d))›

  And then define an algorithm that computes the Echelon Form using such a definition to the gcd.
  This would allow us to prove the correctness over bezout rings, but we would not be able
  to execute over euclidean rings because it is not possible to demonstrate a (code) lemma 
  stating that (gcd_bezout_ring a b) = gcd_eucl a b› (the gcd is not unique over 
  bezout rings and GCD rings).
  
  \item Create a bezout_ring_norm› class and define a gcd normalized over bezout rings:
  definition gcd_bezout_ring_norm a b = gcd_bezout_ring a b div normalisation_factor (gcd_bezout_ring a b)›

  Then, one could demonstrate a (code) lemma stating that: (gcd_bezout_ring_norm a b) 
  = gcd_eucl a b›
  This allows us to execute the gcd function, but with bezout it is not possible.

  \item The third option (and the chosen one) consists of defining the algorithm over bezout domains 
  and parametrizing the algorithm by a bezout› operation which must satisfy 
  suitable properties (i.e @{term "is_bezout_ext bezout"}). Then we can prove the correctness over 
  bezout domains and we will execute over euclidean domains, since we can prove that the 
  operation  @{term "euclid_ext2"} is an executable operation which satisfies 
  @{term "is_bezout_ext euclid_ext2"}.
\end{enumerate}
›

subsubsection‹Definition of the algorithm›

context bezout_ring
begin

definition 
  bezout_matrix :: "'a^'cols^'rows  'rows  'rows  'cols 
                     ('a  'a  ('a × 'a × 'a × 'a × 'a))  'a^'rows^'rows"
  where 
  "bezout_matrix A a b j bezout = (χ x y. 
      (let 
        (p, q, u, v, d) = bezout (A $ a $ j) (A $ b $ j) 
       in
         if x = a  y = a then p else
         if x = a  y = b then q else
         if x = b  y = a then u else
         if x = b  y = b then v else
         if x = y then 1 else 0))"

end

primrec
  bezout_iterate :: "'a::{bezout_ring}^'cols^'rows::{mod_type} 
                      nat  'rows::{mod_type} 
                      'cols  ('a 'a  ('a × 'a × 'a × 'a × 'a))  'a^'cols^'rows::{mod_type}"
where "bezout_iterate A 0 i j bezout = A"
    | "bezout_iterate A (Suc n) i j bezout = 
        (if (Suc n)  to_nat i then A else 
              bezout_iterate (bezout_matrix A i (from_nat (Suc n)) j bezout ** A) n i j bezout)"

text‹If every element in column @{term "k::nat"} over index @{term "i::nat"} are equal to zero,
      the same input is returned. If every element over @{term "i::nat"} 
      is equal to zero, except the pivot, the algorithm does nothing, but pivot @{term "i::nat"}
      is increased in a unit. Finally, if there is a position @{term "n::nat"} 
      whose coefficient is different from zero, its row is interchanged with row 
      @{term "i::nat"} and the bezout coefficients are used to produce a zero in its position.›


definition 
  "echelon_form_of_column_k bezout A' k = 
    (let (A, i) = A' 
     in if (mfrom_nat i. A $ m $ from_nat k = 0)  (i = nrows A) then (A, i) else 
        if (m>from_nat i. A $ m $ from_nat k = 0) then (A, i + 1) else
            let n = (LEAST n. A $ n $ from_nat k  0  from_nat i  n); 
                interchange_A = interchange_rows A (from_nat i) n
           in
            (bezout_iterate (interchange_A) (nrows A - 1) (from_nat i) (from_nat k) bezout, i + 1))"

definition "echelon_form_of_upt_k A k bezout = (fst (foldl (echelon_form_of_column_k bezout) (A,0) [0..<Suc k]))"
definition "echelon_form_of A bezout = echelon_form_of_upt_k A (ncols A - 1) bezout"

subsubsection‹The executable definition:›

context euclidean_space
begin

definition [code_unfold]: "echelon_form_of_euclidean A = echelon_form_of A euclid_ext2"

end

subsubsection‹Properties of the bezout matrix›

lemma bezout_matrix_works1:
  assumes ib: "is_bezout_ext bezout"
  and a_not_b: "a  b"
  shows "(bezout_matrix A a b j bezout ** A) $ a $ j = snd (snd (snd (snd (bezout (A $ a $ j) (A $ b $ j)))))"
proof (unfold matrix_matrix_mult_def bezout_matrix_def Let_def, simp)
  let ?a = "(A $ a $ j)"
  let ?b = "(A $ b $ j)"
  let ?z = "bezout (A $ a $ j) (A $ b $ j)"
  obtain p q u v d where bz: "(p, q, u, v, d) = ?z" by (cases ?z, auto)
  from ib have foo: "(a b. let (p, q, u, v, gcd_a_b) = bezout a b
           in p * a + q * b = gcd_a_b 
              gcd_a_b dvd a 
              gcd_a_b dvd b  (d'. d' dvd a  d' dvd b  d' dvd gcd_a_b)  gcd_a_b * u = - b  gcd_a_b * v = a)"
           using is_bezout_ext_def [of bezout] by simp
  have foo: "p * ?a + q * ?b = d  d dvd ?a 
            d dvd ?b  (d'. d' dvd ?a  d' dvd ?b  d' dvd d)  d * u = - ?b  d * v = ?a" 
            using ib using is_bezout_ext_def using bz [symmetric]
            using foo [of ?a ?b] by fastforce
  have pa_bq_d: "p * ?a + ?b * q = d" using foo by (auto simp add: mult.commute)
  define f where "f k = (if k = a then p
              else if a = a  k = b then q
              else if a = b  k = a then u
              else if a = b  k = b then v
              else if a = k then 1 else 0) * A $ k $ j" for k
  have UNIV_rw: "UNIV = insert b (insert a (UNIV - {a} - {b}))" by auto
  have sum_rw: "sum f (insert a (UNIV - {a} - {b})) = f a + sum f (UNIV - {a} - {b})"
    by (rule sum.insert, auto)
  have sum0: "sum f (UNIV - {a} - {b}) = 0" by (rule sum.neutral, simp add: f_def)
  have "(kUNIV.
       (case bezout (A $ a $ j) (A $ b $ j) of
        (p, q, u, v, d) 
          if k = a then p
          else if a = a  k = b then q
               else if a = b  k = a then u else if a = b  k = b then v else if a = k then 1 else 0) *
       A $ k $ j) = (kUNIV.
       (if k = a then p
          else if a = a  k = b then q
               else if a = b  k = a then u else if a = b  k = b then v else if a = k then 1 else 0) *
       A $ k $ j)" unfolding bz [symmetric] by auto
  also have "... = sum f UNIV" unfolding f_def ..
  also have "sum f UNIV = sum f (insert b (insert a (UNIV - {a} - {b})))" using UNIV_rw by simp
  also have "... = f b + sum f (insert a (UNIV - {a} - {b}))"
    by (rule sum.insert, auto, metis a_not_b)
  also have "... = f b + f a" unfolding sum_rw sum0 by simp
  also have "... = d"
    unfolding f_def using a_not_b bz [symmetric] by (auto, metis add.commute mult.commute pa_bq_d)
  also have "... = snd (snd (snd (snd (bezout (A $ a $ j) (A $ b $ j)))))"
    using bz by (metis snd_conv)
  finally show "(kUNIV.
       (case bezout (A $ a $ j) (A $ b $ j) of
        (p, q, u, v, d) 
          if k = a then p
          else if a = a  k = b then q
               else if a = b  k = a then u else if a = b  k = b then v else if a = k then 1 else 0) *
       A $ k $ j) =
    snd (snd (snd (snd (bezout (A $ a $ j) (A $ b $ j)))))" unfolding f_def by simp
qed

lemma bezout_matrix_not_zero:
  assumes ib: "is_bezout_ext bezout"
  and a_not_b: "a  b"
  and Aaj: "A $ a $ j  0"
  shows "(bezout_matrix A a b j bezout ** A) $ a $ j  0"
proof -
  have "(bezout_matrix A a b j bezout ** A) $ a $ j = snd (snd (snd(snd (bezout (A $ a $ j) (A $ b $ j)))))"
    using bezout_matrix_works1[OF ib a_not_b] .
  also have "... = (λa b. (case bezout a b of (_, _,_ ,_,gcd')  (gcd'))) (A $ a $ j) (A $ b $ j)"
    by (simp add: split_beta)
  also have "...  0" using gcd'_zero[OF is_gcd_is_bezout_ext[OF ib]] Aaj by simp
  finally show ?thesis .
qed

lemma ua_vb_0:
  fixes a::"'a::bezout_domain"
  assumes ib: "is_bezout_ext bezout" and nz: "snd (snd (snd (snd (bezout a b))))  0"
  shows "fst (snd (snd (bezout a b))) * a + fst (snd (snd (snd (bezout a b)))) * b = 0"
proof-
  obtain p q u v d where bz: "(p, q, u, v, d) = bezout a b" by (cases "bezout a b", auto)
  from ib have foo: "(a b. let (p, q, u, v, gcd_a_b) = bezout a b
           in p * a + q * b = gcd_a_b 
              gcd_a_b dvd a 
              gcd_a_b dvd b  (d'. d' dvd a  d' dvd b  d' dvd gcd_a_b)  gcd_a_b * u = - b  gcd_a_b * v = a)"
           using is_bezout_ext_def [of bezout] by simp
  have "p * a + q * b = d  d dvd a 
            d dvd b  (d'. d' dvd a  d' dvd b  d' dvd d)  d * u = - b  d * v = a" 
            using foo [of a b] using bz by fastforce
  hence dub: "d * u = - b" and dva: "d * v = a" by (simp_all)
  hence "d * u * a + d * v * b = 0" 
    using eq_neg_iff_add_eq_0 mult.commute mult_minus_left by auto
  hence "u * a + v * b = 0"
    by (metis (no_types, lifting) dub dva minus_minus mult_minus_left 
          neg_eq_iff_add_eq_0 semiring_normalization_rules(18) semiring_normalization_rules(7))
  thus ?thesis using bz [symmetric]
  by simp 
qed

lemma bezout_matrix_works2:
  fixes A::"'a::bezout_domain^'cols^'rows"
  assumes ib: "is_bezout_ext bezout"
  and a_not_b: "a  b"
  and not_0: "A $ a $ j  0  A $ b $ j  0"
  shows "(bezout_matrix A a b j bezout ** A) $ b $ j = 0"
proof (unfold matrix_matrix_mult_def bezout_matrix_def Let_def, auto)
  let ?a = "(A $ a $ j)"
  let ?b = "(A $ b $ j)"
  let ?z = "bezout (A $ a $ j) (A $ b $ j)"
  from ib have foo: "(a b. let (p, q, u, v, gcd_a_b) = bezout a b
           in p * a + q * b = gcd_a_b 
              gcd_a_b dvd a 
              gcd_a_b dvd b  (d'. d' dvd a  d' dvd b  d' dvd gcd_a_b)  gcd_a_b * u = - b  gcd_a_b * v = a)"
           using is_bezout_ext_def [of bezout] by simp
  obtain p q u v d where bz: "(p, q, u, v, d) = ?z" by (cases ?z, auto)
  hence pib: "p * ?a + q * ?b = d  d dvd ?a 
            d dvd ?b  (d'. d' dvd ?a  d' dvd ?b  d' dvd d)  d * u = - ?b  d * v = ?a" 
            using foo [of ?a ?b] by fastforce
  hence pa_bq_d: "p * ?a + ?b * q = d" by (simp add: mult.commute)
  have d_dvd_a: "d dvd ?a" using pib by auto
  have d_dvd_b: "d dvd -?b" using pib by auto
  have pa_bq_d: "p * ?a + ?b * q = d" using pa_bq_d by (simp add: mult.commute)
  define f where "f k = (if b = a  k = a then p
                 else if b = a  k = b then q
                      else if b = b  k = a then u
                           else if b = b  k = b then v else if b = k then 1 else 0) *
                A $ k $ j" for k
  have UNIV_rw: "UNIV = insert b (insert a (UNIV - {a} - {b}))" by auto
  have sum_rw: "sum f (insert a (UNIV - {a} - {b})) = f a + sum f (UNIV - {a} - {b})"
    by (rule sum.insert, auto)
  have sum0: "sum f (UNIV - {a} - {b}) = 0" by (rule sum.neutral, simp add: f_def)
  have "(kUNIV.
       (case bezout (A $ a $ j) (A $ b $ j) of
        (p, q, u, v, d) 
          if b = a  k = a then p
          else if b = a  k = b then q
               else if b = b  k = a then u else if b = b  k = b then v else if b = k then 1 else 0) *
       A $ k $ j) = sum f UNIV" unfolding f_def bz [symmetric] by simp
  also have "sum f UNIV = sum f (insert b (insert a (UNIV - {a} - {b})))" using UNIV_rw by simp
  also have "... = f b + sum f (insert a (UNIV - {a} - {b}))"
    by (rule sum.insert, auto, metis a_not_b)
  also have "... = f b + f a" unfolding sum_rw sum0 by simp
  also have "... = v * ?b + u * ?a" unfolding f_def using a_not_b by auto
  also have "... = u * ?a + v * ?b" by auto
  also have "... = 0"
    using ua_vb_0 [OF ib] bz
    by (metis fst_conv minus_minus minus_zero mult_eq_0_iff pib snd_conv)
  finally show "(kUNIV.
       (case bezout (A $ a $ j) (A $ b $ j) of
        (p, q, u, v, d) 
          if b = a  k = a then p
          else if b = a  k = b then q
               else if b = b  k = a then u else if b = b  k = b then v else if b = k then 1 else 0) *
       A $ k $ j) =
       0" .
qed

lemma bezout_matrix_preserves_previous_columns:
  assumes ib: "is_bezout_ext bezout"
  and i_not_j: "i  j"
  and Aik: "A $ i $ k  0"
  and b_k: "b<k"
  and i: "is_zero_row_upt_k i (to_nat k) A" and j: "is_zero_row_upt_k j (to_nat k) A"
  shows "(bezout_matrix A i j k bezout ** A) $ a $ b = A $ a $ b"
  unfolding matrix_matrix_mult_def unfolding bezout_matrix_def Let_def 
proof (auto)
  let ?B = "bezout_matrix A i j k bezout"
  let ?i = "(A $ i $ k)"
  let ?j = "(A $ j $ k)"
  let ?z = "bezout (A $ i $ k) (A $ j $ k)"
  from ib have foo: "(a b. let (p, q, u, v, gcd_a_b) = bezout a b
           in p * a + q * b = gcd_a_b 
              gcd_a_b dvd a 
              gcd_a_b dvd b  (d'. d' dvd a  d' dvd b  d' dvd gcd_a_b)  gcd_a_b * u = - b  gcd_a_b * v = a)"
           using is_bezout_ext_def [of bezout] by simp
  obtain p q u v d where bz: "(p, q, u, v, d) = ?z" by (cases ?z, auto)
  have Aib: "A $ i $ b = 0" by (metis b_k i is_zero_row_upt_k_def to_nat_mono)
  have Ajb: "A $ j $ b = 0" by (metis b_k j is_zero_row_upt_k_def to_nat_mono)
  define f where "f ka = (if a = i  ka = i then p
                  else if a = i  ka = j then q
                  else if a = j  ka = i then u
                  else if a = j  ka = j then v else if a = ka then 1 else 0) * A $ ka $ b" for ka
  show "(kaUNIV.
       (case bezout (A $ i $ k) (A $ j $ k) of
        (p, q, u, v, d) 
          if a = i  ka = i then p
          else if a = i  ka = j then q
               else if a = j  ka = i then u else if a = j  ka = j then v else if a = ka then 1 else 0) *
       A $ ka $ b) =
    A $ a $ b"
  proof (cases "a=i")
    case True
    have "(kaUNIV.
       (case bezout (A $ i $ k) (A $ j $ k) of
        (p, q, u, v, d) 
          if a = i  ka = i then p
          else if a = i  ka = j then q
               else if a = j  ka = i then u else if a = j  ka = j then v else if a = ka then 1 else 0) *
       A $ ka $ b) = sum f UNIV" unfolding f_def bz [symmetric] by simp
    also have "sum f UNIV = 0" by (rule sum.neutral, auto simp add: Aib Ajb f_def True i_not_j)
    also have "... = A $ a $ b" unfolding True using Aib by simp
    finally show ?thesis .
  next
    case False note a_not_i=False
    show ?thesis
    proof (cases "a=j")
      case True 
      have "(kaUNIV.
       (case bezout (A $ i $ k) (A $ j $ k) of
        (p, q, u, v, d) 
          if a = i  ka = i then p
          else if a = i  ka = j then q
               else if a = j  ka = i then u else if a = j  ka = j then v else if a = ka then 1 else 0) *
       A $ ka $ b) = sum f UNIV" unfolding f_def bz [symmetric] by simp
    also have "sum f UNIV = 0" by (rule sum.neutral, auto simp add: Aib Ajb f_def True i_not_j)
      also have "... = A $ a $ b" unfolding True using Ajb by simp
      finally show ?thesis .
    next
      case False 
      have UNIV_rw: "UNIV = insert j (insert i (UNIV - {i} - {j}))" by auto
      have UNIV_rw2: "UNIV - {i} - {j} = insert a (UNIV - {i} - {j}-{a})"
        using False a_not_i by auto
      have sum0: "sum f (UNIV - {i} - {j} - {a}) = 0"
        by (rule sum.neutral, simp add: f_def)
      have "(kaUNIV.
       (case bezout (A $ i $ k) (A $ j $ k) of
        (p, q, u, v, d) 
          if a = i  ka = i then p
          else if a = i  ka = j then q
               else if a = j  ka = i then u else if a = j  ka = j then v else if a = ka then 1 else 0) *
       A $ ka $ b) = sum f UNIV" unfolding f_def bz [symmetric] by simp
      also have "sum f UNIV = sum f (insert j (insert i (UNIV - {i} - {j})))"
        using UNIV_rw by simp
      also have "... = f j + sum f (insert i (UNIV - {i} - {j}))"
        by (rule sum.insert, auto, metis i_not_j)
      also have "... = sum f (insert i (UNIV - {i} - {j}))"
        unfolding f_def using False a_not_i by auto
      also have "... = f i + sum f (UNIV - {i} - {j})"  by (rule sum.insert, auto)
      also have "... = sum f (UNIV - {i} - {j})" unfolding f_def using False a_not_i by auto
      also have "... = sum f (insert a (UNIV - {i} - {j} - {a}))" using UNIV_rw2 by simp
      also have "... = f a + sum f (UNIV - {i} - {j} - {a})" by (rule sum.insert, auto)
      also have "... = f a" unfolding sum0 by simp
      also have "... = A $ a $ b" unfolding f_def using False a_not_i by auto
      finally show ?thesis .
    qed
  qed
qed

lemma det_bezout_matrix:
  fixes A::"'a::{bezout_domain}^'cols^'rows::{finite,wellorder}"
  assumes ib: "is_bezout_ext bezout"
  and a_less_b: "a < b"
  and aj: "A $ a $ j  0"
  shows "det (bezout_matrix A a b j bezout) = 1"
proof -
  let ?B = "bezout_matrix A a b j bezout"
  let ?a = "(A $ a $ j)"
  let ?b = "(A $ b $ j)"
  let ?z = "bezout ?a ?b"
  from ib have foo: "(a b. let (p, q, u, v, gcd_a_b) = bezout a b
           in p * a + q * b = gcd_a_b 
              gcd_a_b dvd a 
              gcd_a_b dvd b  (d'. d' dvd a  d' dvd b  d' dvd gcd_a_b)  gcd_a_b * u = - b  gcd_a_b * v = a)"
           using is_bezout_ext_def [of bezout] by simp
  obtain p q u v d where bz: "(p, q, u, v, d) = ?z" by (cases ?z, auto)
  hence pib: "p * ?a + q * ?b = d  d dvd ?a 
            d dvd ?b  (d'. d' dvd ?a  d' dvd ?b  d' dvd d)  d * u = - ?b  d * v = ?a" 
            using foo [of ?a ?b] by fastforce
  hence pa_bq_d: "p * ?a + ?b * q = d" by (simp add: mult.commute)
  have a_not_b: "a  b" using a_less_b by auto
  have d_dvd_a: "d dvd ?a" using pib by auto
  have UNIV_rw: "UNIV = insert b (insert a (UNIV - {a} - {b}))" by auto
  show ?thesis
  proof (cases "p = 0")
    case True note p0=True
    have q_not_0: "q  0"
    proof (rule ccontr, simp)
      assume q: "q = 0"
      have "d = 0" using pib
        by (metis True q add.right_neutral mult.commute mult_zero_right)
      hence "A $ a $ j = 0  A $ b $ j = 0"
        by (metis aj d_dvd_a dvd_0_left_iff)
      thus False using aj by auto
    qed
    have d_not_0: "d  0"
      by (metis aj d_dvd_a dvd_0_left_iff)
    have qb_not_0: "q *(-?b)  0"
      by (metis d_not_0 mult_cancel_left1 neg_equal_0_iff_equal 
          no_zero_divisors p0 pa_bq_d q_not_0 right_minus)
    have "det (interchange_rows ?B a b) = (iUNIV. (interchange_rows ?B a b) $ i $ i)"
    proof (rule det_upperdiagonal)
      fix i ja::'rows assume ja_i: "ja<i"
      show "interchange_rows (bezout_matrix A a b j bezout) a b $ i $ ja = 0"    
        unfolding interchange_rows_def using a_less_b ja_i p0 a_not_b 
        using bz [symmetric]
        unfolding bezout_matrix_def Let_def by auto
    qed
    also have " = -1" 
    proof -
      define f where "f i = interchange_rows (bezout_matrix A a b j bezout) a b $ i $ i" for i
      have prod_rw: "prod f (insert a (UNIV - {a} - {b})) 
        =  f a * prod f (UNIV - {a} - {b})"
        by (rule prod.insert, simp_all)
      have prod1: "prod f (UNIV - {a} - {b}) = 1"
        by (rule prod.neutral) 
            (simp add: f_def interchange_rows_def bezout_matrix_def Let_def)
      have "prod f UNIV = prod f (insert b (insert a (UNIV - {a} - {b})))" 
        using UNIV_rw by simp
      also have "... = f b * prod f (insert a (UNIV - {a} - {b}))"
      proof (rule prod.insert, simp)
        show "b  insert a (UNIV - {a} - {b})" using a_not_b by auto
      qed
      also have "... = f b * f a" unfolding prod_rw prod1 by auto
      also have "... = q * u" 
        using a_not_b 
        using bz [symmetric]
        unfolding f_def interchange_rows_def bezout_matrix_def Let_def by auto
      also have "... = -1"
      proof -
        let ?r = "q * u"
        have du_b: " d * u = -?b" using pib by auto
        hence "q * (-?b) = d * ?r" by (metis mult.left_commute)
        also have "... = (p * ?a + ?b * q) * ?r" unfolding pa_bq_d by auto
        also have "... = ?b * q * ?r" using True by auto
        also have "... = q * (-?b) * (-?r)" by auto
        finally show ?thesis using qb_not_0 
          unfolding mult_cancel_left1 by (metis minus_minus)
      qed
      finally show ?thesis unfolding f_def .
    qed
    finally have det_inter_1: "det (interchange_rows ?B a b) = - 1" .
    have "det (bezout_matrix A a b j bezout) = - 1 * det (interchange_rows ?B a b)"
      unfolding det_interchange_rows using a_not_b by auto
    thus ?thesis unfolding det_inter_1 by simp
  next
    case False
    define mult_b_dp where "mult_b_dp = mult_row ?B b (d * p)"
    define sum_ab where "sum_ab = row_add mult_b_dp b a ?b"
    have "det (sum_ab) = prod (λi. sum_ab $ i $ i) UNIV"
    proof (rule det_upperdiagonal)
      fix i j::'rows 
      assume j_less_i: "j < i"
      have "d * p * u + ?b * p = 0"
        using pib
        by (metis eq_neg_iff_add_eq_0 mult_minus_left semiring_normalization_rules(16))
      thus "sum_ab $ i $ j = 0"
        unfolding sum_ab_def mult_b_dp_def unfolding row_add_def
        unfolding mult_row_def bezout_matrix_def
        using a_not_b j_less_i a_less_b 
        unfolding bz [symmetric] by auto
    qed
    also have "... = d * p"
    proof -
      define f where "f i = sum_ab $ i $ i" for i
      have prod_rw: "prod f (insert a (UNIV - {a} - {b})) 
        =  f a * prod f (UNIV - {a} - {b})"
        by (rule prod.insert, simp_all)
      have prod1: "prod f (UNIV - {a} - {b}) = 1"
        by (rule prod.neutral) (simp add: f_def sum_ab_def row_add_def 
          mult_b_dp_def mult_row_def bezout_matrix_def Let_def)
      have ap_bq_d: "A $ a $ j * p + A $ b $ j * q = d" by (metis mult.commute pa_bq_d)
      have "prod f UNIV = prod f (insert b (insert a (UNIV - {a} - {b})))"
        using UNIV_rw by simp
      also have "... = f b * prod f (insert a (UNIV - {a} - {b}))"
      proof (rule prod.insert, simp)
        show "b  insert a (UNIV - {a} - {b})" using a_not_b by auto
      qed
      also have "... = f b * f a" unfolding prod_rw prod1 by auto
      also have "... = (d * p * v + ?b * q) * p" 
        unfolding f_def sum_ab_def row_add_def mult_b_dp_def mult_row_def bezout_matrix_def
        unfolding bz [symmetric]
        using a_not_b by auto
      also have "... = d * p" 
        using pib ap_bq_d semiring_normalization_rules(16) by auto
      finally show ?thesis unfolding f_def .
    qed
    finally have "det (sum_ab) = d * p" .
    moreover have "det (sum_ab) = d * p * det ?B"
      unfolding sum_ab_def
      unfolding det_row_add'[OF not_sym[OF a_not_b]]
      unfolding mult_b_dp_def unfolding det_mult_row ..
    ultimately show ?thesis
      by (metis (erased, opaque_lifting) False aj d_dvd_a dvd_0_left_iff mult_cancel_left1 mult_eq_0_iff)
  qed
qed

lemma invertible_bezout_matrix:
  fixes A::"'a::{bezout_ring_div}^'cols^'rows::{finite,wellorder}"
  assumes ib: "is_bezout_ext bezout"
  and a_less_b: "a < b"
  and aj: "A $ a $ j  0"
  shows "invertible (bezout_matrix A a b j bezout)"
  unfolding invertible_iff_is_unit
  unfolding det_bezout_matrix[OF assms]
  by simp

lemma echelon_form_upt_k_bezout_matrix:
  fixes A k and i::"'b::mod_type"
  assumes e: "echelon_form_upt_k A k"
  and ib: "is_bezout_ext bezout"
  and Aik_0: "A $ i $ from_nat k  0"
  and zero_i: "is_zero_row_upt_k i k A" 
  and i_less_n: "i<n"
  and k: "k<ncols A"
  shows "echelon_form_upt_k (bezout_matrix A i n (from_nat k) bezout ** A) k"
proof -
  let ?B="(bezout_matrix A i n (from_nat k) bezout ** A)"
  have i_not_n: "i  n" using i_less_n by simp
  have zero_n: "is_zero_row_upt_k n k A" 
    by (metis assms(5) e echelon_form_upt_k_condition1 zero_i)
  have zero_i2: "is_zero_row_upt_k i (to_nat (from_nat k::'c)) A" 
    using zero_i 
    by (metis k ncols_def to_nat_from_nat_id)
  have zero_n2: "is_zero_row_upt_k n (to_nat (from_nat k::'c)) A"
    using zero_n by (metis k ncols_def to_nat_from_nat_id)
  show ?thesis
    unfolding echelon_form_upt_k_def
  proof (auto)
    fix ia j
    assume ia: "is_zero_row_upt_k ia k ?B"
      and ia_j: "ia < j"
    have ia_A: "is_zero_row_upt_k ia k A"
    proof (unfold is_zero_row_upt_k_def, clarify)
      fix j::'c assume j_less_k: "to_nat j < k"
      have "A $ ia $ j = ?B $ ia $ j"
      proof (rule bezout_matrix_preserves_previous_columns
          [symmetric, OF ib i_not_n Aik_0 _ zero_i2 zero_n2])
        show "j < from_nat k" using j_less_k k 
          by (metis from_nat_mono from_nat_to_nat_id ncols_def)
      qed
      also have "... = 0" by (metis ia is_zero_row_upt_k_def j_less_k)
      finally show "A $ ia $ j = 0" .
    qed
    show "is_zero_row_upt_k j k ?B" 
    proof (unfold is_zero_row_upt_k_def, clarify)
      fix ja::'c assume ja_less_k: "to_nat ja < k" 
      have "?B $ j $ ja = A $ j $ ja"
      proof (rule bezout_matrix_preserves_previous_columns[OF ib i_not_n Aik_0 _ zero_i2 zero_n2])
        show "ja < from_nat k" using ja_less_k k 
          by (metis from_nat_mono from_nat_to_nat_id ncols_def)
      qed
      also have "... = 0" 
        by (metis e echelon_form_upt_k_condition1 ia_A ia_j is_zero_row_upt_k_def ja_less_k)
     finally show "?B $ j $ ja = 0" .
   qed
 next
   fix ia j
   assume ia_j: "ia < j"
     and not_zero_ia_B: "¬ is_zero_row_upt_k ia k ?B"
     and not_zero_j_B: "¬ is_zero_row_upt_k j k ?B"
   obtain na where na: "to_nat na < k" and Biana: "?B $ ia $ na  0" 
     using not_zero_ia_B unfolding is_zero_row_upt_k_def by auto
   obtain na2 where na2: "to_nat na2 < k" and Bjna2: "?B $ j $ na2  0" 
     using not_zero_j_B unfolding is_zero_row_upt_k_def by auto
   have na_less_fn: "na < from_nat k" by (metis from_nat_mono from_nat_to_nat_id k na ncols_def)
   have "A $ ia $ na = ?B $ ia $ na"
     by (rule bezout_matrix_preserves_previous_columns
     [symmetric, OF ib i_not_n Aik_0 na_less_fn zero_i2 zero_n2])
   also have "...  0" using Biana by simp
   finally have A: "A $ ia $ na  0" .
   have na_less_fn2: "na2 < from_nat k" by (metis from_nat_mono from_nat_to_nat_id k na2 ncols_def)
   have "A $ j $ na2 = ?B $ j $ na2"
     by (rule bezout_matrix_preserves_previous_columns
     [symmetric, OF ib i_not_n Aik_0 na_less_fn2 zero_i2 zero_n2])
   also have "...  0" using Bjna2 by simp
   finally have A2: "A $ j $ na2  0" .
   have not_zero_ia_A: "¬ is_zero_row_upt_k ia k A"
     unfolding is_zero_row_upt_k_def using na A by auto
   have not_zero_j_A: "¬ is_zero_row_upt_k j k A"
     unfolding is_zero_row_upt_k_def using na2 A2 by auto
   obtain na where A: "A $ ia $ na  0" and na_less_k: "to_nat na<k" 
     using not_zero_ia_A unfolding is_zero_row_upt_k_def by auto
   have na_less_fn: "na<from_nat k" using na_less_k
    by (metis from_nat_mono from_nat_to_nat_id k ncols_def)
   obtain na2 where A2: "A $ j $ na2  0" and na2_less_k: "to_nat na2<k" 
     using not_zero_j_A unfolding is_zero_row_upt_k_def by auto
   have na_less_fn2: "na2<from_nat k" using na2_less_k
    by (metis from_nat_mono from_nat_to_nat_id k ncols_def)
   have least_eq: "(LEAST na. ?B $ ia $ na  0) = (LEAST na. A $ ia $ na  0)"
   proof (rule Least_equality)
     have "?B $ ia $ (LEAST na. A $ ia $ na  0) = A $ ia $ (LEAST na. A $ ia $ na  0)"
     proof (rule bezout_matrix_preserves_previous_columns[OF ib i_not_n Aik_0 _ zero_i2 zero_n2])
       show "(LEAST na. A $ ia $ na  0) < from_nat k" using Least_le A na_less_fn by fastforce
     qed
     also have "...  0" by (metis (mono_tags) A LeastI)
     finally show "?B $ ia $ (LEAST na. A $ ia $ na  0)  0" .
     fix y
     assume B_ia_y: "?B $ ia $ y  0"
     show "(LEAST na. A $ ia $ na  0)  y"
     proof (cases "y<from_nat k")
       case True
       show ?thesis
       proof (rule Least_le)
         have "A $ ia $ y = ?B $ ia $ y"
           by (rule bezout_matrix_preserves_previous_columns[symmetric, 
             OF ib i_not_n Aik_0 True zero_i2 zero_n2])
         also have "...  0" using B_ia_y .
         finally show "A $ ia $ y  0" .
       qed
     next
       case False
       show ?thesis using False
         by (metis (mono_tags) A Least_le dual_order.strict_iff_order
            le_less_trans na_less_fn not_le)
     qed
   qed
   have least_eq2: "(LEAST na. ?B $ j $ na  0) = (LEAST na. A $ j $ na  0)"
   proof (rule Least_equality)
     have "?B $ j $ (LEAST na. A $ j $ na  0) = A $ j $ (LEAST na. A $ j $ na  0)"
     proof (rule bezout_matrix_preserves_previous_columns[OF ib i_not_n Aik_0 _ zero_i2 zero_n2])
       show "(LEAST na. A $ j $ na  0) < from_nat k" using Least_le A2 na_less_fn2 by fastforce
     qed
     also have "...  0" by (metis (mono_tags) A2 LeastI)
     finally show "?B $ j $ (LEAST na. A $ j $ na  0)  0" .
     fix y
     assume B_ia_y: "?B $ j $ y  0"
     show "(LEAST na. A $ j $ na  0)  y" 
     proof (cases "y<from_nat k")
       case True
       show ?thesis
       proof (rule Least_le)
         have "A $ j $ y = ?B $ j $ y"
           by (rule bezout_matrix_preserves_previous_columns[symmetric, 
            OF ib i_not_n Aik_0 True zero_i2 zero_n2])
         also have "...  0" using B_ia_y .
         finally show "A $ j $ y  0" .
       qed
     next
       case False
       show ?thesis using False
         by (metis (mono_tags) A2 Least_le dual_order.strict_iff_order
            le_less_trans na_less_fn2 not_le)
     qed
   qed   
   show "(LEAST na. ?B $ ia $ na  0) < (LEAST na. ?B $ j $ na  0)" unfolding least_eq least_eq2
     by (rule echelon_form_upt_k_condition2[OF e ia_j not_zero_ia_A not_zero_j_A])
 qed
qed

lemma bezout_matrix_preserves_rest:
  assumes ib: "is_bezout_ext bezout"
  and a_not_n: "an"
  and i_not_n: "in"
  and a_not_i: "ai"
  and Aik_0: "A $ i $ k  0"
  and zero_ikA: "is_zero_row_upt_k i (to_nat k) A"
  shows "(bezout_matrix A i n k bezout ** A) $ a $ b = A $ a $ b"
  unfolding matrix_matrix_mult_def unfolding bezout_matrix_def Let_def 
proof (auto simp add: a_not_n i_not_n a_not_i)
  have UNIV_rw: "UNIV = insert a (UNIV - {a})" by auto
  let ?f="(λk. (if a = k then 1 else 0) * A $ k $ b)"
  have sum0: "sum ?f (UNIV - {a}) = 0" by (rule sum.neutral, auto)
  have "sum ?f UNIV = sum ?f (insert a (UNIV - {a}))" using UNIV_rw by simp
  also have "... = ?f a + sum ?f (UNIV - {a})" by (rule sum.insert, simp_all)
  also have "... = ?f a" using sum0 by auto
  also have "... = A $ a $ b" by simp
  finally show "sum ?f UNIV = A $ a $ b" .
qed

text‹Code equations to execute the bezout matrix›

definition "bezout_matrix_row A a b j bezout x
  = (let (p, q, u, v, d) = bezout (A $ a $ j) (A $ b $ j) 
      in
         vec_lambda (λy. if x = a  y = a then p else
                         if x = a  y = b then q else
                         if x = b  y = a then u else
                         if x = b  y = b then v else
                         if x = y then 1 else 0))"

lemma bezout_matrix_row_code [code abstract]:
  "vec_nth (bezout_matrix_row A a b j bezout x) = 
      (let (p, q, u, v, d) = bezout (A $ a $ j) (A $ b $ j)
        in
          (λy. if x = a  y = a then p else
               if x = a  y = b then q else
               if x = b  y = a then u else
               if x = b  y = b then v else
               if x = y then 1 else 0))" unfolding bezout_matrix_row_def 
               by (cases "bezout (A $ a $ j) (A $ b $ j)") auto

lemma [code abstract]: "vec_nth (bezout_matrix A a b j bezout) = bezout_matrix_row A a b j bezout"
  unfolding bezout_matrix_def unfolding bezout_matrix_row_def[abs_def] Let_def 
  by (cases "bezout (A $ a $ j) (A $ b $ j)") auto

subsubsection‹Properties of the bezout iterate function›

lemma bezout_iterate_not_zero:
  assumes Aik_0: "A $ i $ from_nat k  0"
  and n: "n<nrows A"
  and a: "to_nat i  n"
  and ib: "is_bezout_ext bezout"
  shows "bezout_iterate A n i (from_nat k) bezout $ i $ from_nat k  0"
  using Aik_0 n a
proof (induct n arbitrary: A)
  case 0
  have "to_nat i = 0" by (metis "0.prems"(3) le_0_eq)
  hence i0: "i=0" by (metis to_nat_eq_0)
  show ?case using "0.prems"(1) unfolding i0 by auto
next
  case (Suc n)
  show ?case
  proof (cases "Suc n = to_nat i")
    case True show ?thesis unfolding bezout_iterate.simps using True Suc.prems(1) by simp
  next
    case False
    let ?B="(bezout_matrix A i (from_nat (Suc n)) (from_nat k) bezout ** A)"
    have i_le_n: "to_nat i < Suc n" using Suc.prems(3) False by auto
    have "bezout_iterate A (Suc n) i (from_nat k) bezout $ i $ from_nat k 
      = bezout_iterate ?B n i (from_nat k) bezout $ i $ from_nat k" 
      unfolding bezout_iterate.simps using i_le_n by auto
    also have "...  0"
    proof (rule Suc.hyps, rule bezout_matrix_not_zero[OF ib])
      show "i  from_nat (Suc n)" by (metis False Suc.prems(2) nrows_def to_nat_from_nat_id)
      show "A $ i $ from_nat k  0" by (rule Suc.prems(1))
      show "n < nrows ?B" by (metis Suc.prems(2) Suc_lessD nrows_def)
      show "to_nat i  n" using i_le_n by auto
    qed
    finally show ?thesis .
  qed
qed



lemma bezout_iterate_preserves:
  fixes A k and i::"'b::mod_type"
  assumes e: "echelon_form_upt_k A k"
  and ib: "is_bezout_ext bezout"
  and Aik_0: "A $ i $ from_nat k  0"
  and n: "n<nrows A"
  and "b < from_nat k"
  and i_le_n: "to_nat i  n"
  and k: "k<ncols A"
  and zero_upt_k_i: "is_zero_row_upt_k i k A" 
  shows "bezout_iterate A n i (from_nat k) bezout $ a $ b = A $ a $ b"
  using Aik_0 n i_le_n k zero_upt_k_i  e
proof (induct n arbitrary: A)
  case 0
  show ?case unfolding bezout_iterate.simps ..
next
  case (Suc n)
  show ?case
  proof (cases "Suc n = to_nat i")
    case True show ?thesis unfolding bezout_iterate.simps using True  by simp
  next
    case False
    have i_not_fn:" i  from_nat (Suc n)"
      by (metis False Suc.prems(2) nrows_def to_nat_from_nat_id)
    let ?B="(bezout_matrix A i (from_nat (Suc n)) (from_nat k) bezout ** A)"
    have i_le_n: "to_nat i < Suc n" by (metis False Suc.prems(3) le_imp_less_or_eq)
    have "bezout_iterate A (Suc n) i (from_nat k) bezout $ a $ b 
      = bezout_iterate ?B n i (from_nat k) bezout $ a $ b"
      unfolding bezout_iterate.simps using i_le_n by auto
    also have "... = ?B $ a $ b"
    proof (rule Suc.hyps)
      show "?B $ i $ from_nat k  0" 
        by (metis False Suc.prems(1) Suc.prems(2) bezout_matrix_not_zero 
            ib nrows_def to_nat_from_nat_id)   
      show "n < nrows ?B" by (metis Suc.prems(2) Suc_lessD nrows_def)
      show "k < ncols ?B" by (metis Suc.prems(4) ncols_def)
      show "to_nat i  n" using i_le_n by auto
      show "is_zero_row_upt_k i k ?B"
      proof (unfold is_zero_row_upt_k_def, clarify)
        fix j::'c assume j_k: "to_nat j < k"
        have j_k2: "j < from_nat k" by (metis from_nat_mono from_nat_to_nat_id j_k k ncols_def)
        have "?B $ i $ j = A $ i $ j" 
        proof (rule bezout_matrix_preserves_previous_columns[OF ib i_not_fn Suc.prems(1) j_k2], 
            unfold to_nat_from_nat_id[OF Suc.prems(4)[unfolded ncols_def]])
          show "is_zero_row_upt_k i k A" by (rule Suc.prems(5))
          show "is_zero_row_upt_k (from_nat (Suc n)) k A" 
            using echelon_form_upt_k_condition1[OF Suc.prems(6) Suc.prems(5)] 
            by (metis Suc.prems(2) from_nat_mono from_nat_to_nat_id i_le_n nrows_def)
        qed
        also have "... = 0" by (metis Suc.prems(5) is_zero_row_upt_k_def j_k)
        finally show "?B $ i $ j = 0" .
      qed
      show "echelon_form_upt_k ?B k" 
      proof (rule echelon_form_upt_k_bezout_matrix)
        show "echelon_form_upt_k A k" by (metis Suc.prems(6))
        show "is_bezout_ext bezout" by (rule ib)
        show "A $ i $ from_nat k  0" by (rule Suc.prems(1))
        show "is_zero_row_upt_k i k A" by (rule Suc.prems(5))
        have "(from_nat (to_nat i)::'b)from_nat (Suc n)"
          by (rule from_nat_mono'[OF Suc.prems(3) Suc.prems(2)[unfolded nrows_def]])
        thus "i < from_nat (Suc n)" using i_not_fn by auto
        show "k < ncols A" by (rule Suc.prems(4))
      qed
    qed
    also have "... = A $ a $ b"
    proof (rule bezout_matrix_preserves_previous_columns[OF ib])
      show "i  from_nat (Suc n)" by (metis False Suc.prems(2) nrows_def to_nat_from_nat_id)
      show "A $ i $ from_nat k  0" by (rule Suc.prems(1))
      show "b < from_nat k" by (rule assms(5))
      show "is_zero_row_upt_k i (to_nat (