Theory Subresultant_Gcd

section ‹Computing the Gcd via the subresultant PRS›

text ‹This theory now formalizes how the subresultant PRS can be used to calculate the gcd
  of two polynomials. Moreover, it proves the connection between resultants and gcd, namely that
  the resultant is 0 iff the degree of the gcd is non-zero.›

theory Subresultant_Gcd
imports
  Subresultant
  Polynomial_Factorization.Missing_Polynomial_Factorial
begin

subsection ‹Algorithm›

locale div_exp_sound_gcd = div_exp_sound div_exp for 
  div_exp :: "'a :: {semiring_gcd_mult_normalize,factorial_ring_gcd} ⇒ 'a ⇒ nat ⇒ 'a" 
begin
definition gcd_impl_primitive where
  [code del]: "gcd_impl_primitive G1 G2 = normalize (primitive_part (fst (subresultant_prs G1 G2)))" 

definition gcd_impl_main where
  [code del]: "gcd_impl_main G1 G2 = (if G1 = 0 then 0 else if G2 = 0 then normalize G1 else
   smult (gcd (content G1) (content G2))
     (gcd_impl_primitive (primitive_part G1) (primitive_part G2)))"

definition gcd_impl where
  "gcd_impl f g = (if length (coeffs f) ≥ length (coeffs g) then gcd_impl_main f g  else gcd_impl_main g f)"

subsection ‹Soundness Proof for @{term "gcd_impl = gcd"}›
end

locale subresultant_prs_gcd = subresultant_prs_locale2 F n δ f k β G1 G2 for
       F :: "nat ⇒ 'a ::  {factorial_ring_gcd,semiring_gcd_mult_normalize} fract poly"
    and n :: "nat ⇒ nat"
    and δ :: "nat ⇒ nat"
    and f :: "nat ⇒ 'a fract"
    and k :: nat
    and β :: "nat ⇒ 'a fract"
    and G1 G2 :: "'a poly"
begin
text ‹The subresultant PRS computes the gcd up to a scalar multiple.›

context
  fixes div_exp :: "'a ⇒ 'a ⇒ nat ⇒ 'a"
  assumes div_exp_sound: "div_exp_sound div_exp"
begin

interpretation div_exp_sound_gcd div_exp 
  using div_exp_sound by (rule div_exp_sound_gcd.intro)


lemma subresultant_prs_gcd: assumes "subresultant_prs G1 G2 = (Gk, hk)"
  shows "∃ a b. a ≠ 0 ∧ b ≠ 0 ∧ smult a (gcd G1 G2) = smult b (normalize Gk)"
proof -
  from subresultant_prs[OF div_exp_sound assms]
  have Fk: "F k = ffp Gk" and "∀ i. ∃ H. i ≠ 0 ⟶ F i = ffp H"
    and "∀ i. ∃ b. 3 ≤ i ⟶ i ≤ Suc k ⟶ β i = ff b" by auto
  from choice[OF this(2)] choice[OF this(3)] obtain H beta where
    FH: "⋀ i. i ≠ 0 ⟹ F i = ffp (H i)" and
    beta: "⋀ i. 3 ≤ i ⟹ i ≤ Suc k ⟹ β i = ff (beta i)" by auto
  from Fk FH[OF k0] FH[of 1] FH[of 2] FH[of "Suc k"] F0[of "Suc k"] F1 F2
  have border: "H k = Gk" "H 1 = G1" "H 2 = G2" "H (Suc k) = 0" by auto
  have "i ≠ 0 ⟹ i ≤ k ⟹ ∃ a b. a ≠ 0 ∧ b ≠ 0 ∧ smult a (gcd G1 G2) = smult b (gcd (H i) (H (Suc i)))" for i
  proof (induct i rule: less_induct)
    case (less i)
    from less(3) have ik: "i ≤ k" .
    from less(2) have "i = 1 ∨ i ≥ 2" by auto
    thus ?case
    proof
      assume "i = 1"
      thus ?thesis unfolding border[symmetric] by (intro exI[of _ 1], auto simp: numeral_2_eq_2)
    next
      assume i2: "i ≥ 2"
      with ik have "i - 1 < i" "i - 1 ≠ 0" and imk: "i - 1 ≤ k" by auto
      from less(1)[OF this] i2
      obtain a b where a: "a ≠ 0" and b: "b ≠ 0" and IH: "smult a (gcd G1 G2) = smult b (gcd (H (i - 1)) (H i))" by auto
      define M where "M = pseudo_mod (H (i - 1)) (H i)"
      define c where "c = β (Suc i)"
      have M: "pseudo_mod (F (i - 1)) (F i) = ffp M" unfolding to_fract_hom.pseudo_mod_hom[symmetric] M_def
         using i2 FH by auto
      have c: "c ≠ 0" using β0 unfolding c_def .
      from i2 ik have 3: "Suc i ≥ 3" "Suc i ≤ Suc k" by auto
      from pmod[OF 3]
      have pm: "smult c (F (Suc i)) = pseudo_mod (F (i - 1)) (F i)" unfolding c_def by simp
      from beta[OF 3, folded c_def] obtain d where cd: "c = ff d" by auto
      with c have d: "d ≠ 0" by auto
      from pm[unfolded cd M] FH[of "Suc i"]
      have "ffp (smult d (H (Suc i))) = ffp M" by auto
      hence pm: "smult d (H (Suc i)) = M" by (rule map_poly_hom.injectivity)
      from ik F0[of i] i2 FH[of i] have Hi0: "H i ≠ 0" by auto
      from pseudo_mod[OF this, of "H (i - 1)", folded M_def]
      obtain c Q where c: "c ≠ 0" and "smult c (H (i - 1)) = H i * Q + M" by auto
      from this[folded pm] have "smult c (H (i - 1)) = Q * H i + smult d (H (Suc i))" by simp
      from gcd_add_mult[of "H i" Q "smult d (H (Suc i))", folded this]
      have "gcd (H i) (smult c (H (i - 1))) = gcd (H i) (smult d (H (Suc i)))" .
      with gcd_smult_ex[OF c, of "H (i - 1)" "H i"] obtain e where
        e: "e ≠ 0" and "gcd (H i) (smult d (H (Suc i))) = smult e (gcd (H i) (H (i - 1)))"
        unfolding gcd.commute[of "H i"] by auto
      with gcd_smult_ex[OF d, of "H (Suc i)" "H i"] obtain c where
        c: "c ≠ 0" and "smult c (gcd (H i) (H (Suc i))) = smult e (gcd (H (i - 1)) (H i))"
        unfolding gcd.commute[of "H i"] by auto
      from arg_cong[OF this(2), of "smult b"] arg_cong[OF IH, of "smult e"]
      have "smult (e * a) (gcd G1 G2) = smult (b * c) (gcd (H i) (H (Suc i)))" unfolding smult_smult
        by (simp add: ac_simps)
      moreover have "e * a ≠ 0" "b * c ≠ 0" using a b c e by auto
      ultimately show ?thesis by blast
    qed
  qed
  from this[OF k0 le_refl, unfolded border]
  obtain a b where "a ≠ 0" "b ≠ 0" and "smult a (gcd G1 G2) = smult b (normalize Gk)" by auto
  thus ?thesis by auto
qed


lemma gcd_impl_primitive: assumes "primitive_part G1 = G1" and "primitive_part G2 = G2"
shows "gcd_impl_primitive G1 G2 = gcd G1 G2"
proof -
  let ?pp = primitive_part
  let ?c = "content"
  let ?n = normalize
  from F2 F0[of 2] k2 have G2: "G2 ≠ 0" by auto
  obtain Gk hk where sub: "subresultant_prs G1 G2 = (Gk, hk)" by force
  have impl: "gcd_impl_primitive G1 G2 = ?n (?pp Gk)" unfolding gcd_impl_primitive_def sub by auto
  from subresultant_prs_gcd[OF sub]
  obtain a b where a: "a ≠ 0" and b: "b ≠ 0" and id: "smult a (gcd G1 G2) = smult b (?n Gk)"
    by auto
  define c where "c = unit_factor (gcd G1 G2)"
  define d where "d = smult (unit_factor a) c"
  from G2 have c: "is_unit c" unfolding c_def by auto
  from arg_cong[OF id, of ?pp, unfolded primitive_part_smult primitive_part_gcd assms
     primitive_part_normalize c_def[symmetric]]
  have id: "d * gcd G1 G2 = smult (unit_factor b) (?n (?pp Gk))" unfolding d_def by simp
  have d: "is_unit d" unfolding d_def using c a
    by (simp add: is_unit_smult_iff)
  from is_unitE[OF d]
  obtain e where e: "is_unit e" and de: "d * e = 1" by metis
  define a where "a = smult (unit_factor b) e"
  from arg_cong[OF id, of "λ x. e * x"]
  have "(d * e) * gcd G1 G2 = a * (?n (?pp Gk))" by (simp add: ac_simps a_def)
  hence id: "gcd G1 G2 = a * (?n (?pp Gk))" using de by simp
  have a: "is_unit a" unfolding a_def using b e
    by (simp add: is_unit_smult_iff)
  define b where "b = unit_factor (?pp Gk)"
  have "Gk ≠ 0" using subresultant_prs[OF div_exp_sound sub] F0[OF k0] by auto
  hence b: "is_unit b" unfolding b_def by auto
  from is_unitE[OF b]
  obtain c where c: "is_unit c" and bc: "b * c = 1" by metis
  obtain d where d: "is_unit d" and dac: "d = a * c" using c a by auto
  have "gcd G1 G2 = d * (b * ?n (?pp Gk))"
    unfolding id dac using bc by (simp add: ac_simps)
  also have "b * ?n (?pp Gk) = ?pp Gk" unfolding b_def by simp
  finally have "gcd G1 G2 = d * ?pp Gk" by simp
  from arg_cong[OF this, of ?n]
  have "gcd G1 G2 = ?n (d * ?pp Gk)" by simp
  also have "… = ?n (?pp Gk)" using d
    unfolding normalize_mult by (simp add: is_unit_normalize)
  finally show ?thesis unfolding impl ..
qed
end
end

context div_exp_sound_gcd
begin

lemma gcd_impl_main: assumes len: "length (coeffs G1) ≥ length (coeffs G2)"
  shows "gcd_impl_main G1 G2 = gcd G1 G2"
proof (cases "G1 = 0")
  case G1: False
  show ?thesis
  proof (cases "G2 = 0")
    case G2: False
    let ?pp = "primitive_part"
    from G2 have G2: "?pp G2 ≠ 0" and id: "(G2 = 0) = False" by auto
    from len have len: "length (coeffs (?pp G1)) ≥ length (coeffs (?pp G2))" by simp
    from enter_subresultant_prs[OF len G2] obtain F n d f k b
      where "subresultant_prs_locale2 F n d f k b (?pp G1) (?pp G2)" by auto
    interpret subresultant_prs_locale2 F n d f k b "?pp G1" "?pp G2" by fact
    interpret subresultant_prs_gcd F n d f k b "?pp G1" "?pp G2" ..
    show ?thesis unfolding gcd_impl_main_def gcd_poly_decompose[of G1] id if_False using G1
      by (subst gcd_impl_primitive, auto intro: div_exp_sound_axioms)
  next
    case True
    thus ?thesis unfolding gcd_impl_main_def by simp
  qed
next
  case True
  with len have "G2 = 0" by auto
  thus ?thesis using True unfolding gcd_impl_main_def by simp
qed


theorem gcd_impl[simp]: "gcd_impl = gcd"
proof (intro ext)
  fix f g :: "'a poly"
  show "gcd_impl f g = gcd f g"
  proof (cases "length (coeffs f) ≥ length (coeffs g)")
    case True
    thus ?thesis unfolding gcd_impl_def gcd_impl_main[OF True] by auto
  next
    case False
    hence "length (coeffs g) ≥ length (coeffs f)" by auto
    from gcd_impl_main[OF this]
    show ?thesis unfolding gcd_impl_def gcd.commute[of f g] using False by auto
  qed
qed


text ‹The implementation also reveals an important connection between resultant and gcd.›

lemma resultant_0_gcd: "resultant (f :: 'a poly) g = 0 ⟷ degree (gcd f g) ≠ 0"
proof -
  {
    fix f g :: "'a poly"
    assume len: "length (coeffs f) ≥ length (coeffs g)"
    {
      assume g: "g ≠ 0"
      with len have f: "f ≠ 0" by auto
      let ?f = "primitive_part f"
      let ?g = "primitive_part g"
      let ?c = "content"
      from len have len: "length (coeffs ?f) ≥ length (coeffs ?g)" by simp
      obtain Gk hk where sub: "subresultant_prs ?f ?g = (Gk,hk)" by force
      have cf: "?c f ≠ 0" and cg: "?c g ≠ 0" using f g by auto
      {
        from g have "?g ≠ 0" by auto
        from enter_subresultant_prs[OF len this] obtain F n d f k b
          where "subresultant_prs_locale2 F n d f k b ?f ?g" by auto
        interpret subresultant_prs_locale2 F n d f k b ?f ?g by fact
        from subresultant_prs[OF div_exp_sound_axioms sub] have "h k = ff hk" by auto
        with h0[OF le_refl] have "hk ≠ 0" by auto
      } note hk0 = this
      have "resultant f g = 0 ⟷ resultant (smult (?c f) ?f) (smult (?c g) ?g) = 0" by simp
      also have "… ⟷ resultant ?f ?g = 0" unfolding resultant_smult_left[OF cf] resultant_smult_right[OF cg]
        using cf cg by auto
      also have "… ⟷ resultant_impl_main ?f ?g = 0" 
        unfolding resultant_impl[symmetric] resultant_impl_def resultant_impl_main_def 
        using len by auto
      also have "… ⟷ (degree Gk ≠ 0)"
        unfolding resultant_impl_main_def sub split using g hk0 by auto
      also have "degree Gk = degree (gcd_impl_primitive ?f ?g)"
        unfolding gcd_impl_primitive_def sub by simp
      also have "… = degree (gcd_impl_main f g)"
        unfolding gcd_impl_main_def using f g by auto
      also have "… = degree (gcd f g)" unfolding gcd_impl[symmetric] gcd_impl_def using len by auto
      finally have "(resultant f g = 0) = (degree (gcd f g) ≠ 0)" .
    }
    moreover
    {
      assume g: "g = 0" and f: "degree f ≠ 0"
      have "(resultant f g = 0) = (degree (gcd f g) ≠ 0)"
        unfolding g using f by auto
    }
    moreover
    {
      assume g: "g = 0" and f: "degree f = 0"
      have "(resultant f g = 0) = (degree (gcd f g) ≠ 0)"
        unfolding g using f by (auto simp: resultant_def sylvester_mat_def sylvester_mat_sub_def)
    }
    ultimately have "(resultant f g = 0) = (degree (gcd f g) ≠ 0)" by blast
  } note main = this
  show ?thesis
  proof (cases "length (coeffs f) ≥ length (coeffs g)")
    case True
    from main[OF True] show ?thesis .
  next
    case False
    hence "length (coeffs g) ≥ length (coeffs f)" by auto
    from main[OF this] show ?thesis
      unfolding gcd.commute[of g f] resultant_swap[of g f] by (simp split: if_splits)
  qed
qed


subsection ‹Code Equations›

definition "gcd_impl_rec = subresultant_prs_main_impl fst"
definition "gcd_impl_start = subresultant_prs_impl fst"

lemma gcd_impl_rec_code:
  "gcd_impl_rec Gi_1 Gi ni_1 d1_1 hi_2 = (
    let pmod = pseudo_mod Gi_1 Gi
     in
     if pmod = 0 then Gi
        else let
           ni = degree Gi;
           d1 = ni_1 - ni;
           gi_1 = lead_coeff Gi_1;
           hi_1 = (if d1_1 = 1 then gi_1 else div_exp gi_1 hi_2 d1_1);
           divisor = if d1 = 1 then gi_1 * hi_1 else if even d1 then - gi_1 * hi_1 ^ d1 else gi_1 * hi_1 ^ d1;
           Gi_p1 = sdiv_poly pmod divisor
       in gcd_impl_rec Gi Gi_p1 ni d1 hi_1)"
  unfolding gcd_impl_rec_def subresultant_prs_main_impl.simps[of _ Gi_1] split Let_def
  unfolding gcd_impl_rec_def[symmetric]
  by (rule if_cong, auto)

lemma gcd_impl_start_code:
  "gcd_impl_start G1 G2 =
     (let pmod = pseudo_mod G1 G2
         in if pmod = 0 then G2
            else let
                 n2 = degree G2;
                 n1 = degree G1;
                 d1 = n1 - n2;
                 G3 = if even d1 then - pmod else pmod;
                 pmod = pseudo_mod G2 G3
                 in if pmod = 0
                    then G3
                    else let
                           g2 = lead_coeff G2;
                           n3 = degree G3;
                           h2 = (if d1 = 1 then g2 else g2 ^ d1);
                           d2 = n2 - n3;
                           divisor = (if d2 = 1 then g2 * h2 else if even d2 then - g2 * h2 ^ d2 else g2 * h2 ^ d2);
                           G4 = sdiv_poly pmod divisor
                         in gcd_impl_rec G3 G4 n3 d2 h2)"
proof -
  obtain d1 where d1: "degree G1 - degree G2 = d1" by auto
  have id1: "(if even d1 then - pmod else pmod) = (-1)^ (d1 + 1) * (pmod :: 'a poly)" for pmod by simp
  show ?thesis
    unfolding gcd_impl_start_def subresultant_prs_impl_def gcd_impl_rec_def[symmetric] Let_def split
    unfolding d1
    unfolding id1
    by (rule if_cong, auto)
qed

lemma gcd_impl_main_code:
  "gcd_impl_main G1 G2 = (if G1 = 0 then 0 else if G2 = 0 then normalize G1 else
    let c1 = content G1;
      c2 = content G2;
      p1 = map_poly (λ x. x div c1) G1;
      p2 = map_poly (λ x. x div c2) G2
     in smult (gcd c1 c2) (normalize (primitive_part (gcd_impl_start p1 p2))))"
  unfolding gcd_impl_main_def Let_def primitive_part_def gcd_impl_start_def gcd_impl_primitive_def
    subresultant_prs_impl by simp

lemmas gcd_code_lemmas = 
  gcd_impl_main_code
  gcd_impl_start_code
  gcd_impl_rec_code
  gcd_impl_def

corollary gcd_via_subresultant: "gcd = gcd_impl" by simp
end

global_interpretation div_exp_Lazard_gcd: div_exp_sound_gcd "dichotomous_Lazard :: 'a :: {semiring_gcd_mult_normalize,factorial_ring_gcd} ⇒ _" 
  defines 
    gcd_impl_Lazard = div_exp_Lazard_gcd.gcd_impl and
    gcd_impl_main_Lazard = div_exp_Lazard_gcd.gcd_impl_main and
    gcd_impl_start_Lazard = div_exp_Lazard_gcd.gcd_impl_start and
    gcd_impl_rec_Lazard = div_exp_Lazard_gcd.gcd_impl_rec
  by (simp add: Subresultant.dichotomous_Lazard div_exp_sound_gcd_def)

declare div_exp_Lazard_gcd.gcd_code_lemmas[code]

lemmas resultant_0_gcd = div_exp_Lazard_gcd.resultant_0_gcd

thm div_exp_Lazard_gcd.gcd_via_subresultant

text ‹Note that we did not activate @{thm div_exp_Lazard_gcd.gcd_via_subresultant} as code-equation, since according to our experiments,
  the subresultant-gcd algorithm is not always more efficient than the currently active equation.
  In particular, on @{typ "int poly"} @{const gcd_impl_Lazard} performs worse, but on multi-variate polynomials,
  e.g., @{typ "int poly poly poly"}, @{const gcd_impl_Lazard} is preferable.›

end