# 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
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
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
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;
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
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

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
```