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

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" 
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"}

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"
text ‹The subresultant PRS computes the gcd up to a scalar multiple.›

  fixes div_exp :: "'a  'a  nat  'a"
  assumes div_exp_sound: "div_exp_sound div_exp"

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
      assume "i = 1"
      thus ?thesis unfolding border[symmetric] by (intro exI[of _ 1], auto simp: numeral_2_eq_2)
      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
  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

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

context div_exp_sound_gcd

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)
    case True
    thus ?thesis unfolding gcd_impl_main_def by simp
  case True
  with len have "G2 = 0" by auto
  thus ?thesis using True unfolding gcd_impl_main_def by simp

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

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)" .
      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
      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 .
    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)

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

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 = 

corollary gcd_via_subresultant: "gcd = gcd_impl" by simp

global_interpretation div_exp_Lazard_gcd: div_exp_sound_gcd "dichotomous_Lazard :: 'a :: {semiring_gcd_mult_normalize,factorial_ring_gcd}  _" 
    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.›