Theory Kronecker_Factorization

section ‹Kronecker Factorization›

text ‹This theory contains Kronecker's factorization algorithm to factor 
  integer or rational polynomials.›

theory Kronecker_Factorization
imports 
  Polynomial_Interpolation.Polynomial_Interpolation
  Sqrt_Babylonian.Sqrt_Babylonian_Auxiliary
  Missing_List
  Prime_Factorization
  Precomputation
  Gauss_Lemma
  Dvd_Int_Poly
begin

subsection ‹Definitions›  
context 
  fixes df :: "int ⇒ int list"
  and dp :: "int ⇒ int list"
  and bnd :: nat
begin

definition kronecker_samples :: "nat ⇒ int list" where
  "kronecker_samples n ≡ let min = - int (n div 2) in [min .. min + int n]"

lemma kronecker_samples_0: "0 ∈ set (kronecker_samples n)" unfolding kronecker_samples_def by auto

text ‹Since 0 is always a samples value, we make a case analysis: 
   we only take positive divisors of $p(0)$, and consider all divisors for other $p(j)$.›
definition kronecker_factorization_main :: "int poly ⇒ int poly option" where
  "kronecker_factorization_main p ≡ if degree p ≤ 1 then None else let 
     p = primitive_part p;
     js = kronecker_samples bnd;
     cjs = map (λ j. (poly p j, j)) js
   in (case map_of cjs 0 of 
     Some j ⇒ Some ([:- j, 1 :])
   | None ⇒ let djs = map (λ (v,j). map (Pair j) (if j = 0 then dp v else df v)) cjs in 
     map_option the (find_map_filter newton_interpolation_poly_int
     (λ go. case go of None ⇒ False | Some g ⇒ dvd_int_poly_non_0 g p ∧ degree g ≥ 1) 
       (concat_lists djs)))"

definition kronecker_factorization_rat_main :: "rat poly ⇒ rat poly option" where
  "kronecker_factorization_rat_main p ≡ map_option (map_poly of_int) 
     (kronecker_factorization_main (snd (rat_to_normalized_int_poly p)))"
end

definition kronecker_factorization :: "int poly ⇒ int poly option" where
  "kronecker_factorization p = 
     kronecker_factorization_main divisors_int divisors_int_pos (degree p div 2) p"

definition kronecker_factorization_rat :: "rat poly ⇒ rat poly option" where
  "kronecker_factorization_rat p = 
     kronecker_factorization_rat_main divisors_int divisors_int_pos (degree p div 2) p"

subsection ‹Code setup for divisors›

definition "divisors_nat_copy n ≡ if n = 0 then [] else remdups_adj (sort (map prod_list (subseqs (prime_factorization_nat n))))"

lemma divisors_nat_copy[simp]: "divisors_nat_copy = divisors_nat" 
  unfolding divisors_nat_def[abs_def] divisors_nat_copy_def[abs_def] ..

definition "memo_divisors_nat ≡ memo_nat 0 100 divisors_nat_copy"

lemma memo_divisors_nat[code_unfold]: "divisors_nat = memo_divisors_nat"  
  unfolding memo_divisors_nat_def by simp


subsection ‹Proofs›

context
begin

lemma rat_to_int_poly_of_int: assumes rp: "rat_to_int_poly (map_poly of_int p) = (c,q)"
  shows "c = 1" "q = p"
proof -
  define xs where "xs = map (snd ∘ quotient_of) (coeffs (map_poly rat_of_int p))"
  have xs: "set xs ⊆ {1}" unfolding xs_def by auto
  from assms[unfolded rat_to_int_poly_def Let_def]
  have c: "c = fst (common_denom (coeffs (map_poly rat_of_int p)))" by auto
  also have "… = list_lcm xs" 
    unfolding common_denom_def Let_def xs_def by (simp add: o_assoc)
  also have "… = 1" using xs
    by (induct xs, auto)
  finally show c: "c = 1" by auto
  from rat_to_int_poly[OF rp, unfolded c] show "q = p" by auto
qed

lemma rat_to_normalized_int_poly_of_int: assumes "rat_to_normalized_int_poly (map_poly of_int p) = (c,q)"
  shows "c ∈ ℤ" "p ≠ 0 ⟹ c = of_int (content p) ∧ q = primitive_part p"
proof -
  obtain d r where ri: "rat_to_int_poly (map_poly rat_of_int p) = (d,r)" by force
  from rat_to_int_poly_of_int[OF ri]
    assms[unfolded rat_to_normalized_int_poly_def ri split] 
  show "c ∈ ℤ" "p ≠ 0 ⟹ c = of_int (content p) ∧ q = primitive_part p" 
    by (auto split: if_splits)
qed

lemma dvd_poly_int_content_1: assumes c_x: "content x = 1"
  shows "(x dvd y) = (map_poly rat_of_int x dvd map_poly of_int y)"
proof -
  let ?r = "rat_of_int"
  let ?rp = "map_poly ?r"
  show ?thesis
  proof
    assume "x dvd y"
    then obtain z where "y = x * z" unfolding dvd_def by auto
    from arg_cong[OF this, of ?rp]
    show "?rp x dvd ?rp y" by auto
  next
    assume dvd: "?rp x dvd ?rp y"
    show "x dvd y"
    proof (cases "y = 0")
      case True
      thus ?thesis by auto
    next
      case False note y0 = this
      hence "?rp y ≠ 0" by simp
      hence rx0: "?rp x ≠ 0" using dvd by auto
      hence x0: "x ≠ 0" by simp
      from dvd obtain z where prod: "?rp y = ?rp x * z" unfolding dvd_def by auto
      obtain cx xx where x: "rat_to_normalized_int_poly (?rp x) = (cx, xx)" by force
      from rat_to_int_factor_explicit[OF prod x] obtain z where y: "y = xx * smult (content y) z" by auto
      from rat_to_normalized_int_poly[OF x] rx0 have xx: "?rp x = smult cx (?rp xx)" 
        and cxx: "content xx = 1" and cx0: "cx > 0" by auto
      obtain cn cd where quot: "quotient_of cx = (cn,cd)" by force
      from quotient_of_div[OF quot] have cx: "cx = ?r cn / ?r cd" by auto
      from quotient_of_denom_pos[OF quot] have cd0: "cd > 0" by auto
      with cx cx0 have cn0: "cn > 0" by (simp add: zero_less_divide_iff)
      from arg_cong[OF xx, of "smult (?r cd)"] have "smult (?r cd) (?rp x) = smult (?r cn) (?rp xx)"
        unfolding cx using cd0 by (auto simp: field_simps)
      from this have id: "smult cd x = smult cn xx" by (fold hom_distribs, unfold of_int_poly_hom.eq_iff)
      from arg_cong[OF this, of content, unfolded content_smult_int cxx] cn0 cd0 
      have cn: "cn = cd * content x" by auto
      from quotient_of_coprime[OF quot, unfolded cn] cd0 have "cd = 1" by auto
      with cx have cx: "cx = ?r cn" by auto
      from xx[unfolded this] have x: "x = smult cn xx" by (fold hom_distribs, simp)
      from arg_cong[OF this, of content, unfolded content_smult_int c_x cxx] cn0 have "cn = 1" by auto
      with x have xx: "xx = x" by auto
      show "x dvd y" using y[unfolded xx] unfolding dvd_def by blast
    qed
  qed
qed

lemma content_x_minus_const_int[simp]: "content [: c, 1 :] = (1 :: int)"
  unfolding content_def by auto

lemma length_upto_add_nat[simp]: "length [a .. a + int n] = Suc n"
proof (induct n arbitrary: a)
  case (0 a)
  show ?case using upto.simps[of a a] by auto
next
  case (Suc n a)
  from Suc[of "a + 1"]
  show ?case using upto.simps[of a "a + int (Suc n)"] by (auto simp: ac_simps)
qed

lemma kronecker_samples: "distinct (kronecker_samples n)" "length (kronecker_samples n) = Suc n"
  unfolding kronecker_samples_def Let_def length_upto_add_nat by auto


lemma dvd_int_poly_non_0_degree_1[simp]: "degree q ≥ 1 ⟹ dvd_int_poly_non_0 q p = (q dvd p)"
  by (intro dvd_int_poly_non_0, auto)


context fixes df dp :: "int ⇒ int list"
  and bnd :: nat
begin

lemma kronecker_factorization_main_sound: assumes some: "kronecker_factorization_main df dp bnd p = Some q"
  and bnd: "degree p ≥ 2 ⟹ bnd ≥ 1"
  shows "degree q ≥ 1" "degree q ≤ bnd" "q dvd p"
proof -
  let ?r = "rat_of_int"
  let ?rp = "map_poly ?r"
  note res = some[unfolded kronecker_factorization_main_def Let_def]
  from res have dp: "degree p ≥ 2" and "(degree p ≤ 1) = False" by (auto split: if_splits)
  note res = res[unfolded this if_False]
  note bnd = bnd[OF dp]
  define P where "P = primitive_part p"
  have degP: "degree P = degree p" unfolding P_def by simp
  define js where "js = kronecker_samples bnd"
  define filt where "filt = (case_option False (λg. dvd_int_poly_non_0 g P ∧ 1 ≤ degree g))"
  define tests where "tests = concat_lists (map (λ(v, j). map (Pair j) (if j = 0 then dp v else df v)) (map (λj. (poly P j, j)) js))"
  note res = res[folded P_def, folded js_def filt_def, folded tests_def]
  let ?zero = "map (λj. (poly P j, j)) js"
  from res have res: "(case map_of ?zero 0 of
     None ⇒ map_option the (find_map_filter newton_interpolation_poly_int filt tests) | Some j ⇒ Some [:- j, 1:]) =
     Some q" by auto
  have "degree q ≥ 1 ∧ degree q ≤ bnd ∧ q dvd P"
  proof (cases "map_of ?zero 0")
    case (Some j)
    with res have q: "q = [: - j, 1 :]" by auto
    from map_of_SomeD[OF Some] have 0: "poly P j = 0" by auto
    hence "poly (?rp P) (?r j) = 0" by simp
    hence "[: - ?r j, 1 :] dvd ?rp P" using poly_eq_0_iff_dvd by blast
    also have "[: - ?r j, 1 :] = ?rp q" unfolding q by simp
    finally have dvd: "?rp q dvd ?rp P" .
    have "q dvd P"
      by (subst dvd_poly_int_content_1, insert dvd q, auto)
    with q dp bnd show ?thesis by auto
  next
    case None
    from res[unfolded None] 
    have res: "map_option the (find_map_filter newton_interpolation_poly_int filt tests) = Some q" by auto
    then obtain qq where 
      res: "find_map_filter newton_interpolation_poly_int filt tests = Some qq" and q: "q = the qq" 
      by (auto split: option.splits)
    from find_map_filter_Some[OF res] 
    have filt: "filt qq" and tests: "qq ∈ newton_interpolation_poly_int ` set tests" by auto
    from filt[unfolded filt_def] q obtain g where dvd: "g dvd P" and dg: "1 ≤ degree g" and qq: "qq = Some g" 
      by (cases qq, auto)
    from q qq have gq: "g = q" by auto
    from tests obtain t where t: "t ∈ set tests" and l: "newton_interpolation_poly_int t = Some g" unfolding qq 
      by auto
    from t[unfolded tests_def]
    have lent: "length t = length js" and "⋀ i. i < length js ⟹ map fst t ! i = js ! i" by auto
    hence id: "map fst t = js"
      by (intro nth_equalityI, auto)
    have dist: "distinct js" and lenj: "length js = Suc bnd" unfolding js_def degP
      using kronecker_samples by auto
    from newton_interpolation_poly_int_Some[OF dist[folded id] l, unfolded lent lenj]
    have "degree g ≤ bnd" by auto
    with dvd dg show ?thesis unfolding gq by auto
  qed note main = this
  thus "degree q ≥ 1" "degree q ≤ bnd" by auto
  from content_times_primitive_part[of p] have "p = smult (content p) P" unfolding P_def by auto
  with main show "q dvd p" by (metis dvd_smult)
qed

lemma kronecker_factorization_rat_main_sound: assumes 
  some: "kronecker_factorization_rat_main df dp bnd p = Some q"
  and bnd: "degree p ≥ 2 ⟹ bnd ≥ 1"
  shows "degree q ≥ 1" "degree q ≤ bnd" "q dvd p"
proof -
  let ?r = "rat_of_int"
  let ?rp = "map_poly ?r"
  let ?p = "rat_to_normalized_int_poly p"
  obtain a P where rp: "?p = (a,P)" by force
  from rat_to_normalized_int_poly[OF this] have p: "p = smult a (?rp P)" and a: "a ≠ 0" 
    and deg: "degree P = degree p" by auto
  from some[unfolded kronecker_factorization_rat_main_def rp]
  obtain Q where some: "kronecker_factorization_main df dp bnd P = Some Q" and q: "q = ?rp Q" by auto
  from kronecker_factorization_main_sound[OF some bnd] have dQ: "1 ≤ degree Q" 
    "degree Q ≤ bnd"
    and dvd: "Q dvd P" unfolding deg by auto
  from dvd obtain R where PQR: "P = Q * R" unfolding dvd_def by auto
  from p[unfolded arg_cong[OF this, of ?rp]]
  have "p = q * smult a (?rp R)" unfolding q by (auto simp: hom_distribs)
  thus "q dvd p" unfolding dvd_def by blast
  from q dQ show "degree q ≥ 1" "degree q ≤ bnd" by auto
qed


context
  assumes df: "divisors_fun df" and dpf: "divisors_pos_fun dp"
begin

lemma kronecker_factorization_main_complete: assumes 
  none: "kronecker_factorization_main df dp bnd p = None"
  and dp: "degree p ≥ 2"
  shows "¬ (∃ q. 1 ≤ degree q ∧ degree q ≤ bnd ∧ q dvd p)"
proof -
  let ?r = "rat_of_int"
  let ?rp = "map_poly ?r"
  from dp have "(degree p ≤ 1) = False" by auto
  note res = none[unfolded kronecker_factorization_main_def Let_def this if_False]
  define P where "P = primitive_part p"
  have degP: "degree P = degree p" unfolding P_def by simp
  define js where "js = kronecker_samples bnd"
  define filt where "filt = (case_option False (λg. dvd_int_poly_non_0 g P ∧ 1 ≤ degree g))"
  define tests where "tests = concat_lists (map (λ(v, j). map (Pair j) (if j = 0 then dp v else df v)) (map (λj. (poly P j, j)) js))"
  note res = res[folded P_def, folded js_def filt_def, folded tests_def]
  let ?zero = "map (λj. (poly P j, j)) js"
  from res have res: "(case map_of ?zero 0 of
     None ⇒ map_option the (find_map_filter newton_interpolation_poly_int filt tests) | Some j ⇒ Some [:- j, 1:]) =
     None" by auto
  hence zero: "map_of ?zero 0 = None" by (auto split: option.splits)
  with res have res: "find_map_filter newton_interpolation_poly_int filt tests = None" by auto
  {
    fix qq
    assume qq: "1 ≤ degree qq" "degree qq ≤ bnd" and dvd: "qq dvd p"
    define q' where "q' = primitive_part qq"
    define q where "q = (if poly q' 0 > 0 then q' else -q')"
    from qq have q': "1 ≤ degree q'" "degree q' ≤ bnd" unfolding q'_def by auto
    hence q: "1 ≤ degree q" "degree q ≤ bnd" unfolding q_def by auto
    from dvd have "qq dvd (smult (content p) P)" 
      using content_times_primitive_part[of p] unfolding P_def by simp
    from dvd_smult_int[OF _ this] dp have "q' dvd P" unfolding q'_def
      by force
    hence dvd: "q dvd P" unfolding q_def by auto      
    then obtain r where P: "P = q * r" unfolding dvd_def by auto
    {
      fix j
      assume j: "j ∈ set js"
      from P have id: "poly P j = poly q j * poly r j" by auto
      hence dvd: "poly q j dvd poly P j" by auto
      from j have "(poly P j, j) ∈ set ?zero" by auto
      with zero have zero: "poly P j ≠ 0" unfolding map_of_eq_None_iff by force
      with id have "poly q j ≠ 0" by auto
      hence "j = 0 ⟹ poly q j > 0" unfolding q_def by auto
      from divisors_funD[OF df zero dvd] divisors_pos_funD[OF dpf zero dvd this] 
      have "poly q j ∈ set (df (poly P j))" "j = 0 ⟹ poly q j ∈ set (dp (poly P j))" .
    } note mem1 = this
    define t where "t = map (λ j. (j, poly q j)) js"
    have t: "t ∈ set tests" unfolding tests_def concat_lists_listset listset length_map map_map o_def
    proof (rule, intro conjI allI impI)
      show "length t = length js" unfolding t_def by simp
      fix i
      assume i: "i < length js"
      hence jsi: "js ! i ∈ set js" by auto
      have ti: "t ! i = (js ! i, poly q (js ! i))" unfolding t_def using i by auto
      let ?f = "(λx. set (case (poly P x, x) of (v, j) ⇒ map (Pair j) (if j = 0 then dp v else df v)))"
      show "t ! i ∈ map ?f js ! i"
        unfolding ti nth_map[OF i] split using mem1[OF jsi] by auto
    qed
    have dist: "distinct js" and lenj: "length js = Suc bnd" unfolding js_def degP
      using kronecker_samples by auto
    have map_fst: "map fst t = js" unfolding t_def
      by (rule nth_equalityI, auto)
    with dist have dist: "distinct (map fst t)" by simp
    from lenj q degP have degq: "degree q < length t" unfolding t_def by auto
    from find_map_filter_None[OF res] t 
    have nfilt: "¬ filt (newton_interpolation_poly_int t)" by auto
    have qt: "⋀x y. (x, y) ∈ set t ⟹ poly q x = y" unfolding t_def by auto
    from interpolation_poly_int_None[OF dist _ qt degq, of Newton] have
      "newton_interpolation_poly_int t ≠ None" by auto
    then obtain g where lt: "newton_interpolation_poly_int t = Some g" by auto
    from newton_interpolation_poly_int_Some[OF dist lt] 
    have gt: "⋀ x y. (x, y) ∈ set t ⟹ poly g x = y" and degg: "degree g < length t" 
      using degq by auto
    from uniqueness_of_interpolation_point_list[OF dist qt degq gt degg]
    have g: "g = q" by auto
    from nfilt[unfolded lt g] have "¬ filt (Some q)" .
    from this[unfolded filt_def] q dvd have False by auto
  } note main = this
  thus ?thesis by auto
qed

lemma kronecker_factorization_rat_main_complete: assumes 
  none: "kronecker_factorization_rat_main df dp bnd p = None"
  and dp: "degree p ≥ 2"
  shows "¬ (∃ q. 1 ≤ degree q ∧ degree q ≤ bnd ∧ q dvd p)"
proof 
  assume "∃ q. 1 ≤ degree q ∧ degree q ≤ bnd ∧ q dvd p"
  then obtain q where q: "1 ≤ degree q" "degree q ≤ bnd" and dvd: "q dvd p" by auto
  from dvd obtain r where prod: "p = q * r" unfolding dvd_def by auto
  let ?r = "rat_of_int"
  let ?rp = "map_poly ?r"
  let ?p = "rat_to_normalized_int_poly p"
  obtain a P where rp: "?p = (a,P)" by force
  from rat_to_normalized_int_poly[OF this] have deg: "degree P = degree p" by auto
  from rat_to_int_factor_normalized_int_poly[OF prod rp]
    obtain g' where dvd: "g' dvd P" and dg: "degree g' = degree q" by (auto intro: dvdI)
  have "kronecker_factorization_main df dp bnd P = None" 
      using none[unfolded kronecker_factorization_rat_main_def rp] by auto
  from kronecker_factorization_main_complete[OF this dp[folded deg]] dg dvd q show False by auto
qed
end
end

lemma kronecker_factorization:
  "kronecker_factorization p = Some q ⟹ 
    degree q ≥ 1 ∧ degree q < degree p ∧ q dvd p"
  "kronecker_factorization p = None ⟹ degree p ≥ 1 ⟹ irreducible⇩d p" 
proof -
  note d = kronecker_factorization_def
  {
    assume "kronecker_factorization p = Some q"
    from kronecker_factorization_main_sound[OF this[unfolded d]]
    show "degree q ≥ 1 ∧ degree q < degree p ∧ q dvd p" by auto linarith
  }
  assume kf: "kronecker_factorization p = None" and deg: "degree p ≥ 1"
  show "irreducible⇩d p"
  proof (cases "degree p = 1")
    case True 
    thus ?thesis by (rule linear_irreducible⇩d)
  next
    case False
    with deg have "degree p ≥ 2" by auto
    with kronecker_factorization_main_complete[OF divisors_fun_int divisors_pos_fun_int kf[unfolded d] this]
    show ?thesis
      by (intro irreducible⇩dI2, auto)
  qed
qed

lemma kronecker_factorization_rat:
  "kronecker_factorization_rat p = Some q ⟹ 
    degree q ≥ 1 ∧ degree q < degree p ∧ q dvd p"
  "kronecker_factorization_rat p = None ⟹ degree p ≥ 1 ⟹ irreducible⇩d p" 
proof -
  note d = kronecker_factorization_rat_def
  {
    assume "kronecker_factorization_rat p = Some q"
    from kronecker_factorization_rat_main_sound[OF this[unfolded d]]
    show "degree q ≥ 1 ∧ degree q < degree p ∧ q dvd p" by auto linarith
  }
  assume kf: "kronecker_factorization_rat p = None" and deg: "degree p ≥ 1"
  show "irreducible⇩d p"
  proof (cases "degree p = 1")
    case True 
    thus ?thesis by (rule linear_irreducible⇩d)
  next
    case False
    with deg have "degree p ≥ 2" by auto
    with kronecker_factorization_rat_main_complete[OF divisors_fun_int divisors_pos_fun_int kf[unfolded d] this]
    show ?thesis
      by (intro irreducible⇩dI2, auto)
  qed
qed

end  
end