Theory Berlekamp_Type_Based

(*
    Authors:      Jose Divasón
                  Sebastiaan Joosten
                  René Thiemann
                  Akihisa Yamada
*)
section ‹The Berlekamp Algorithm›

theory Berlekamp_Type_Based
imports
  Jordan_Normal_Form.Matrix_Kernel
  Jordan_Normal_Form.Gauss_Jordan_Elimination
  Jordan_Normal_Form.Missing_VectorSpace
  Polynomial_Factorization.Square_Free_Factorization
  Polynomial_Factorization.Missing_Multiset
  Finite_Field
  Chinese_Remainder_Poly
  Poly_Mod_Finite_Field
  "HOL-Computational_Algebra.Field_as_Ring"
begin

hide_const (open) up_ring.coeff up_ring.monom Modules.module subspace
  Modules.module_hom


subsection ‹Auxiliary lemmas›

context
  fixes g :: "'b  'a :: comm_monoid_mult"
begin
lemma prod_list_map_filter: "prod_list (map g (filter f xs)) * prod_list (map g (filter (λ x. ¬ f x) xs)) 
  = prod_list (map g xs)"
  by (induct xs, auto simp: ac_simps)

lemma prod_list_map_partition: 
  assumes "List.partition f xs = (ys, zs)"
  shows "prod_list (map g xs) = prod_list (map g ys) * prod_list (map g zs)"
  using assms  by (subst prod_list_map_filter[symmetric, of _ f], auto simp: o_def)
end

lemma coprime_id_is_unit:
  fixes a::"'b::semiring_gcd"
  shows "coprime a a  is_unit a"
  using dvd_unit_imp_unit by auto

lemma dim_vec_of_list[simp]: "dim_vec (vec_of_list x) = length x"
  by (transfer, auto)

lemma length_list_of_vec[simp]: "length (list_of_vec A) = dim_vec A"
  by (transfer', auto)

lemma list_of_vec_vec_of_list[simp]: "list_of_vec (vec_of_list a) = a"
proof -
  {
  fix aa :: "'a list"
  have "map (λn. if n < length aa then aa ! n else undef_vec (n - length aa)) [0..<length aa]
    = map ((!) aa) [0..<length aa]"
    by simp
  hence "map (λn. if n < length aa then aa ! n else undef_vec (n - length aa)) [0..<length aa] = aa"
    by (simp add: map_nth)
  }
  thus ?thesis by (transfer, simp add: mk_vec_def)
qed


context
assumes "SORT_CONSTRAINT('a::finite)"
begin

lemma inj_Poly_list_of_vec': "inj_on (Poly  list_of_vec) {v. dim_vec v = n}"
proof (rule comp_inj_on)
  show "inj_on list_of_vec {v. dim_vec v = n}"
    by (auto simp add: inj_on_def, transfer, auto simp add: mk_vec_def)
  show "inj_on Poly (list_of_vec ` {v. dim_vec v = n})"
  proof (auto simp add: inj_on_def)
    fix x y::"'c vec" assume "n = dim_vec x" and dim_xy: "dim_vec y = dim_vec x"
    and Poly_eq: "Poly (list_of_vec x) = Poly (list_of_vec y)"
    note [simp del] = nth_list_of_vec
    show "list_of_vec x = list_of_vec y"  
    proof (rule nth_equalityI, auto simp: dim_xy)
      have length_eq: "length (list_of_vec x ) = length (list_of_vec y)"
        using dim_xy by (transfer, auto)
      fix i assume "i < dim_vec x"
      thus "list_of_vec x ! i = list_of_vec y ! i" using Poly_eq unfolding poly_eq_iff coeff_Poly_eq
        using dim_xy unfolding nth_default_def by (auto, presburger)
     qed
  qed
qed

corollary inj_Poly_list_of_vec: "inj_on (Poly  list_of_vec) (carrier_vec n)"
  using inj_Poly_list_of_vec' unfolding carrier_vec_def .

lemma list_of_vec_rw_map: "list_of_vec m = map (λn. m $ n) [0..<dim_vec m]"
    by (transfer, auto simp add: mk_vec_def)

lemma degree_Poly':
assumes xs: "xs  []"
shows "degree (Poly xs) < length xs"
using xs
by (induct xs, auto intro: Poly.simps(1))

lemma vec_of_list_list_of_vec[simp]: "vec_of_list (list_of_vec a) = a"
by (transfer, auto simp add: mk_vec_def)

lemma row_mat_of_rows_list:
assumes b: "b < length A"
and nc: "i. i < length A  length (A ! i) = nc"
shows "(row (mat_of_rows_list nc A) b) = vec_of_list (A ! b)"
proof (auto simp add: vec_eq_iff)
  show "dim_col (mat_of_rows_list nc A) = length (A ! b)"
    unfolding mat_of_rows_list_def using b nc by auto
  fix i assume i: "i < length (A ! b)"
  show "row (mat_of_rows_list nc A) b $ i = vec_of_list (A ! b) $ i"
    using i b nc
    unfolding mat_of_rows_list_def row_def
    by (transfer, auto simp add: mk_vec_def mk_mat_def)
qed

lemma degree_Poly_list_of_vec:
assumes n: "x  carrier_vec n"
and n0: "n > 0"
shows "degree (Poly (list_of_vec x)) < n"
proof -
  have x_dim: "dim_vec x = n" using n by auto
  have l: "(list_of_vec x)  []"
    by (auto simp add: list_of_vec_rw_map vec_of_dim_0[symmetric] n0 n x_dim)
  have "degree (Poly (list_of_vec x)) < length (list_of_vec x)" by (rule degree_Poly'[OF l])
  also have "... = n" using x_dim by auto
  finally show ?thesis .
qed

lemma list_of_vec_nth:
  assumes i: "i < dim_vec x"
  shows "list_of_vec x ! i = x $ i"
  using i
  by (transfer, auto simp add: mk_vec_def)

lemma coeff_Poly_list_of_vec_nth':
 assumes i: "i < dim_vec x"
 shows "coeff (Poly (list_of_vec x)) i = x $ i"
 using i
 by (auto simp add: list_of_vec_nth nth_default_def)

lemma list_of_vec_row_nth:
assumes x: "x < dim_col A"
shows "list_of_vec (row A i) ! x = A $$ (i, x)"
using x unfolding row_def by (transfer', auto simp add: mk_vec_def)

lemma coeff_Poly_list_of_vec_nth:
assumes x: "x < dim_col A"
shows "coeff (Poly (list_of_vec (row A i))) x = A $$ (i, x)"
proof -
  have "coeff (Poly (list_of_vec (row A i))) x  = nth_default 0 (list_of_vec (row A i)) x"
    unfolding coeff_Poly_eq  by simp
  also have "... = A $$ (i, x)" using x list_of_vec_row_nth
    unfolding nth_default_def by (auto simp del: nth_list_of_vec)
  finally show ?thesis .
qed

lemma inj_on_list_of_vec: "inj_on list_of_vec (carrier_vec n)"
 unfolding inj_on_def unfolding list_of_vec_rw_map by auto

lemma vec_of_list_carrier[simp]: "vec_of_list x  carrier_vec (length x)"
  unfolding carrier_vec_def by simp

lemma card_carrier_vec: "card (carrier_vec n:: 'b::finite vec set) = CARD('b) ^ n"
proof -
  let ?A = "UNIV::'b set"
  let ?B = "{xs. set xs  ?A  length xs = n}"
  let ?C = "(carrier_vec n:: 'b::finite vec set)"
  have "card ?C = card ?B"
  proof -
    have "bij_betw (list_of_vec) ?C ?B"
    proof (unfold bij_betw_def, auto)
      show "inj_on list_of_vec (carrier_vec n)" by (rule inj_on_list_of_vec)
      fix x::"'b list"
      assume n: "n = length x"
      thus "x  list_of_vec ` carrier_vec (length x)"
        unfolding image_def
        by auto (rule bexI[of _ "vec_of_list x"], auto)
    qed
    thus ?thesis using bij_betw_same_card by blast
  qed
  also have "... = card ?A ^ n"
    by (rule card_lists_length_eq, simp)
  finally show ?thesis .
qed


lemma finite_carrier_vec[simp]: "finite (carrier_vec n:: 'b::finite vec set)"
  by (rule card_ge_0_finite, unfold card_carrier_vec, auto)


lemma row_echelon_form_dim0_row:
assumes "A  carrier_mat 0 n"
shows "row_echelon_form A"
using assms
unfolding row_echelon_form_def pivot_fun_def Let_def by auto

lemma row_echelon_form_dim0_col:
assumes "A  carrier_mat n 0"
shows "row_echelon_form A"
using assms
unfolding row_echelon_form_def pivot_fun_def Let_def by auto

lemma row_echelon_form_one_dim0[simp]: "row_echelon_form (1m 0)"
  unfolding row_echelon_form_def pivot_fun_def Let_def by auto

lemma Poly_list_of_vec_0[simp]: "Poly (list_of_vec (0v 0)) = [:0:]"
  by (simp add: poly_eq_iff nth_default_def)

lemma monic_normalize:
assumes "(p :: 'b :: {field,euclidean_ring_gcd} poly)  0" shows "monic (normalize p)"
by (simp add: assms normalize_poly_old_def)


lemma exists_factorization_prod_list:
fixes P::"'b::field poly list"
assumes "degree (prod_list P) > 0"
  and " u. u  set P  degree u > 0  monic u"
  and "square_free (prod_list P)"
shows "Q. prod_list Q = prod_list P  length P  length Q
            (u. u  set Q  irreducible u  monic u)"
using assms
proof (induct P)
  case Nil
  thus ?case by auto
next
  case (Cons x P)
  have sf_P: "square_free (prod_list P)"
    by (metis Cons.prems(3) dvd_triv_left prod_list.Cons mult.commute square_free_factor)
  have deg_x: "degree x > 0" using Cons.prems by auto
  have distinct_P: "distinct P"
    by (meson Cons.prems(2) Cons.prems(3) distinct.simps(2) square_free_prod_list_distinct)
  have "A. finite A  x = A  A  {q. irreducible q  monic q}"
    proof (rule monic_square_free_irreducible_factorization)
      show "monic x" by (simp add: Cons.prems(2))
      show "square_free x"
        by (metis Cons.prems(3) dvd_triv_left prod_list.Cons square_free_factor)
    qed
    from this obtain A where fin_A: "finite A"
    and xA: "x = A"
    and A: "A  {q. irreducibled q  monic q}"
    by auto
    obtain A' where s: "set A' = A" and length_A': "length A' = card A"
      using finite A distinct_card finite_distinct_list by force
  have A_not_empty: "A  {}" using xA deg_x by auto
  have x_prod_list_A': "x = prod_list A'"
  proof -
    have "x = A" using xA by simp
    also have "... = prod id A" by simp
    also have "... = prod id (set A')" unfolding s by simp
    also have "... = prod_list (map id A')"
      by (rule prod.distinct_set_conv_list, simp add: card_distinct length_A' s)
    also have "... =  prod_list A'" by auto
    finally show ?thesis .
  qed
  show ?case
  proof (cases "P = []")
    case True
    show ?thesis
    proof (rule exI[of _ "A'"], auto simp add: True)
      show "prod_list A' = x" using x_prod_list_A' by simp
      show "Suc 0  length A'" using A_not_empty using s length_A'
        by (simp add: Suc_leI card_gt_0_iff fin_A)
      show "u. u  set A'  irreducible u" using s A by auto
      show "u. u  set A'  monic u" using s A by auto
    qed
 next
  case False
  have hyp: "Q. prod_list Q = prod_list P
     length P  length Q  (u. u  set Q  irreducible u  monic u)"
  proof (rule Cons.hyps[OF _ _ sf_P])
    have set_P: "set P  {}" using False by auto
    have "prod_list P = prod_list (map id P)" by simp
    also have "... = prod id (set P)"
      using prod.distinct_set_conv_list[OF distinct_P, of id] by simp
    also have "... = (set P)" by simp
    finally have "prod_list P = (set P)" .
    hence "degree (prod_list P) = degree ((set P))" by simp
    also have "... = degree (prod id (set P))" by simp
    also have "... = (i(set P). degree (id i))"
    proof (rule degree_prod_eq_sum_degree)
      show "iset P. id i  0" using Cons.prems(2) by force
    qed
    also have "... > 0"
      by (metis Cons.prems(2) List.finite_set set_P gr0I id_apply insert_iff list.set(2) sum_pos)
    finally show "degree (prod_list P) > 0" by simp
    show "u. u  set P  degree u > 0  monic u" using Cons.prems by auto
  qed
  from this obtain Q where QP: "prod_list Q = prod_list P" and length_PQ: "length P  length Q"
  and monic_irr_Q: "(u. u  set Q  irreducible u  monic u)" by blast
  show ?thesis
  proof (rule exI[of _ "A' @ Q"], auto simp add: monic_irr_Q)
    show "prod_list A' * prod_list Q = x * prod_list P" unfolding QP x_prod_list_A' by auto
    have "length A'  0" using A_not_empty using s length_A' by auto
    thus "Suc (length P)  length A' + length Q" using QP length_PQ by linarith
    show "u. u  set A'  irreducible u" using s A by auto
    show "u. u  set A'  monic u" using s A by auto
  qed
qed
qed

lemma normalize_eq_imp_smult:
  fixes p :: "'b :: {euclidean_ring_gcd} poly"
  assumes n: "normalize p = normalize q"
  shows " c. c  0  q = smult c p"
proof(cases "p = 0")
  case True with n show ?thesis by (auto intro:exI[of _ 1])
next
  case p0: False
  have degree_eq: "degree p = degree q" using n degree_normalize by metis
  hence q0:  "q  0" using p0 n by auto
  have p_dvd_q: "p dvd q" using n by (simp add: associatedD1)
  from p_dvd_q obtain k where q: "q = k * p" unfolding dvd_def by (auto simp: ac_simps)
  with q0 have "k  0" by auto
  then have "degree k = 0"
    using degree_eq degree_mult_eq p0 q by fastforce
  then obtain c where k: "k = [: c :]" by (metis degree_0_id)
  with k  0 have "c  0" by auto
  have "q = smult c p" unfolding q k by simp
  with c  0 show ?thesis by auto
qed

lemma prod_list_normalize: 
  fixes P :: "'b :: {idom_divide,normalization_semidom_multiplicative} poly list"
  shows "normalize (prod_list P) = prod_list (map normalize P)"
proof (induct P)
  case Nil
  show ?case by auto
next
  case (Cons p P)
  have "normalize (prod_list (p # P)) = normalize p * normalize (prod_list P)"
    using normalize_mult by auto
  also have "... = normalize p * prod_list (map normalize P)" using Cons.hyps by auto
  also have "... = prod_list (normalize p # (map normalize P))" by auto
  also have "... = prod_list (map normalize (p # P))" by auto
  finally show ?case .
qed


lemma prod_list_dvd_prod_list_subset:
fixes A::"'b::comm_monoid_mult list"
assumes dA: "distinct A"
  and dB: "distinct B" (*Maybe this condition could be avoided*)
  and s: "set A  set B"
shows "prod_list A dvd prod_list B"
proof -
  have "prod_list A = prod_list (map id A)" by auto
  also have "... = prod id (set A)"
    by (rule prod.distinct_set_conv_list[symmetric, OF dA])
  also have "... dvd prod id (set B)"
    by (rule prod_dvd_prod_subset[OF _ s], auto)
  also have "... = prod_list (map id B)"
    by (rule prod.distinct_set_conv_list[OF dB])
  also have "... = prod_list B" by simp
  finally show ?thesis .
qed

end

lemma gcd_monic_constant:
  "gcd f g  {1, f}" if "monic f" and "degree g = 0"
    for f g :: "'a :: {field_gcd} poly"
proof (cases "g = 0")
  case True
  moreover from monic f have "normalize f = f"
    by (rule normalize_monic)
  ultimately show ?thesis
    by simp
next
  case False
  with degree g = 0 have "is_unit g"
    by simp
  then have "Rings.coprime f g"
    by (rule is_unit_right_imp_coprime)
  then show ?thesis
    by simp
qed

lemma distinct_find_base_vectors:
fixes A::"'a::field mat"
assumes ref: "row_echelon_form A"
  and A: "A  carrier_mat nr nc"
shows "distinct (find_base_vectors A)"
proof -
  note non_pivot_base = non_pivot_base[OF ref A]
  let ?pp = "set (pivot_positions A)"
  from A have dim: "dim_row A = nr" "dim_col A = nc" by auto
  {
    fix j j'
    assume j: "j < nc" "j  snd ` ?pp" and j': "j' < nc" "j'  snd ` ?pp" and neq: "j'  j"
    from non_pivot_base(2)[OF j] non_pivot_base(4)[OF j' j neq]
    have "non_pivot_base A (pivot_positions A) j  non_pivot_base A (pivot_positions A) j'" by auto
  }
  hence inj: "inj_on (non_pivot_base A (pivot_positions A))
     (set [j[0..<nc] . j  snd ` ?pp])" unfolding inj_on_def by auto
  thus ?thesis unfolding  find_base_vectors_def Let_def unfolding distinct_map dim by auto
qed

lemma length_find_base_vectors:
fixes A::"'a::field mat"
assumes ref: "row_echelon_form A"
  and A: "A  carrier_mat nr nc"
shows "length (find_base_vectors A) = card (set (find_base_vectors A))"
using  distinct_card[OF distinct_find_base_vectors[OF ref A]] by auto


subsection ‹Previous Results›

definition power_poly_f_mod :: "'a::field poly  'a poly  nat  'a poly" where
  "power_poly_f_mod modulus = (λa n. a ^ n mod modulus)"

lemma power_poly_f_mod_binary: "power_poly_f_mod m a n = (if n = 0 then 1 mod m
    else let (d, r) = Euclidean_Rings.divmod_nat n 2;
       rec = power_poly_f_mod m ((a * a) mod m) d in
    if r = 0 then rec else (rec * a) mod m)"
  for m a :: "'a :: {field_gcd} poly"
proof -
  note d = power_poly_f_mod_def
  show ?thesis
  proof (cases "n = 0")
    case True
    thus ?thesis unfolding d by simp
  next
    case False
    obtain q r where div: "Euclidean_Rings.divmod_nat n 2 = (q,r)" by force
    hence n: "n = 2 * q + r" and r: "r = 0  r = 1" unfolding Euclidean_Rings.divmod_nat_def by auto
    have id: "a ^ (2 * q) = (a * a) ^ q"
      by (simp add: power_mult_distrib semiring_normalization_rules)
    show ?thesis
    proof (cases "r = 0")
      case True
      show ?thesis
        using power_mod [of "a * a" m q]
        by (auto simp add: Euclidean_Rings.divmod_nat_def Let_def True n d div id)
    next
      case False
      with r have r: "r = 1" by simp
      show ?thesis
        by (auto simp add: d r div Let_def mod_simps)
          (simp add: n r mod_simps ac_simps power_mult_distrib power_mult power2_eq_square)
    qed
  qed
qed


fun power_polys where
  "power_polys mul_p u curr_p (Suc i) = curr_p #
      power_polys mul_p u ((curr_p * mul_p) mod u) i"
| "power_polys mul_p u curr_p 0 = []"

context
assumes "SORT_CONSTRAINT('a::prime_card)"
begin

lemma fermat_theorem_mod_ring [simp]:
  fixes a::"'a mod_ring"
  shows "a ^ CARD('a) = a"
proof (cases "a = 0")
  case True
  then show ?thesis by auto
next
  case False
  then show ?thesis
  proof transfer
    fix a
    assume "a  {0..<int CARD('a)}" and "a  0"
    then have a: "1  a" "a < int CARD('a)"
      by simp_all
    then have [simp]: "a mod int CARD('a) = a"
      by simp
    from a have "¬ int CARD('a) dvd a"
      by (auto simp add: zdvd_not_zless)
    then have "¬ CARD('a) dvd nat ¦a¦"
      by simp
    with a have "¬ CARD('a) dvd nat a"
      by simp
    with prime_card have "[nat a ^ (CARD('a) - 1) = 1] (mod CARD('a))"
      by (rule fermat_theorem)
    with a have "int (nat a ^ (CARD('a) - 1) mod CARD('a)) = 1"
      by (simp add: cong_def)
    with a have "a ^ (CARD('a) - 1) mod CARD('a) = 1"
      by (simp add: of_nat_mod)
    then have "a * (a ^ (CARD('a) - 1) mod int CARD('a)) = a"
      by simp
    then have "(a * (a ^ (CARD('a) - 1) mod int CARD('a))) mod int CARD('a) = a mod int CARD('a)"
      by (simp only:)
    then show "a ^ CARD('a) mod int CARD('a) = a"
      by (simp add: mod_simps semiring_normalization_rules(27))
  qed
qed


lemma mod_eq_dvd_iff_poly: "((x::'a mod_ring poly) mod n = y mod n) = (n dvd x - y)"
proof
  assume H: "x mod n = y mod n"
  hence "x mod n - y mod n = 0" by simp
  hence "(x mod n - y mod n) mod n = 0" by simp
  hence "(x - y) mod n = 0" by (simp add: mod_diff_eq)
  thus "n dvd x - y" by (simp add: dvd_eq_mod_eq_0)
next
  assume H: "n dvd x - y"
  then obtain k where k: "x-y = n*k" unfolding dvd_def by blast
  hence "x = n*k + y" using diff_eq_eq by blast
  hence "x mod n = (n*k + y) mod n" by simp
  thus "x mod n = y mod n" by (simp add: mod_add_left_eq)
qed

lemma cong_gcd_eq_poly:
  "gcd a m = gcd b m" if "[(a::'a mod_ring poly) = b] (mod m)"
  using that by (simp add: cong_def) (metis gcd_mod_left mod_by_0)

lemma coprime_h_c_poly:
fixes h::"'a mod_ring poly"
assumes "c1  c2"
shows "coprime (h - [:c1:]) (h - [:c2:])"
proof (intro coprimeI)
  fix d assume "d dvd h - [:c1:]"
  and "d dvd h - [:c2:]"
  hence "h mod d = [:c1:] mod d" and "h mod d = [:c2:] mod d"
    using mod_eq_dvd_iff_poly by simp+
  hence "[:c1:] mod d = [:c2:] mod d" by simp
  hence "d dvd [:c2 - c1:]"
    by (metis (no_types) mod_eq_dvd_iff_poly diff_pCons right_minus_eq)
  thus "is_unit d"
    by (metis (no_types) assms dvd_trans is_unit_monom_0 monom_0 right_minus_eq)
qed


lemma coprime_h_c_poly2:
fixes h::"'a mod_ring poly"
assumes "coprime (h - [:c1:]) (h - [:c2:])"
and "¬ is_unit (h - [:c1:])"
shows "c1  c2"
using assms coprime_id_is_unit by blast


lemma degree_minus_eq_right:
fixes p::"'b::ab_group_add poly"
shows "degree q < degree p  degree (p - q) = degree p"
using degree_add_eq_left[of "-q" p] degree_minus by auto

lemma coprime_prod:
  fixes A::"'a mod_ring set" and g::"'a mod_ring  'a mod_ring poly"
  assumes "xA. coprime (g a) (g x)"
  shows "coprime (g a) (prod (λx. g x) A)"
proof -
  have f: "finite A" by simp
  show ?thesis
  using f using assms
  proof (induct A)
    case (insert x A)
    have "(cinsert x A. g c) = (g x) * (cA. g c)"
      by (simp add: insert.hyps(2))
    with insert.prems show ?case
      by (auto simp: insert.hyps(3) prod_coprime_right)
  qed auto
qed


lemma coprime_prod2:
  fixes A::"'b::semiring_gcd set"
  assumes "xA. coprime (a) (x)" and f: "finite A"
  shows "coprime (a) (prod (λx. x) A)"
  using f using assms
proof (induct A)
  case (insert x A)
  have "(cinsert x A. c) = (x) * (cA. c)"
    by (simp add: insert.hyps)
  with insert.prems show ?case
    by (simp add: insert.hyps prod_coprime_right)
qed auto



lemma divides_prod:
  fixes g::"'a mod_ring  'a mod_ring poly"
  assumes "c1 c2. c1  A  c2  A  c1  c2  coprime (g c1) (g c2)"
  assumes "c A. g c dvd f"
  shows "(cA. g c) dvd f"
proof -
  have finite_A: "finite A" using finite[of A] .
  thus ?thesis using assms
  proof (induct A)
    case (insert x A)
    have "(cinsert x A. g c) =  g x * (c A. g c)"
      by (simp add: insert.hyps(2))
    also have "... dvd f"
    proof (rule divides_mult)
      show "g x dvd f" using insert.prems by auto
      show "prod g A dvd f" using insert.hyps(3) insert.prems by auto
      from insert show "Rings.coprime (g x) (prod g A)"
        by (auto intro: prod_coprime_right)
    qed
    finally show ?case .
   qed auto
qed

(*
  Proof of equation 9 in the book by Knuth
  x^p - x = (x-0)(x-1)...(x-(p-1))  (mod p)
*)

lemma poly_monom_identity_mod_p:
  "monom (1::'a mod_ring) (CARD('a)) - monom 1 1 = prod (λx. [:0,1:] - [:x:]) (UNIV::'a mod_ring set)"
  (is "?lhs = ?rhs")
proof -
  let ?f="(λx::'a mod_ring. [:0,1:] - [:x:])"
  have "?rhs dvd ?lhs"
  proof (rule divides_prod)
    {
    fix a::"'a mod_ring"
    have "poly ?lhs a = 0"
      by (simp add: poly_monom)
    hence "([:0,1:] - [:a:]) dvd ?lhs"
      using poly_eq_0_iff_dvd by fastforce
    }
    thus "xUNIV::'a mod_ring set. [:0, 1:] - [:x:] dvd monom 1 CARD('a) - monom 1 1" by fast
    show "c1 c2. c1  UNIV  c2  UNIV  c1  (c2 :: 'a mod_ring)  coprime ([:0, 1:] - [:c1:]) ([:0, 1:] - [:c2:])"
      by (auto dest!: coprime_h_c_poly[of _ _ "[:0,1:]"])
  qed
  from this obtain g where g: "?lhs = ?rhs * g" using dvdE by blast
  have degree_lhs_card: "degree ?lhs = CARD('a)"
  proof -
    have "degree (monom (1::'a mod_ring) 1) = 1" by (simp add: degree_monom_eq)
    moreover have d_c: "degree (monom (1::'a mod_ring) CARD('a)) = CARD('a)"
      by (simp add: degree_monom_eq)
    ultimately have "degree (monom (1::'a mod_ring) 1) < degree (monom (1::'a mod_ring) CARD('a))"
      using prime_card unfolding prime_nat_iff by auto
    hence "degree ?lhs = degree (monom (1::'a mod_ring) CARD('a))"
      by (rule degree_minus_eq_right)
    thus ?thesis unfolding d_c .
  qed
  have degree_rhs_card: "degree ?rhs = CARD('a)"
  proof -
    have "degree (prod ?f UNIV) = sum (degree  ?f) UNIV
       coeff (prod ?f UNIV) (sum (degree  ?f) UNIV) = 1"
      by (rule degree_prod_sum_monic, auto)
    moreover have "sum (degree  ?f) UNIV = CARD('a)" by auto
    ultimately show ?thesis by presburger
  qed
  have monic_lhs: "monic ?lhs" using degree_lhs_card by auto
  have monic_rhs: "monic ?rhs" by (rule monic_prod, simp)
  have degree_eq: "degree ?rhs = degree ?lhs" unfolding degree_lhs_card degree_rhs_card ..
  have g_not_0: "g  0" using g monic_lhs by auto
  have degree_g0: "degree g = 0"
  proof -
    have "degree (?rhs * g) = degree ?rhs + degree g"
      by (rule degree_monic_mult[OF monic_rhs g_not_0])
    thus ?thesis using degree_eq g by simp
  qed
  have monic_g: "monic g" using monic_factor g monic_lhs monic_rhs by auto
  have "g = 1" using monic_degree_0[OF monic_g] degree_g0 by simp
  thus ?thesis using g by auto
qed


(*
  Proof of equation 10 in the book by Knuth
  v(x)^p - v(x) = (v(x)-0)(v(x)-1)...(v(x)-(p-1))  (mod p)
*)


lemma poly_identity_mod_p:
  "v^(CARD('a)) - v = prod (λx. v - [:x:]) (UNIV::'a mod_ring set)"
 proof -
  have id: "monom 1 1 p v = v" "[:0, 1:] p v = v" unfolding pcompose_def
    apply (auto)
    by (simp add: fold_coeffs_def)
  have id2: "monom 1 (CARD('a)) p v = v ^ (CARD('a))" by (metis id(1) pcompose_hom.hom_power x_pow_n)
  show ?thesis using arg_cong[OF poly_monom_identity_mod_p, of "λ f. f p v"]
    unfolding pcompose_hom.hom_minus pcompose_hom.hom_prod id pcompose_const id2 .
qed



lemma coprime_gcd:
  fixes h::"'a mod_ring poly"
  assumes "Rings.coprime (h-[:c1:]) (h-[:c2:])"
  shows "Rings.coprime (gcd f(h-[:c1:])) (gcd f (h-[:c2:]))"
  using assms coprime_divisors by blast


lemma divides_prod_gcd:
  fixes h::"'a mod_ring poly"
  assumes "c1 c2. c1  A  c2  A  c1  c2 coprime (h-[:c1:]) (h-[:c2:])"
  shows "(cA. gcd f (h - [:c:])) dvd f"
proof -
  have finite_A: "finite A" using finite[of A] .
  thus ?thesis using assms
  proof (induct A)
    case (insert x A)
    have "(cinsert x A. gcd f (h - [:c:])) =  (gcd f (h - [:x:])) * (c A. gcd f (h - [:c:]))"
      by (simp add: insert.hyps(2))
    also have "... dvd f"
    proof (rule divides_mult)
      show "gcd f (h - [:x:]) dvd f" by simp
      show "(cA. gcd f (h - [:c:])) dvd f" using insert.hyps(3) insert.prems by auto
      show "Rings.coprime (gcd f (h - [:x:])) (cA. gcd f (h - [:c:]))"
        by (rule prod_coprime_right)
          (metis Berlekamp_Type_Based.coprime_h_c_poly coprime_gcd coprime_iff_coprime insert.hyps(2))
    qed
    finally show ?case .
   qed auto
qed

lemma monic_prod_gcd:
assumes f: "finite A" and f0: "(f :: 'b :: {field_gcd} poly)  0"
shows "monic (cA. gcd f (h - [:c:]))"
using f
proof (induct A)
  case (insert x A)
  have rw: "(cinsert x A. gcd f (h - [:c:]))
    = (gcd f (h - [:x:])) * (c A. gcd f (h - [:c:]))"
   by (simp add: insert.hyps)
  show ?case
  proof (unfold rw, rule monic_mult)
    show "monic (gcd f (h - [:x:]))"
      using poly_gcd_monic[of f] f0
      using insert.prems insert_iff by blast
    show "monic (cA. gcd f (h - [:c:]))"
      using insert.hyps(3) insert.prems by blast
  qed
qed auto

lemma coprime_not_unit_not_dvd:
fixes a::"'b::semiring_gcd"
assumes "a dvd b"
and "coprime b c"
and "¬ is_unit a"
shows "¬ a dvd c"
using assms coprime_divisors coprime_id_is_unit by fastforce

lemma divides_prod2:
  fixes A::"'b::semiring_gcd set"
  assumes f: "finite A"
  and "aA. a dvd c"
  and "a1 a2. a1  A  a2  A  a1  a2  coprime a1 a2"
  shows "A dvd c"
using assms
proof (induct A)
  case (insert x A)
  have "(insert x A) = x * A" by (simp add: insert.hyps(1) insert.hyps(2))
  also have "... dvd c"
  proof (rule divides_mult)
    show "x dvd c" by (simp add: insert.prems)
    show "A dvd c" using insert by auto
    from insert show "Rings.coprime x (A)"
      by (auto intro: prod_coprime_right)
  qed
  finally show ?case .
qed auto


lemma coprime_polynomial_factorization:
  fixes a1 :: "'b :: {field_gcd} poly"
  assumes  irr: "as  {q. irreducible q  monic q}"
  and "finite as" and a1: "a1  as" and a2: "a2  as" and a1_not_a2: "a1  a2"
  shows "coprime a1 a2"
proof (rule ccontr)
  assume not_coprime: "¬ coprime a1 a2"
  let ?b= "gcd a1 a2"
  have b_dvd_a1: "?b dvd a1" and b_dvd_a2: "?b dvd a2" by simp+
  have irr_a1: "irreducible a1" using a1 irr by blast
  have irr_a2: "irreducible a2" using a2 irr by blast
  have a2_not0: "a2  0" using a2 irr by auto
  have degree_a1: "degree a1  0" using irr_a1 by auto
  have degree_a2: "degree a2  0" using irr_a2 by auto
  have not_a2_dvd_a1: "¬ a2 dvd a1"
  proof (rule ccontr, simp)
    assume a2_dvd_a1: "a2 dvd a1"
    from this obtain k where k: "a1 = a2 * k" unfolding dvd_def by auto
    have k_not0: "k  0" using degree_a1 k by auto
    show False
    proof (cases "degree a2 = degree a1")
      case False
      have "degree a2 < degree a1"
        using False a2_dvd_a1 degree_a1 divides_degree
        by fastforce
      hence "¬ irreducible a1"
        using degree_a2 a2_dvd_a1 degree_a2
        by (metis degree_a1 irreducibledD(2) irreducibled_multD irreducible_connect_field k neq0_conv)
      thus False using irr_a1 by contradiction
    next
      case True
      have "degree a1 = degree a2 + degree k"
        unfolding k using degree_mult_eq[OF a2_not0 k_not0] by simp
      hence "degree k = 0" using True by simp
      hence "k = 1" using monic_factor a1 a2 irr k monic_degree_0 by auto
      hence "a1 = a2" using k by simp
      thus False using a1_not_a2 by contradiction
    qed
  qed
  have b_not0: "?b  0" by (simp add: a2_not0)
  have degree_b: "degree ?b > 0"
    using not_coprime[simplified] b_not0 is_unit_gcd is_unit_iff_degree by blast
  have "degree ?b < degree a2"
    by (meson b_dvd_a1 b_dvd_a2 irreducibleD' dvd_trans gcd_dvd_1 irr_a2 not_a2_dvd_a1 not_coprime)
  hence "¬ irreducibled a2" using degree_a2 b_dvd_a2 degree_b
    by (metis degree_smult_eq irreducibled_dvd_smult less_not_refl3)
  thus False using irr_a2 by auto
qed

(*
  Proof of equation 14 in the book by Knuth
*)
theorem Berlekamp_gcd_step:
fixes f::"'a mod_ring poly" and h::"'a mod_ring poly"
assumes hq_mod_f: "[h^(CARD('a)) = h] (mod f)" and monic_f: "monic f" and sf_f: "square_free f"
shows "f = prod (λc. gcd f (h - [:c:])) (UNIV::'a mod_ring set)"  (is "?lhs = ?rhs")
proof (cases "f=0")
  case True
  thus ?thesis using coeff_0 monic_f zero_neq_one by auto
  next
  case False note f_not_0 = False
  show ?thesis
  proof (rule poly_dvd_antisym)
    show "?rhs dvd f"
      using coprime_h_c_poly by (intro divides_prod_gcd, auto)
    have "monic ?rhs" by (rule monic_prod_gcd[OF _ f_not_0], simp)
    thus "coeff f (degree f) = coeff ?rhs (degree ?rhs)"
      using monic_f by auto
    next
    show "f dvd ?rhs"
    proof -
      let ?p = "CARD('a)"
      obtain P  where finite_P: "finite P"
      and f_desc_square_free: "f = (aP. a)"
      and P: "P  {q. irreducible q  monic q}"
        using monic_square_free_irreducible_factorization[OF monic_f sf_f] by auto
      have f_dvd_hqh: "f dvd (h^?p - h)" using hq_mod_f unfolding cong_def
        using mod_eq_dvd_iff_poly by blast
      also have hq_h_rw: "... = prod (λc. h - [:c:]) (UNIV::'a mod_ring set)"
        by (rule poly_identity_mod_p)
      finally have f_dvd_hc: "f dvd prod (λc. h - [:c:]) (UNIV::'a mod_ring set)" by simp
      have "f = P" using f_desc_square_free by simp
      also have "... dvd ?rhs"
      proof (rule divides_prod2[OF finite_P])
        show "a1 a2. a1  P  a2  P  a1  a2  coprime a1 a2"
          using coprime_polynomial_factorization[OF P finite_P] by simp
        show "aP. a dvd (cUNIV. gcd f (h - [:c:]))"
        proof
          fix fi assume fi_P: "fi  P"
          show "fi dvd ?rhs"
          proof (rule dvd_prod, auto)
            show "fi dvd f" using f_desc_square_free fi_P
             using dvd_prod_eqI finite_P by blast
            hence "fi dvd (h^?p - h)" using dvd_trans f_dvd_hqh by auto
            also have "... = prod (λc. h - [:c:]) (UNIV::'a mod_ring set)"
              unfolding hq_h_rw by simp
            finally have fi_dvd_prod_hc: "fi dvd prod (λc. h - [:c:]) (UNIV::'a mod_ring set)" .
            have irr_fi: "irreducible (fi)" using fi_P P by blast
            have fi_not_unit: "¬ is_unit fi" using irr_fi by (simp add: irreducibledD(1) poly_dvd_1)
            have fi_dvd_hc: "cUNIV::'a mod_ring set. fi dvd (h-[:c:])"
              by (rule irreducible_dvd_prod[OF _ fi_dvd_prod_hc], simp add: irr_fi)
            thus "c. fi dvd h - [:c:]" by simp
          qed
        qed
      qed
      finally show "f dvd ?rhs" .
    qed
  qed
qed


(******* Implementation of Berlekamp's algorithm (type-based version) *******)
subsection ‹Definitions›

definition berlekamp_mat :: "'a mod_ring poly  'a mod_ring mat" where
  "berlekamp_mat u = (let n = degree u;
    mul_p = power_poly_f_mod u [:0,1:] (CARD('a));
    xks = power_polys mul_p u 1 n
   in
    mat_of_rows_list n (map (λ cs. let coeffs_cs = (coeffs cs);
                                        k = n - length (coeffs cs)
                                   in (coeffs cs) @ replicate k 0) xks))"


definition berlekamp_resulting_mat :: "('a mod_ring) poly  'a mod_ring mat" where
"berlekamp_resulting_mat u = (let Q = berlekamp_mat u;
    n = dim_row Q;
    QI = mat n n (λ (i,j). if i = j then Q $$ (i,j) - 1 else Q $$ (i,j))
    in (gauss_jordan_single (transpose_mat QI)))"

definition berlekamp_basis :: "'a mod_ring poly  'a mod_ring poly list" where
  "berlekamp_basis u = (map (Poly o list_of_vec) (find_base_vectors (berlekamp_resulting_mat u)))"

lemma berlekamp_basis_code[code]: "berlekamp_basis u =
  (map (poly_of_list o list_of_vec) (find_base_vectors (berlekamp_resulting_mat u)))"
  unfolding berlekamp_basis_def poly_of_list_def ..

primrec berlekamp_factorization_main :: "nat  'a mod_ring poly list  'a mod_ring poly list  nat  'a mod_ring poly list" where
  "berlekamp_factorization_main i divs (v # vs) n = (if v = 1 then berlekamp_factorization_main i divs vs n else
    if length divs = n then divs else
    let facts = [ w . u  divs, s  [0 ..< CARD('a)], w  [gcd u (v - [:of_int s:])], w  1];
      (lin,nonlin) = List.partition (λ q. degree q = i) facts
      in lin @ berlekamp_factorization_main i nonlin vs (n - length lin))"
  | "berlekamp_factorization_main i divs [] n = divs"
  
definition berlekamp_monic_factorization :: "nat  'a mod_ring poly  'a mod_ring poly list" where
  "berlekamp_monic_factorization d f = (let
     vs = berlekamp_basis f;
     n = length vs;
     fs = berlekamp_factorization_main d [f] vs n
    in fs)"

subsection ‹Properties›

lemma power_polys_works:
fixes u::"'b::unique_euclidean_semiring"
assumes i: "i < n" and c: "curr_p = curr_p mod u" (*Equivalent to degree curr_p < degree u*)
shows "power_polys mult_p u curr_p n ! i = curr_p * mult_p ^ i mod u"
using i c
proof (induct n arbitrary: curr_p i)
  case 0 thus ?case by simp
next
  case (Suc n)
  have p_rw: "power_polys mult_p u curr_p (Suc n) ! i
      = (curr_p # power_polys mult_p u (curr_p * mult_p mod u) n) ! i"
    by simp
  show ?case
  proof (cases "i=0")
    case True
    show ?thesis using Suc.prems unfolding p_rw True by auto
  next
    case False note i_not_0 = False
    show ?thesis
    proof (cases "i < n")
      case True note i_less_n = True
      have "power_polys mult_p u curr_p (Suc n) ! i = power_polys mult_p u (curr_p * mult_p mod u) n ! (i - 1)"
        unfolding p_rw using nth_Cons_pos False by auto
      also have "... = (curr_p * mult_p mod u) * mult_p ^ (i-1) mod u"
        by (rule Suc.hyps) (auto simp add: i_less_n less_imp_diff_less)
      also have "... = curr_p * mult_p ^ i mod u"
        using False by (cases i) (simp_all add: algebra_simps mod_simps)
      finally show ?thesis .
    next
      case False
      hence i_n: "i = n" using Suc.prems by auto
      have "power_polys mult_p u curr_p (Suc n) ! i = power_polys mult_p u (curr_p * mult_p mod u) n ! (n - 1)"
          unfolding p_rw using nth_Cons_pos i_n i_not_0 by auto
      also have "... = (curr_p * mult_p mod u) * mult_p ^ (n-1) mod u"
      proof (rule Suc.hyps)
        show "n - 1 < n" using i_n i_not_0 by linarith
        show "curr_p * mult_p mod u = curr_p * mult_p mod u mod u" by simp
      qed
      also have "... = curr_p * mult_p ^ i mod u"
        using i_n [symmetric] i_not_0 by (cases i) (simp_all add: algebra_simps mod_simps)
      finally show ?thesis .
    qed
  qed
qed


lemma length_power_polys[simp]: "length (power_polys mult_p u curr_p n) = n"
  by (induct n arbitrary: curr_p, auto)


(*
  Equation 12
*)

lemma Poly_berlekamp_mat:
assumes k: "k < degree u"
shows "Poly (list_of_vec (row (berlekamp_mat u) k)) = [:0,1:]^(CARD('a) * k) mod u"
proof -
  let ?map ="(map (λcs. coeffs cs @ replicate (degree u - length (coeffs cs)) 0)
              (power_polys (power_poly_f_mod u [:0, 1:] (nat (int CARD('a)))) u 1 (degree u)))"
  have "row (berlekamp_mat u) k = row (mat_of_rows_list (degree u) ?map) k"
    by (simp add: berlekamp_mat_def Let_def)
  also have "... = vec_of_list (?map ! k)"
  proof-
    {
      fix i assume i: "i < degree u"
      then have u  0
        by auto
      let ?c= "power_polys (power_poly_f_mod u [:0, 1:] CARD('a)) u 1 (degree u) ! i"
      let ?coeffs_c="(coeffs ?c)"
      have "?c = 1*([:0, 1:] ^ CARD('a) mod u)^i mod u"
      proof (unfold power_poly_f_mod_def, rule power_polys_works[OF i])
        show "1 = 1 mod u" using k mod_poly_less by force
      qed
      also have "... = [:0, 1:] ^ (CARD('a) * i) mod u" by (simp add: power_mod power_mult)
      finally have c_rw: "?c = [:0, 1:] ^ (CARD('a) * i) mod u" .
      have "length ?coeffs_c  degree u"
      proof -
        show ?thesis
        proof (cases "?c = 0")
          case True thus ?thesis by auto
          next
          case False
          have "length ?coeffs_c = degree (?c) + 1" by (rule length_coeffs[OF False])
          also have "... = degree ([:0, 1:] ^ (CARD('a) * i) mod u) + 1" using c_rw by simp
          also have "...  degree u"
            using i < degree u u  0 degree_mod_less [of u pCons 0 1 ^ (CARD('a) * i)]
            by auto
          finally show ?thesis .
        qed
      qed
      then have "length ?coeffs_c + (degree u - length ?coeffs_c) = degree u" by auto
    }
    with k show ?thesis by (intro row_mat_of_rows_list, auto)
  qed
  finally have row_rw: "row (berlekamp_mat u) k = vec_of_list (?map ! k)" .
  have "Poly (list_of_vec (row (berlekamp_mat u) k)) = Poly (list_of_vec (vec_of_list (?map ! k)))"
    unfolding row_rw ..
  also have "... = Poly (?map ! k)" by simp
  also have "... = [:0,1:]^(CARD('a) * k) mod u"
  proof -
    let ?cs = "(power_polys (power_poly_f_mod u [:0, 1:] (nat (int CARD('a)))) u 1 (degree u)) ! k"
    let ?c = "coeffs ?cs @ replicate (degree u - length (coeffs ?cs)) 0"
    have map_k_c: "?map ! k = ?c" by (rule nth_map, simp add: k)
    have "(Poly (?map ! k)) = Poly (coeffs ?cs)" unfolding map_k_c Poly_append_replicate_0 ..
    also have "... = ?cs" by simp
    also have "... = power_polys ([:0, 1:] ^ CARD('a) mod u) u 1 (degree u) ! k"
      by (simp add: power_poly_f_mod_def)
    also have "... = 1* ([:0,1:]^(CARD('a)) mod u) ^ k mod u"
    proof (rule power_polys_works[OF k])
      show "1 = 1 mod u" using k mod_poly_less by force
    qed
    also have "... = ([:0,1:]^(CARD('a)) mod u) ^ k mod u" by auto
    also have "... = [:0,1:]^(CARD('a) * k) mod u" by (simp add: power_mod power_mult)
    finally show ?thesis .
  qed
    finally show ?thesis .
qed

corollary Poly_berlekamp_cong_mat:
assumes k: "k < degree u"
shows "[Poly (list_of_vec (row (berlekamp_mat u) k)) = [:0,1:]^(CARD('a) * k)] (mod u)"
using Poly_berlekamp_mat[OF k] unfolding cong_def by auto

lemma mat_of_rows_list_dim[simp]:
  "mat_of_rows_list n vs  carrier_mat (length vs) n"
  "dim_row (mat_of_rows_list n vs) = length vs"
  "dim_col (mat_of_rows_list n vs) = n"
  unfolding mat_of_rows_list_def by auto

lemma berlekamp_mat_closed[simp]:
  "berlekamp_mat u  carrier_mat (degree u) (degree u)"
  "dim_row (berlekamp_mat u) = degree u"
  "dim_col (berlekamp_mat u) = degree u"
 unfolding carrier_mat_def berlekamp_mat_def Let_def by auto


lemma vec_of_list_coeffs_nth:
assumes i: "i  {..degree h}" and h_not0: "h  0"
shows "vec_of_list (coeffs h) $ i = coeff h i"
proof -
  have "vec_of_list (map (coeff h) [0..<degree h] @ [coeff h (degree h)]) $ i = coeff h i"
      using i
      by (transfer', auto simp add: mk_vec_def)
         (metis (no_types, lifting) Cons_eq_append_conv coeffs_def coeffs_nth degree_0
         diff_zero length_upt less_eq_nat.simps(1) list.simps(8) list.simps(9) map_append
         nth_Cons_0 upt_Suc upt_eq_Nil_conv)
  thus "vec_of_list (coeffs h) $ i = coeff h i"
    using i h_not0
    unfolding coeffs_def by simp
qed



lemma poly_mod_sum:
  fixes x y z :: "'b::field poly"
  assumes f: "finite A"
  shows "sum f A mod z = sum (λi. f i mod z) A"
using f
by (induct, auto simp add: poly_mod_add_left)


lemma prime_not_dvd_fact:
assumes kn: "k < n" and prime_n: "prime n"
shows "¬ n dvd fact k"
using kn
proof (induct k)
  case 0
  thus ?case using prime_n unfolding prime_nat_iff by auto
next
  case (Suc k)
  show ?case
  proof (rule ccontr, unfold not_not)
    assume "n dvd fact (Suc k)"
    also have "... = Suc k * {1..k}" unfolding fact_Suc unfolding fact_prod by simp
    finally have "n dvd Suc k * {1..k}" .
    hence "n dvd Suc k  n dvd {1..k}" using prime_dvd_mult_eq_nat[OF prime_n] by blast
    moreover have  "¬ n dvd Suc k" by (simp add: Suc.prems(1) nat_dvd_not_less)
    moreover hence "¬ n dvd {1..k}" using Suc.hyps Suc.prems
      using Suc_lessD fact_prod[of k] by (metis of_nat_id)
    ultimately show False by simp
  qed
qed


lemma dvd_choose_prime:
assumes kn: "k < n" and k: "k  0" and n: "n  0" and prime_n: "prime n"
shows "n dvd (n choose k)"
proof -
  have "n dvd (fact n)" by (simp add: fact_num_eq_if n)
  moreover have "¬ n dvd (fact k * fact (n-k))"
  proof (rule ccontr, simp)
    assume "n dvd fact k * fact (n - k)"
    hence "n dvd fact k  n dvd fact (n - k)" using prime_dvd_mult_eq_nat[OF prime_n] by simp
    moreover have "¬ n dvd (fact k)" by (rule prime_not_dvd_fact[OF kn prime_n])
    moreover have "¬ n dvd fact (n - k)" using  prime_not_dvd_fact[OF _ prime_n] kn k by simp
    ultimately show False by simp
  qed
  moreover have "(fact n::nat) = fact k * fact (n-k) * (n choose k)"
    using binomial_fact_lemma kn by auto
  ultimately show ?thesis using prime_n
    by (auto simp add: prime_dvd_mult_iff)
qed



lemma add_power_poly_mod_ring:
fixes x :: "'a mod_ring poly"
shows "(x + y) ^ CARD('a) = x ^ CARD('a) + y ^ CARD('a)"
proof -
  let ?A="{0..CARD('a)}"
  let ?f="λk. of_nat (CARD('a) choose k) * x ^ k * y ^ (CARD('a) - k)"
  have A_rw: "?A = insert CARD('a) (insert 0 (?A - {0} - {CARD('a)}))"
    by fastforce
  have sum0: "sum ?f (?A - {0} - {CARD('a)}) = 0"
  proof (rule sum.neutral, rule)
    fix xa assume xa: "xa  {0..CARD('a)} - {0} - {CARD('a)}"
    have card_dvd_choose: "CARD('a) dvd  (CARD('a) choose xa)"
    proof (rule dvd_choose_prime)
      show "xa < CARD('a)" using xa by simp
      show "xa  0" using xa by simp
      show "CARD('a)  0" by simp
      show "prime CARD('a)" by (rule prime_card)
    qed
    hence rw0: "of_int (CARD('a) choose xa) = (0 :: 'a mod_ring)"
      by transfer simp
    have "of_nat (CARD('a) choose xa) = [:of_int (CARD('a) choose xa) :: 'a mod_ring:]"
      by (simp add: of_nat_poly)
    also have "... = [:0:]" using rw0 by simp
    finally show "of_nat (CARD('a) choose xa) * x ^ xa * y ^ (CARD('a) - xa) = 0" by auto
  qed
  have "(x + y)^CARD('a)
    = (k = 0..CARD('a). of_nat (CARD('a) choose k) * x ^ k * y ^ (CARD('a) - k))"
    unfolding binomial_ring by (rule sum.cong, auto)
  also have "... = sum ?f (insert CARD('a) (insert 0 (?A - {0} - {CARD('a)})))"
    using A_rw by simp
  also have "... = ?f 0 + ?f CARD('a) + sum ?f (?A - {0} - {CARD('a)})" by auto
  also have "... = x^CARD('a) + y^CARD('a)" unfolding sum0 by auto
  finally show ?thesis .
qed


lemma power_poly_sum_mod_ring:
fixes f :: "'b  'a mod_ring poly"
assumes f: "finite A"
shows "(sum f A) ^ CARD('a) = sum (λi. (f i) ^ CARD('a)) A"
using f by (induct, auto simp add: add_power_poly_mod_ring)


lemma poly_power_card_as_sum_of_monoms:
  fixes h :: "'a mod_ring poly"
  shows "h ^ CARD('a) = (idegree h. monom (coeff h i) (CARD('a)*i))"
proof -
  have "h ^ CARD('a) = (idegree h. monom (coeff h i) i) ^ CARD('a)"
    by (simp add: poly_as_sum_of_monoms)
  also have "... = (idegree h. (monom (coeff h i) i) ^ CARD('a))"
    by (simp add: power_poly_sum_mod_ring)
  also have "... = (idegree h. monom (coeff h i) (CARD('a)*i))"
  proof (rule sum.cong, rule)
    fix x assume x: "x  {..degree h}"
    show "monom (coeff h x) x ^ CARD('a) = monom (coeff h x) (CARD('a) * x)"
      by (unfold poly_eq_iff, auto simp add: monom_power)
  qed
  finally show ?thesis .
qed


lemma degree_Poly_berlekamp_le:
assumes i: "i < degree u"
shows "degree (Poly (list_of_vec (row (berlekamp_mat u) i))) < degree u"
by (metis Poly_berlekamp_mat degree_0 degree_mod_less gr_implies_not0 i linorder_neqE_nat)


(*
  Equation 12: alternative statement.
*)

lemma monom_card_pow_mod_sum_berlekamp:
assumes i: "i < degree u"
shows "monom 1 (CARD('a) * i) mod u = (j<degree u. monom ((berlekamp_mat u) $$ (i,j)) j)"
proof -
  let ?p = "Poly (list_of_vec (row (berlekamp_mat u) i))"
  have degree_not_0: "degree u  0" using i by simp
  hence set_rw: "{..degree u - 1} = {..<degree u}" by auto
  have degree_le: "degree ?p < degree u"
    by (rule degree_Poly_berlekamp_le[OF i])
  hence degree_le2: "degree ?p  degree u - 1" by auto
  have "monom 1 (CARD('a) * i) mod u = [:0, 1:] ^ (CARD('a) * i) mod u"
    using x_as_monom x_pow_n by metis
  also have "... = ?p"
    unfolding Poly_berlekamp_mat[OF i] by simp
  also have "... = (idegree u - 1. monom (coeff ?p i) i)"
        using degree_le2 poly_as_sum_of_monoms' by fastforce
  also have "... = (i<degree u. monom (coeff ?p i) i)" using set_rw by auto
  also have "... = (j<degree u. monom ((berlekamp_mat u) $$ (i,j)) j)"
  proof (rule sum.cong, rule)
    fix x assume x: "x  {..<degree u}"
    have "coeff ?p x = berlekamp_mat u $$ (i, x)"
    proof (rule coeff_Poly_list_of_vec_nth)
      show "x < dim_col (berlekamp_mat u)" using x by auto
    qed
    thus "monom (coeff ?p x) x = monom (berlekamp_mat u $$ (i, x)) x"
      by (simp add: poly_eq_iff)
  qed
  finally show ?thesis .
qed



lemma col_scalar_prod_as_sum:
assumes "dim_vec v = dim_row A"
shows "col A j  v = (i = 0..<dim_vec v. A $$ (i,j) * v $ i)"
  using assms
  unfolding col_def scalar_prod_def
  by transfer' (rule sum.cong, transfer', auto simp add: mk_vec_def mk_mat_def )

lemma row_transpose_scalar_prod_as_sum:
assumes j: "j < dim_col A" and dim_v: "dim_vec v = dim_row A"
shows "row (transpose_mat A) j  v = (i = 0..<dim_vec v. A $$ (i,j) * v $ i)"
proof -
  have "row (transpose_mat A) j  v = col A j  v" using j row_transpose by auto
  also have "... = (i = 0..<dim_vec v. A $$ (i,j) * v $ i)"
    by (rule col_scalar_prod_as_sum[OF dim_v])
  finally show ?thesis .
qed


lemma poly_as_sum_eq_monoms:
assumes ss_eq: "(i<n. monom (f i) i) = (i<n. monom (g i) i)"
and a_less_n: "a<n"
shows "f a = g a"
proof -
  let ?f="λi. if i = a then f i else 0"
  let ?g="λi. if i = a then g i else 0"
  have sum_f_0: "sum ?f ({..<n} - {a}) = 0" by (rule sum.neutral, auto)
  have "coeff (i<n. monom (f i) i) a = coeff (i<n. monom (g i) i) a"
    using ss_eq unfolding poly_eq_iff by simp
  hence "(i<n. coeff (monom (f i) i) a) = (i<n. coeff (monom (g i) i) a)"
    by (simp add: coeff_sum)
  hence 1: "(i<n. if i = a then f i else 0) = (i<n. if i = a then g i else 0)"
    unfolding coeff_monom by auto
  have set_rw: "{..<n} = (insert a ({..<n} - {a}))" using a_less_n by auto
  have "(i<n. if i = a then f i else 0) = sum ?f (insert a ({..<n} - {a}))"
    using set_rw by auto
  also have "... = ?f a + sum ?f ({..<n} - {a})"
    by (simp add: sum.insert_remove)
  also have "... = ?f a" using sum_f_0 by simp
  finally have 2: "(i<n. if i = a then f i else 0) = ?f a" .
  have "sum ?g {..<n} = sum ?g (insert a ({..<n} - {a}))"
    using set_rw by auto
  also have "... = ?g a + sum ?g ({..<n} - {a})"
    by (simp add: sum.insert_remove)
  also have "... = ?g a" using sum_f_0 by simp
  finally have 3: "(i<n. if i = a then g i else 0) = ?g a" .
  show ?thesis using 1 2 3 by auto
qed



lemma dim_vec_of_list_h:
assumes "degree h < degree u"
shows "dim_vec (vec_of_list ((coeffs h) @ replicate (degree u - length (coeffs h)) 0)) = degree u"
proof -
  have "length (coeffs h)  degree u"
    by (metis Suc_leI assms coeffs_0_eq_Nil degree_0 length_coeffs_degree
        list.size(3) not_le_imp_less order.asym)
  thus ?thesis by simp
qed




lemma vec_of_list_coeffs_nth':
assumes i: "i  {..degree h}" and h_not0: "h  0"
assumes "degree h < degree u"
shows "vec_of_list ((coeffs h) @ replicate (degree u - length (coeffs h)) 0) $ i = coeff h i"
using assms
by (transfer', auto simp add: mk_vec_def coeffs_nth length_coeffs_degree nth_append)


lemma vec_of_list_coeffs_replicate_nth_0:
assumes i: "i  {..<degree u}"
shows "vec_of_list (coeffs 0 @ replicate (degree u - length (coeffs 0)) 0) $ i = coeff 0 i"
using assms
by (transfer', auto simp add: mk_vec_def)


lemma vec_of_list_coeffs_replicate_nth:
assumes i: "i  {..<degree u}"
assumes "degree h < degree u"
shows "vec_of_list ((coeffs h) @ replicate (degree u - length (coeffs h)) 0) $ i = coeff h i"
proof (cases "h = 0")
  case True
  thus ?thesis using vec_of_list_coeffs_replicate_nth_0 i by auto
next
  case False note h_not0 = False
  show ?thesis
  proof (cases "i {..degree h}")
    case True thus ?thesis using assms vec_of_list_coeffs_nth' h_not0 by simp
  next
    case False
    have c0: "coeff h i = 0" using False le_degree by auto
    thus ?thesis
      using assms False h_not0
      by (transfer', auto simp add: mk_vec_def length_coeffs_degree nth_append c0)
  qed
qed


(*
  Equation 13
*)

lemma equation_13:
  fixes u