Theory Subresultant
section ‹Subresultants and the subresultant PRS›
text ‹This theory contains most of the soundness proofs of the subresultant PRS algorithm,
where we closely follow the papers of Brown \<^cite>‹"Brown"› and Brown and Traub \<^cite>‹"BrownTraub"›.
This is in contrast to a similar Coq formalization of Mahboubi \<^cite>‹"Mahboubi06"› which
is based on polynomial determinants.
Whereas the current file only contains an algorithm to compute the resultant of two polynomials
efficiently, there is another theory ``Subresultant-Gcd'' which also contains the algorithm
to compute the GCD of two polynomials via the subresultant algorithm. In both algorithms we
integrate Lazard's optimization in the dichotomous version,
but not the second optmization described by Ducos \<^cite>‹"Ducos"›.›
theory Subresultant
imports
Resultant_Prelim
Dichotomous_Lazard
Binary_Exponentiation
More_Homomorphisms
Coeff_Int
begin
subsection ‹Algorithm›
locale div_exp_param =
fixes div_exp :: "'a :: idom_divide ⇒ 'a ⇒ nat ⇒ 'a"
begin
partial_function(tailrec) subresultant_prs_main where
"subresultant_prs_main f g c = (let
m = degree f;
n = degree g;
lf = lead_coeff f;
lg = lead_coeff g;
δ = m - n;
d = div_exp lg c δ;
h = pseudo_mod f g
in if h = 0 then (g,d)
else subresultant_prs_main g (sdiv_poly h ((-1) ^ (δ + 1) * lf * (c ^ δ))) d)"
definition subresultant_prs where
"subresultant_prs f g = (let
h = pseudo_mod f g;
δ = (degree f - degree g);
d = lead_coeff g ^ δ
in if h = 0 then (g,d)
else subresultant_prs_main g ((- 1) ^ (δ + 1) * h) d)"
definition resultant_impl_main where
"resultant_impl_main G1 G2 = (if G2 = 0 then (if degree G1 = 0 then 1 else 0) else
case subresultant_prs G1 G2 of
(Gk,hk) ⇒ (if degree Gk = 0 then hk else 0))"
definition resultant_impl where
"resultant_impl f g =
(if length (coeffs f) ≥ length (coeffs g) then resultant_impl_main f g
else let res = resultant_impl_main g f in
if even (degree f) ∨ even (degree g) then res else - res)"
end
locale div_exp_sound = div_exp_param +
assumes div_exp: "⋀ x y n.
(to_fract x)^n / (to_fract y)^(n-1) ∈ range to_fract
⟹ to_fract (div_exp x y n) = (to_fract x)^n / (to_fract y)^(n-1)"
definition basic_div_exp :: "'a :: idom_divide ⇒ 'a ⇒ nat ⇒ 'a" where
"basic_div_exp x y n = x^n div y^(n-1)"
text ‹We have an instance for arbitrary integral domains.›
lemma basic_div_exp: "div_exp_sound basic_div_exp"
by (unfold_locales, unfold basic_div_exp_def, rule sym, rule div_divide_to_fract, auto simp: hom_distribs)
text ‹Lazard's optimization is only proven for factorial rings.›
lemma dichotomous_Lazard: "div_exp_sound (dichotomous_Lazard :: 'a :: factorial_ring_gcd ⇒ _)"
by (unfold_locales, rule dichotomous_Lazard)
subsection ‹Soundness Proof for @{term "div_exp_param.resultant_impl div_exp = resultant"}›
abbreviation pdivmod :: "'a::field poly ⇒ 'a poly ⇒ 'a poly × 'a poly"
where
"pdivmod p q ≡ (p div q, p mod q)"
lemma even_sum_list: assumes "⋀ x. x ∈ set xs ⟹ even (f x) = even (g x)"
shows "even (sum_list (map f xs)) = even (sum_list (map g xs))"
using assms by (induct xs, auto)
lemma for_all_Suc: "P i ⟹ (∀ j ≥ Suc i. P j) = (∀ j ≥ i. P j)" for P
by (metis (full_types) Suc_le_eq less_le)
lemma pseudo_mod_left_0[simp]: "pseudo_mod 0 x = 0"
unfolding pseudo_mod_def pseudo_divmod_def
by (cases "x = 0"; cases "length (coeffs x)", auto)
lemma pseudo_mod_right_0[simp]: "pseudo_mod x 0 = x"
unfolding pseudo_mod_def pseudo_divmod_def by simp
lemma snd_pseudo_divmod_main_cong:
assumes "a1 = b1" "a3 = b3" "a4 = b4" "a5 = b5" "a6 = b6"
shows "snd (pseudo_divmod_main a1 a2 a3 a4 a5 a6) = snd (pseudo_divmod_main b1 b2 b3 b4 b5 b6)"
using assms snd_pseudo_divmod_main by metis
lemma snd_pseudo_mod_smult_invar_right:
shows "(snd (pseudo_divmod_main (x * lc) q r (smult x d) dr n))
= snd (pseudo_divmod_main lc q' (smult (x^n) r) d dr n)"
proof(induct n arbitrary: q q' r dr)
case (Suc n)
let ?q = "smult (x * lc) q + monom (coeff r dr) n"
let ?r = "smult (x * lc) r - (smult x (monom (coeff r dr) n * d))"
let ?dr = "dr - 1"
let ?rec_lhs = "pseudo_divmod_main (x * lc) ?q ?r (smult x d) ?dr n"
let ?rec_rhs = "pseudo_divmod_main lc q' (smult (x^n) ?r) d ?dr n"
have [simp]: "⋀ n. x ^ n * (x * lc) = lc * (x * x ^ n)"
"⋀ n c. x ^ n * (x * c) = x * x ^ n * c"
"⋀ n. x * x ^ n * lc = lc * (x * x ^ n)"
by (auto simp: ac_simps)
have "snd (pseudo_divmod_main (x * lc) q r (smult x d) dr (Suc n)) = snd ?rec_lhs"
by (auto simp:Let_def)
also have "… = snd ?rec_rhs" using Suc by auto
also have "… = snd (pseudo_divmod_main lc q' (smult (x^Suc n) r) d dr (Suc n))"
unfolding pseudo_divmod_main.simps Let_def
proof(rule snd_pseudo_divmod_main_cong,goal_cases)
case 2
show ?case by (auto simp:smult_add_right smult_diff_right smult_monom smult_monom_mult)
qed auto
finally show ?case by auto
qed auto
lemma snd_pseudo_mod_smult_invar_left:
shows "snd (pseudo_divmod_main lc q (smult x r) d dr n)
= smult x (snd (pseudo_divmod_main lc q' r d dr n))"
proof(induct n arbitrary:x lc q q' r d dr)
case (Suc n)
have sm:"smult lc (smult x r) - monom (coeff (smult x r) dr) n * d
=smult x (smult lc r - monom (coeff r dr) n * d) "
by (auto simp:smult_diff_right smult_monom smult_monom_mult mult.commute[of lc x])
let ?q' = "smult lc q' + monom (coeff r dr) n"
show ?case unfolding pseudo_divmod_main.simps Let_def Suc(1)[of lc _ _ _ _ _ ?q'] sm by auto
qed auto
lemma snd_pseudo_mod_smult_left[simp]:
shows "snd (pseudo_divmod (smult (x::'a::idom) p) q) = (smult x (snd (pseudo_divmod p q)))"
unfolding pseudo_divmod_def
by (auto simp:snd_pseudo_mod_smult_invar_left[of _ _ _ _ _ _ _ 0] Polynomial.coeffs_smult)
lemma pseudo_mod_smult_right:
assumes "(x::'a::idom)≠0" "q≠0"
shows "(pseudo_mod p (smult (x::'a::idom) q)) = (smult (x^(Suc (length (coeffs p)) - length (coeffs q))) (pseudo_mod p q))"
unfolding pseudo_divmod_def pseudo_mod_def
by (auto simp:snd_pseudo_mod_smult_invar_right[of _ _ _ _ _ _ _ 0]
snd_pseudo_mod_smult_invar_left[of _ _ _ _ _ _ _ 0] Polynomial.coeffs_smult assms)
lemma pseudo_mod_zero[simp]:
"pseudo_mod 0 f = (0::'a :: {idom} poly)"
"pseudo_mod f 0 = f"
unfolding pseudo_mod_def snd_pseudo_mod_smult_left[of 0 _ f,simplified]
unfolding pseudo_divmod_def by auto
lemma prod_combine:
assumes "j ≤ i"
shows "f i * (∏l←[j..<i]. (f l :: 'a::comm_monoid_mult)) = prod_list (map f [j..<Suc i])"
proof(subst prod_list_map_remove1[of i "[j..<Suc i]" f],goal_cases)
case 2
have "remove1 i ([j..<i] @ [i]) = [j..<i]" by (simp add: remove1_append)
thus ?case by auto
qed (insert assms, auto)
lemma prod_list_minus_1_exp: "prod_list (map (λ i. (-1)^(f i)) xs)
= (-1)^(sum_list (map f xs))"
by (induct xs, auto simp: power_add)
lemma minus_1_power_even: "(- (1 :: 'b :: comm_ring_1))^ k = (if even k then 1 else (-1))"
by auto
lemma minus_1_even_eqI: assumes "even k = even l" shows
"(- (1 :: 'b :: comm_ring_1))^k = (- 1)^l"
unfolding minus_1_power_even assms by auto
lemma (in comm_monoid_mult) prod_list_multf:
"(∏x←xs. f x * g x) = prod_list (map f xs) * prod_list (map g xs)"
by (induct xs) (simp_all add: algebra_simps)
lemma inverse_prod_list: "inverse (prod_list xs) = prod_list (map inverse (xs :: 'a :: field list))"
by (induct xs, auto)
definition pow_int :: "'a :: field ⇒ int ⇒ 'a" where
"pow_int x e = (if e < 0 then 1 / (x ^ (nat (-e))) else x ^ (nat e))"
lemma pow_int_0[simp]: "pow_int x 0 = 1" unfolding pow_int_def by auto
lemma pow_int_1[simp]: "pow_int x 1 = x" unfolding pow_int_def by auto
lemma exp_pow_int: "x ^ n = pow_int x n"
unfolding pow_int_def by auto
lemma pow_int_add: assumes x: "x ≠ 0" shows "pow_int x (a + b) = pow_int x a * pow_int x b"
proof -
have *:
"¬ a + b < 0 ⟹ a < 0 ⟹ nat b = nat (a + b) + nat (-a)"
"¬ a + b < 0 ⟹ b < 0 ⟹ nat a = nat (a + b) + nat (-b)"
"a + b < 0 ⟹ ¬ a < 0 ⟹ nat (-b) = nat a + nat (-a -b) "
"a + b < 0 ⟹ ¬ b < 0 ⟹ nat (-a) = nat b + nat (-a -b) "
by auto
have pow_eq: "l = m ⟹ (x ^ l = x ^ m)" for l m by auto
from x show ?thesis unfolding pow_int_def
by (auto split: if_splits simp: power_add[symmetric] simp: * intro!: pow_eq, auto simp: power_add)
qed
lemma pow_int_mult: "pow_int (x * y) a = pow_int x a * pow_int y a"
unfolding pow_int_def by (cases "a < 0", auto simp: power_mult_distrib)
lemma pow_int_base_1[simp]: "pow_int 1 a = 1"
unfolding pow_int_def by (cases "a < 0", auto)
lemma pow_int_divide: "a / pow_int x b = a * pow_int x (-b)"
unfolding pow_int_def by (cases b rule: linorder_cases[of _ 0], auto)
lemma divide_prod_assoc: "x / (y * z :: 'a :: field) = x / y / z" by (simp add: field_simps)
lemma minus_1_inverse_pow[simp]: "x / (-1)^n = (x :: 'a :: field) * (-1)^n"
by (simp add: minus_1_power_even)
definition subresultant_mat :: "nat ⇒ 'a :: comm_ring_1 poly ⇒ 'a poly ⇒ 'a poly mat" where
"subresultant_mat J F G = (let
dg = degree G; df = degree F; f = coeff_int F; g = coeff_int G; n = (df - J) + (dg - J)
in mat n n (λ (i,j). if j < dg - J then
if i = n - 1 then monom 1 (dg - J - 1 - j) * F else [: f (df - int i + int j) :]
else let jj = j - (dg - J) in
if i = n - 1 then monom 1 (df - J - 1 - jj) * G else [: g (dg - int i + int jj) :]))"
lemma subresultant_mat_dim[simp]:
fixes j p q
defines "S ≡ subresultant_mat j p q"
shows "dim_row S = (degree p - j) + (degree q - j)" and "dim_col S = (degree p - j) + (degree q - j)"
unfolding S_def subresultant_mat_def Let_def by auto
definition subresultant'_mat :: "nat ⇒ nat ⇒ 'a :: comm_ring_1 poly ⇒ 'a poly ⇒ 'a mat" where
"subresultant'_mat J l F G = (let
γ = degree G; φ = degree F; f = coeff_int F; g = coeff_int G; n = (φ - J) + (γ - J)
in mat n n (λ (i,j). if j < γ - J then
if i = n - 1 then (f (l - int (γ - J - 1) + int j)) else (f (φ - int i + int j))
else let jj = j - (γ - J) in
if i = n - 1 then (g (l - int (φ - J - 1) + int jj)) else (g (γ - int i + int jj))))"
lemma subresultant_index_mat:
fixes F G
assumes i: "i < (degree F - J) + (degree G - J)" and j: "j < (degree F - J) + (degree G - J)"
shows "subresultant_mat J F G $$ (i,j) =
(if j < degree G - J then
if i = (degree F - J) + (degree G - J) - 1 then monom 1 (degree G - J - 1 - j) * F else ([: coeff_int F ( degree F - int i + int j) :])
else let jj = j - (degree G - J) in
if i = (degree F - J) + (degree G - J) - 1 then monom 1 ( degree F - J - 1 - jj) * G else ([: coeff_int G (degree G - int i + int jj) :]))"
unfolding subresultant_mat_def Let_def
unfolding index_mat(1)[OF i j] split by auto
definition subresultant :: "nat ⇒ 'a :: comm_ring_1 poly ⇒ 'a poly ⇒ 'a poly" where
"subresultant J F G = det (subresultant_mat J F G)"
lemma subresultant_smult_left: assumes "(c :: 'a :: {comm_ring_1, semiring_no_zero_divisors}) ≠ 0"
shows "subresultant J (smult c f) g = smult (c ^ (degree g - J)) (subresultant J f g)"
proof -
let ?df = "degree f"
let ?dg = "degree g"
let ?n = "(?df - J) + (?dg - J)"
let ?m = "?dg - J"
let ?M = "mat ?n ?n (λ (i,j). if i = j then if i < ?m then [:c:] else 1 else 0)"
from ‹c ≠ 0› have deg: "degree (smult c f) = ?df" by simp
let ?S = "subresultant_mat J f g"
let ?cS = "subresultant_mat J (smult c f) g"
have dim: "dim_row ?S = ?n" "dim_col ?S = ?n" "dim_row ?cS = ?n" "dim_col ?cS = ?n" using deg by auto
hence C: "?S ∈ carrier_mat ?n ?n" "?cS ∈ carrier_mat ?n ?n" "?M ∈ carrier_mat ?n ?n" by auto
have dim': "dim_row (?S * ?M) = ?n" "dim_col (?S * ?M) = ?n" using dim (1,2) by simp_all
define S where "S = ?S"
have "?cS = ?S * ?M"
proof (rule eq_matI, unfold dim' dim)
fix i j
assume ij: "i < ?n" "j < ?n"
have "(?S * ?M) $$ (i,j) = row ?S i ∙ col ?M j"
by (rule index_mult_mat, insert ij dim, auto)
also have "… = (∑k = 0..<?n. row S i $ k * col ?M j $ k)" unfolding scalar_prod_def S_def[symmetric]
by simp
also have "… = (∑k = 0..<?n. S $$ (i,k) * ?M $$ (k,j))"
by (rule sum.cong, insert ij, auto simp: S_def)
also have "… = S $$ (i,j) * ?M $$ (j,j) + sum (λ k. S $$ (i,k) * ?M $$ (k,j)) ({0..<?n} - {j})"
by (rule sum.remove, insert ij, auto)
also have "… = S $$ (i,j) * ?M $$ (j,j)"
by (subst sum.neutral, insert ij, auto)
also have "… = ?cS $$ (i,j)" unfolding subresultant_index_mat[OF ij] S_def
by (subst subresultant_index_mat, unfold deg, insert ij, auto)
finally show "?cS $$ (i,j) = (?S * ?M) $$ (i,j)" by simp
qed auto
from arg_cong[OF this, of det] det_mult[OF C(1) C(3)]
have "subresultant J (smult c f) g = subresultant J f g * det ?M"
unfolding subresultant_def by auto
also have "det ?M = [:c ^ ?m :]"
proof (subst det_upper_triangular[OF _ C(3)])
show "upper_triangular ?M"
by (rule upper_triangularI, auto)
have "prod_list (diag_mat ?M) = (∏k = 0..<?n. (?M $$ (k,k)))"
unfolding prod_list_diag_prod by simp
also have "… = (∏k = 0..<?m. ?M $$ (k,k)) * (∏k = ?m..<?n. ?M $$ (k,k))"
by (subst prod.union_disjoint[symmetric], (auto)[3], rule prod.cong, auto)
also have "(∏k = 0..<?m. ?M $$ (k,k)) = (∏k = 0..<?m. [: c :])"
by (rule prod.cong, auto)
also have "(∏k = 0..<?m. [: c :]) = [: c :] ^ ?m" by simp
also have "(∏k = ?m..<?n. ?M $$ (k,k)) = (∏k = ?m..<?n. 1)"
by (rule prod.cong, auto)
finally show "prod_list (diag_mat ?M) = [: c^?m :]" unfolding poly_const_pow by simp
qed
finally show ?thesis by simp
qed
lemma subresultant_swap:
shows "subresultant J f g = smult ((- 1) ^ ((degree f - J) * (degree g - J))) (subresultant J g f)"
proof -
let ?A = "subresultant_mat J f g"
let ?k = "degree f - J"
let ?n = "degree g - J"
have nk: "?n + ?k = ?k + ?n" by simp
have change: "j < ?k + ?n ⟹ ((if j < ?k then j + ?n else j - ?k) < ?n)
= (¬ (j < ?k))" for j by auto
have "subresultant J f g = det ?A" unfolding subresultant_def by simp
also have "… = (-1)^(?k * ?n) * det (mat (?k + ?n) (?k + ?n) (λ (i,j).
?A $$ (i,(if j < ?k then j + ?n else j - ?k))))" (is "_ = _ * det ?B")
by (rule det_swap_cols, auto simp: subresultant_mat_def Let_def)
also have "?B = subresultant_mat J g f"
unfolding subresultant_mat_def Let_def
by (rule eq_matI, unfold dim_row_mat dim_col_mat nk index_mat split,
subst index_mat, (auto)[2], unfold split, subst change, force,
unfold if_conn, rule if_cong[OF refl if_cong if_cong], auto)
also have "det … = subresultant J g f" unfolding subresultant_def ..
also have "(-1)^(?k * ?n) * … = [: (-1)^(?k * ?n) :] * … " by (unfold hom_distribs, simp)
also have "… = smult ((-1)^(?k * ?n)) (subresultant J g f)" by simp
finally show ?thesis .
qed
lemma subresultant_smult_right:assumes "(c :: 'a :: {comm_ring_1, semiring_no_zero_divisors}) ≠ 0"
shows "subresultant J f (smult c g) = smult (c ^ (degree f - J)) (subresultant J f g)"
unfolding subresultant_swap[of _ f] subresultant_smult_left[OF assms]
degree_smult_eq using assms by (simp add: ac_simps)
lemma coeff_subresultant: "coeff (subresultant J F G) l =
(if degree F - J + (degree G - J) = 0 ∧ l ≠ 0 then 0 else det (subresultant'_mat J l F G))"
proof (cases "degree F - J + (degree G - J) = 0")
case True
show ?thesis unfolding subresultant_def subresultant_mat_def subresultant'_mat_def Let_def True
by simp
next
case False
let ?n = "degree F - J + (degree G - J)"
define n where "n = ?n"
from False have n: "n ≠ 0" unfolding n_def by auto
hence id: "{0..<n} = insert (n - 1) {0..< (n - 1)}" by (cases n, auto)
have idn: "(x = x) = True" for x :: nat by simp
let ?M = "subresultant_mat J F G"
define M where "M = ?M"
let ?L = "subresultant'_mat J l F G"
define L where "L = ?L"
{
fix p
assume p: "p permutes {0..<n}"
from n p have n1: "n - 1 < n" "p (n - 1) < n" by auto
have "coeff_int (∏i = 0..<n. M $$ (i, p i)) l =
(∏i = 0 ..< (n - 1). coeff_int (M $$ (i, p i)) 0) * coeff_int (M $$ (n - 1, p (n - 1))) l"
unfolding id
proof (rule coeff_int_prod_const, (auto)[2])
fix i
assume "i ∈ {0 ..< n - 1}"
with p have i: "i ≠ n - 1" and "i < n" "p i < n" by (auto simp: n_def)
note id = subresultant_index_mat[OF this(2-3)[unfolded n_def], folded M_def n_def]
show "degree (M $$ (i, p i)) = 0" unfolding id Let_def using i
by (simp split: if_splits)
qed
also have "(∏i = 0 ..< (n - 1). coeff_int (M $$ (i, p i)) 0)
= (∏i = 0 ..< (n - 1). L $$ (i, p i))"
proof (rule prod.cong[OF refl])
fix i
assume "i ∈ {0 ..< n - 1}"
with p have i: "i ≠ n - 1" and ii: "i < n" "p i < n" by (auto simp: n_def)
note id = subresultant_index_mat[OF this(2-3)[unfolded n_def], folded M_def n_def]
note id' = L_def[unfolded subresultant'_mat_def Let_def, folded n_def] index_mat[OF ii]
show "coeff_int (M $$ (i, p i)) 0 = L $$ (i, p i)"
unfolding id id' split using i proof (simp add: if_splits Let_def)
qed
qed
also have "coeff_int (M $$ (n - 1, p (n - 1))) l =
(if p (n - 1) < degree G - J then
coeff_int (monom 1 (degree G - J - 1 - p (n - 1)) * F) l
else coeff_int (monom 1 (degree F - J - 1 - (p (n - 1) - (degree G - J))) * G) l)"
using subresultant_index_mat[OF n1[unfolded n_def], folded M_def n_def, unfolded idn if_True Let_def]
by simp
also have "… = (if p (n - 1) < degree G - J
then coeff_int F (int l - int (degree G - J - 1 - p (n - 1)))
else coeff_int G (int l - int (degree F - J - 1 - (p (n - 1) - (degree G - J)))))"
unfolding coeff_int_monom_mult by simp
also have "… = (if p (n - 1) < degree G - J
then coeff_int F (int l - int (degree G - J - 1) + p (n - 1))
else coeff_int G (int l - int (degree F - J - 1) + (p (n - 1) - (degree G - J))))"
proof (cases "p (n - 1) < degree G - J")
case True
hence "int (degree G - J - 1 - p (n - 1)) = int (degree G - J - 1) - p (n - 1)" by simp
hence id: "int l - int (degree G - J - 1 - p (n - 1)) = int l - int (degree G - J - 1) + p (n - 1)" by simp
show ?thesis using True unfolding id by simp
next
case False
from n1 False have "degree F - J - 1 ≥ p (n - 1) - (degree G - J)"
unfolding n_def by linarith
hence "int (degree F - J - 1 - (p (n - 1) - (degree G - J))) = int (degree F - J - 1) - (p (n - 1) - (degree G - J))"
by linarith
hence id: "int l - int (degree F - J - 1 - (p (n - 1) - (degree G - J))) =
int l - int (degree F - J - 1) + (p (n - 1) - (degree G - J))" by simp
show ?thesis unfolding id using False by simp
qed
also have "… = L $$ (n - 1, p (n - 1))"
unfolding L_def subresultant'_mat_def Let_def n_def[symmetric] using n1 by simp
also have "(∏i = 0..<n - 1. L $$ (i, p i)) * … = (∏i = 0..<n. L $$ (i, p i))"
unfolding id by simp
finally have "coeff_int (∏i = 0..<n. M $$ (i, p i)) (int l) = (∏i = 0..<n. L $$ (i, p i))" .
} note * = this
have "coeff_int (subresultant J F G) l =
(∑p∈{p. p permutes {0..<n}}. signof p * coeff_int (∏i = 0..<n. M $$ (i, p i)) l)"
unfolding subresultant_def det_def subresultant_mat_dim idn if_True n_def[symmetric] M_def
coeff_int_sum coeff_int_signof_mult by simp
also have "… = (∑p∈{p. p permutes {0..<n}}. signof p * (∏i = 0..<n. L $$ (i, p i)))"
by (rule sum.cong[OF refl], insert *, simp)
also have "… = det L"
proof -
have id: "dim_row (subresultant'_mat J l F G) = n"
"dim_col (subresultant'_mat J l F G) = n" unfolding subresultant'_mat_def Let_def n_def
by auto
show ?thesis unfolding det_def L_def id by simp
qed
finally show ?thesis unfolding L_def coeff_int_def using False by auto
qed
lemma subresultant'_zero_ge: assumes "(degree f - J) + (degree g - J) ≠ 0" and "k ≥ degree f + (degree g - J)"
shows "det (subresultant'_mat J k f g) = 0"
proof -
obtain dg where dg: "degree g - J = dg" by simp
obtain df where df: "degree f - J = df" by simp
obtain ddf where ddf: "degree f = ddf" by simp
note * = assms(2)[unfolded ddf dg] assms(1)
define M where "M = (λ i j. if j < dg
then coeff_int f (degree f - int i + int j)
else coeff_int g (degree g - int i + int (j - dg)))"
let ?M = "subresultant'_mat J k f g"
have M: "det ?M = det (mat (df + dg) (df + dg)
(λ(i, j).
if i = df + dg - 1 then
if j < dg
then coeff_int f (int k - int (dg - 1) + int j)
else coeff_int g (int k - int (df - 1) + int (j - dg))
else M i j))" (is "_ = det ?N")
unfolding subresultant'_mat_def Let_def M_def
by (rule arg_cong[of _ _ det], rule eq_matI, auto simp: df dg)
also have "?N = mat (df + dg) (df + dg)
(λ(i, j).
if i = df + dg - 1 then 0
else M i j)"
by (rule cong_mat[OF refl refl], unfold split, rule if_cong[OF refl _ refl],
auto simp add: coeff_int_def df dg ddf intro!: coeff_eq_0, insert *(1),
unfold ddf[symmetric] dg[symmetric] df[symmetric], linarith+)
also have "… = mat⇩r (df + dg) (df + dg) (λi. if i = df + dg - 1 then 0⇩v (df + dg) else
vec (df + dg) (λ j. M i j))"
by (rule eq_matI, auto)
also have "det … = 0"
by (rule det_row_0, insert *, auto simp: df[symmetric] dg[symmetric] ddf[symmetric])
finally show ?thesis .
qed
lemma subresultant'_zero_lt: assumes
J: "J ≤ degree f" "J ≤ degree g" "J < k"
and k: "k < degree f + (degree g - J)"
shows "det (subresultant'_mat J k f g) = 0"
proof -
obtain dg where dg: "dg = degree g - J" by simp
obtain df where df: "df = degree f - J" by simp
note * = assms[folded df dg]
define M where "M = (λ i j. if j < dg
then coeff_int f (degree f - int i + int j)
else coeff_int g (degree g - int i + int (j - dg)))"
define N where "N = (λ j. if j < dg
then coeff_int f (int k - int (dg - 1) + int j)
else coeff_int g (int k - int (df - 1) + int (j - dg)))"
let ?M = "subresultant'_mat J k f g"
have M: "?M = mat (df + dg) (df + dg)
(λ(i, j).
if i = df + dg - 1 then N j
else M i j)"
unfolding subresultant'_mat_def Let_def
by (rule eq_matI, auto simp: df dg M_def N_def)
also have "… = mat (df + dg) (df + dg)
(λ(i, j).
if i = df + dg - 1 then N j
else if i = degree f + dg - 1 - k then N j else M i j)" (is "_ = ?N")
unfolding N_def
by (rule cong_mat[OF refl refl], unfold split, rule if_cong[OF refl refl], unfold M_def N_def,
insert J k, auto simp: df dg intro!: arg_cong[of _ _ "coeff_int _"])
finally have id: "?M = ?N" .
have deg: "degree f + dg - 1 - k < df + dg" "df + dg - 1 < df + dg"
using k J unfolding df dg by auto
have id: "row ?M (degree f + dg - 1 - k) = row ?M (df + dg - 1)"
unfolding arg_cong[OF id, of row]
by (rule eq_vecI, insert deg, auto)
show ?thesis
by (rule det_identical_rows[OF _ _ _ _ id, of "df + dg"], insert deg assms,
auto simp: subresultant'_mat_def Let_def df dg)
qed
lemma subresultant'_mat_sylvester_mat: "transpose_mat (subresultant'_mat 0 0 f g) = sylvester_mat f g"
proof -
obtain dg where dg: "degree g = dg" by simp
obtain df where df: "degree f = df" by simp
let ?M = "transpose_mat (subresultant'_mat 0 0 f g)"
let ?n = "degree f + degree g"
have dim: "dim_row ?M = ?n" "dim_col ?M = ?n" by (auto simp: subresultant'_mat_def Let_def)
show ?thesis
proof (rule eq_matI, unfold sylvester_mat_dim dim df dg, goal_cases)
case ij: (1 i j)
have "?M $$ (i,j) = (if i < dg
then if j = df + dg - 1
then coeff_int f (- int (dg - 1) + int i)
else coeff_int f (int df - int j + int i)
else if j = df + dg - 1
then coeff_int g (- int (df - 1) + int (i - dg))
else coeff_int g (int dg - int j + int (i - dg)))"
using ij unfolding subresultant'_mat_def Let_def by (simp add: if_splits df dg)
also have "… = (if i < dg
then coeff_int f (int df - int j + int i)
else coeff_int g (int dg - int j + int (i - dg)))"
proof -
have cong: "(b ⟹ x = z) ⟹ (¬ b ⟹ y = z) ⟹ (if b then coeff_int f x else coeff_int f y) = coeff_int f z"
for b x y z and f :: "'a poly" by auto
show ?thesis
by (rule if_cong[OF refl cong cong], insert ij, auto)
qed
also have "… = sylvester_mat f g $$ (i,j)"
proof -
have *: "i ≤ j ⟹ j - i ≤ df ⟹ nat (int df - int j + int i) = df - (j - i)" for j i df
by simp
show ?thesis unfolding sylvester_index_mat[OF ij[folded df dg]] df dg
proof (rule if_cong[OF refl])
assume i: "i < dg"
have "int df - int j + int i < 0 ⟶ ¬ j - i ≤ df" by auto
thus "coeff_int f (int df - int j + int i) = (if i ≤ j ∧ j - i ≤ df then coeff f (df + i - j) else 0)"
using i ij by (simp add: coeff_int_def *, intro impI coeff_eq_0[of f, unfolded df], linarith)
next
assume i: "¬ i < dg"
hence **: "i - dg ≤ j ⟹ dg - (j + dg - i) = i - j" using ij by linarith
have "int dg - int j + int (i - dg) < 0 ⟶ ¬ j ≤ i" by auto
thus "coeff_int g (int dg - int j + int (i - dg)) = (if i - dg ≤ j ∧ j ≤ i then coeff g (i - j) else 0)"
using ij i by (simp add: coeff_int_def * **, intro impI coeff_eq_0[of g, unfolded dg], linarith)
qed
qed
finally show ?case .
qed auto
qed
lemma coeff_subresultant_0_0_resultant: "coeff (subresultant 0 f g) 0 = resultant f g"
proof -
let ?M = "transpose_mat (subresultant'_mat 0 0 f g)"
have "det (subresultant'_mat 0 0 f g) = det ?M"
by (subst det_transpose, auto simp: subresultant'_mat_def Let_def)
also have "?M = sylvester_mat f g"
by (rule subresultant'_mat_sylvester_mat)
finally show ?thesis by (simp add: coeff_subresultant resultant_def)
qed
lemma subresultant_zero_ge: assumes "k ≥ degree f + (degree g - J)"
and "(degree f - J) + (degree g - J) ≠ 0"
shows "coeff (subresultant J f g) k = 0"
unfolding coeff_subresultant
by (subst subresultant'_zero_ge[OF assms(2,1)], simp)
lemma subresultant_zero_lt: assumes "k < degree f + (degree g - J)"
and "J ≤ degree f" "J ≤ degree g" "J < k"
shows "coeff (subresultant J f g) k = 0"
unfolding coeff_subresultant
by (subst subresultant'_zero_lt[OF assms(2,3,4,1)], simp)
lemma subresultant_resultant: "subresultant 0 f g = [: resultant f g :]"
proof (cases "degree f + degree g = 0")
case True
thus ?thesis unfolding subresultant_def subresultant_mat_def resultant_def Let_def
sylvester_mat_def sylvester_mat_sub_def
by simp
next
case 0: False
show ?thesis
proof (rule poly_eqI)
fix k
show "coeff (subresultant 0 f g) k = coeff [:resultant f g:] k"
proof (cases "k = 0")
case True
thus ?thesis using coeff_subresultant_0_0_resultant[of f g] by auto
next
case False
hence "0 < k ∧ k < degree f + degree g ∨ k ≥ degree f + degree g" by auto
thus ?thesis using subresultant_zero_ge[of f g 0 k] 0
subresultant_zero_lt[of k f g 0] 0 False by (cases k, auto)
qed
qed
qed
lemma (in inj_comm_ring_hom) subresultant_hom:
"map_poly hom (subresultant J f g) = subresultant J (map_poly hom f) (map_poly hom g)"
proof -
note d = subresultant_mat_def Let_def
interpret p: map_poly_inj_comm_ring_hom hom..
show ?thesis unfolding subresultant_def unfolding p.hom_det[symmetric]
proof (rule arg_cong[of _ _ det])
show "p.mat_hom (subresultant_mat J f g) =
subresultant_mat J (map_poly hom f) (map_poly hom g)"
proof (rule eq_matI, goal_cases)
case (1 i j)
hence ij: "i < degree f - J + (degree g - J)" "j < degree f - J + (degree g - J)"
unfolding d degree_map_poly by auto
show ?case
by (auto simp add: coeff_int_def d map_mat_def index_mat(1)[OF ij] hom_distribs)
qed (auto simp: d)
qed
qed
text ‹We now derive properties of the resultant via the connection to subresultants.›
lemma resultant_smult_left: assumes "(c :: 'a :: idom) ≠ 0"
shows "resultant (smult c f) g = c ^ degree g * resultant f g"
unfolding coeff_subresultant_0_0_resultant[symmetric] subresultant_smult_left[OF assms] coeff_smult
by simp
lemma resultant_smult_right: assumes "(c :: 'a :: idom) ≠ 0"
shows "resultant f (smult c g) = c ^ degree f * resultant f g"
unfolding coeff_subresultant_0_0_resultant[symmetric] subresultant_smult_right[OF assms] coeff_smult
by simp
lemma resultant_swap: "resultant f g = (-1)^(degree f * degree g) * (resultant g f)"
unfolding coeff_subresultant_0_0_resultant[symmetric]
unfolding arg_cong[OF subresultant_swap[of 0 f g], of "λ x. coeff x 0"] coeff_smult by simp
text ‹The following equations are taken from Brown-Traub ``On Euclid's Algorithm and
the Theory of Subresultant'' (BT)›
lemma fixes F B G H :: "'a :: idom poly" and J :: nat
defines df: "df ≡ degree F"
and dg: "dg ≡ degree G"
and dh: "dh ≡ degree H"
and db: "db ≡ degree B"
defines
n: "n ≡ (df - J) + (dg - J)"
and f: "f ≡ coeff_int F"
and b: "b ≡ coeff_int B"
and g: "g ≡ coeff_int G"
and h: "h ≡ coeff_int H"
assumes FGH: "F + B * G = H"
and dfg: "df ≥ dg"
and choice: "dg > dh ∨ H = 0 ∧ F ≠ 0 ∧ G ≠ 0"
shows BT_eq_18: "subresultant J F G = smult ((-1)^((df - J) * (dg - J))) (det (mat n n
(λ (i,j).
if j < df - J
then if i = n - 1 then monom 1 ((df - J) - 1 - j) * G
else [:g (int dg - int i + int j):]
else if i = n - 1 then monom 1 ((dg - J) - 1 - (j - (df - J))) * H
else [:h (int df - int i + int (j - (df - J))):])))"
(is "_ = smult ?m1 ?right")
and BT_eq_19: "dh ≤ J ⟹ J < dg ⟹ subresultant J F G = smult (
(-1)^((df - J) * (dg - J)) * lead_coeff G ^ (df - J) * coeff H J ^ (dg - J - 1)) H"
(is "_ ⟹ _ ⟹ _ = smult (_ * ?G * ?H) H")
and BT_lemma_1_12: "J < dh ⟹ subresultant J F G = smult (
(-1)^((df - J) * (dg - J)) * lead_coeff G ^ (df - dh)) (subresultant J G H)"
and BT_lemma_1_13': "J = dh ⟹ dg > dh ∨ H ≠ 0 ⟹ subresultant dh F G = smult (
(-1)^((df - dh) * (dg - dh)) * lead_coeff G ^ (df - dh) * lead_coeff H ^ (dg - dh - 1)) H"
and BT_lemma_1_14: "dh < J ⟹ J < dg - 1 ⟹ subresultant J F G = 0"
and BT_lemma_1_15': "J = dg - 1 ⟹ dg > dh ∨ H ≠ 0 ⟹ subresultant (dg - 1) F G = smult (
(-1)^(df - dg + 1) * lead_coeff G ^ (df - dg + 1)) H"
proof -
define dfj where "dfj = df - J"
define dgj where "dgj = dg - J"
note d = df dg dh db
have F0: "F ≠ 0" using dfg choice df by auto
have G0: "G ≠ 0" using choice dg by auto
have dgh: "dg ≥ dh" using choice unfolding dh by auto
have B0: "B ≠ 0" using FGH dfg dgh choice F0 G0 unfolding d by auto
have dfh: "df ≥ dh" using dfg dgh by auto
have "df = degree (B * G)"
proof (cases "H = 0")
case False
with choice dfg have dfh: "df > dh" by auto
show ?thesis using dfh[folded arg_cong[OF FGH, of degree, folded dh]] choice
unfolding df by (metis ‹degree (F + B * G) < df› degree_add_eq_left degree_add_eq_right df nat_neq_iff)
next
case True
have "F = - B * G" using arg_cong[OF FGH[unfolded True], of "λ x. x - B * G"] by auto
thus ?thesis using F0 G0 B0 unfolding df by simp
qed
hence dfbg: "df = db + dg" using degree_mult_eq[OF B0 G0] by (simp add: d)
hence dbfg: "db = df - dg" by simp
let ?dfj = "df - J"
let ?dgj = "dg - J"
have norm: "?dgj + ?dfj = ?dfj + ?dgj" by simp
let ?bij = "λ i j. b (db - int i + int (j - dfj))"
let ?M = "mat n n (λ (i,j). if i = j then 1 else if j < dfj then 0 else if i < j
then [: ?bij i j :] else 0)"
let ?GF = "λ i j.
if j < dfj
then if i = n - 1 then monom 1 (dfj - 1 - j) * G
else [:g (int dg - int i + int j):]
else if i = n - 1 then monom 1 (dgj - 1 - (j - dfj)) * F
else [:f (int df - int i + int (j - dfj)):]"
let ?G_F = "mat n n (λ (i,j). ?GF i j)"
let ?GH = "λ i j.
if j < dfj
then if i = n - 1 then monom 1 (dfj - 1 - j) * G
else [:g (int dg - int i + int j):]
else if i = n - 1 then monom 1 (dgj - 1 - (j - dfj)) * H
else [:h (int df - int i + int (j - dfj)):]"
let ?G_H = "mat n n (λ (i,j). ?GH i j)"
have hfg: "h i = f i + coeff_int (B * G) i" for i
unfolding FGH[symmetric] f g h unfolding coeff_int_def by simp
have dM1: "det ?M = 1"
by (subst det_upper_triangular, (auto)[2], subst prod_list_diag_prod, auto)
have "subresultant J F G = smult ?m1 (subresultant J G F)"
unfolding subresultant_swap[of _ F] d by simp
also have "subresultant J G F = det ?G_F"
unfolding subresultant_def n norm subresultant_mat_def g f Let_def d[symmetric] dfj_def dgj_def by simp
also have "… = det (?G_F * ?M)"
by (subst det_mult[of _ n], unfold dM1, auto)
also have "?G_F * ?M = ?G_H"
proof (rule eq_matI, unfold dim_col_mat dim_row_mat)
fix i j
assume i: "i < n" and j: "j < n"
have "(?G_F * ?M) $$ (i,j) = row (?G_F) i ∙ col ?M j"
using i j by simp
also have "… = ?GH i j"
proof (cases "j < dfj")
case True
have id: "col ?M j = unit_vec n j"
by (rule eq_vecI, insert True i j, auto)
show ?thesis unfolding id using True i j by simp
next
case False
define d where "d = j - dfj"
from False have jd: "j = d + dfj" unfolding d_def by auto
hence idj: "{0 ..< j} = {0 ..< dfj} ∪ {dfj ..< dfj + d}" by auto
let ?H = "(if i = n - 1 then monom 1 (dgj - Suc d) * H else [:h (int df - int i + int d):])"
have idr: "?GH i j = ?H" unfolding d_def using jd by auto
let ?bi = "λ i. b (db - int i + int d)"
let ?m = "λ i. if i = j then 1 else if i < j then [:?bij i j:] else 0"
let ?P = "λ k. (?GF i k * ?m k)"
let ?Q = "λ k. ?GF i k * [: ?bi k :]"
let ?G = "λ k. if i = n - 1 then monom 1 (dfj - 1 - k) * G else [:g (int dg - int i + int k):]"
let ?Gb = "λ k. ?G k * [:?bi k:]"
let ?off = "- (int db - int dfj + 1 + int d)"
have off0: "?off ≥ 0" using False dfg j unfolding dfj_def d_def dbfg n by simp
from nat_0_le[OF this]
obtain off where off: "int off = ?off" by blast
have "int off ≤ int dfj" unfolding off by auto
hence "off ≤ dfj" by simp
hence split1: "{0 ..< dfj} = {0 ..< off} ∪ {off ..< dfj}" by auto
have "int off + Suc db ≤ dfj" unfolding off by auto
hence split2: "{off ..< dfj} = {off .. off + db} ∪ {off + Suc db ..< dfj} " by auto
let ?g_b = "λk. (if i = n - 1 then monom 1 k * G else [:g (int dg - int i + int (dfj - Suc k)):]) *
[:b (k - int off):]"
let ?gb = "λk. (if i = n - 1 then monom 1 (k + off) * G else [:g (int dg - int i + int (dfj - Suc k - off)):]) *
[:coeff B k:]"
let ?F = "λ k. if i = n - 1 then monom 1 (dgj - 1 - (k - dfj)) * F
else [:f (int df - int i + int (k - dfj)):]"
let ?Fb = "λ k. ?F k * [:?bi k:]"
let ?Pj = "if i = n - 1 then monom 1 (dgj - Suc d) * F else [:f (int df - int i + int d):]"
from False have id: "col ?M j = vec n ?m"
using j i by (intro eq_vecI, auto)
have "row ?G_F i ∙ col ?M j = sum ?P {0 ..< n}"
using i j unfolding id by (simp add: scalar_prod_def)
also have "{0 ..< n} = {0 ..< j} ∪ {j} ∪ {Suc j ..< n}" using j by auto
also have "sum ?P … = sum ?P {0 ..< j} + ?P j + sum ?P {Suc j ..< n}"
by (simp add: sum.union_disjoint)
also have "sum ?P {Suc j ..< n} = 0" by (rule sum.neutral, auto)
also have "?P j = ?Pj"
unfolding d_def using jd by simp
also have "sum ?P {0 ..< j} = sum ?Q {0 ..< j}"
by (rule sum.cong[OF refl], unfold d_def, insert jd, auto)
also have "sum ?Q {0 ..< j} = sum ?Q {0 ..< dfj} + sum ?Q {dfj ..< dfj+d}" unfolding idj
by (simp add: sum.union_disjoint)
also have "sum ?Q {0 ..< dfj} = sum ?Gb {0 ..< dfj}"
by (rule sum.cong, auto)
also have "sum ?Q {dfj ..< dfj+d} = sum ?Fb {dfj ..< dfj+d}"
by (rule sum.cong, auto)
also have "… = 0"
proof (rule sum.neutral, intro ballI)
fix k
assume k: "k ∈ {dfj ..< dfj+d}"
hence k: "db + d < k" using k j False unfolding n db[symmetric] dfbg dfj_def d_def by auto
let ?k = "(int db - int k + int d)"
have "?k < 0" using k by auto
hence "b ?k = 0" unfolding b by (intro coeff_int_eq_0, auto)
thus "?Fb k = 0" by simp
qed
also have "sum ?Gb {0 ..< dfj} = sum ?g_b {0 ..< dfj}"
proof (rule sum.reindex_cong[of "λ k. dfj - Suc k"], (auto simp: inj_on_def off)[2], goal_cases)
case (1 k)
hence "k = dfj - (Suc (dfj - Suc k))" and "(dfj - Suc k) ∈ {0..<dfj}" by auto
thus ?case by blast
next
case (2 k)
hence [simp]: "dfj - Suc (dfj - Suc k) = k"
"int db - int (dfj - Suc k) + int d = int k - off" by (auto simp: off)
show ?case by auto
qed
also have "… = sum ?g_b {0 ..< off} + sum ?g_b {off ..< dfj}" unfolding split1
by (simp add: sum.union_disjoint)
also have "sum ?g_b {0 ..< off} = 0"
by (rule sum.neutral, intro ballI, auto simp: b coeff_int_def)
also have "sum ?g_b {off ..< dfj} = sum ?g_b {off .. off + db} + sum ?g_b {off + Suc db ..< dfj}"
unfolding split2 by (rule sum.union_disjoint, auto)
also have "sum ?g_b {off + Suc db ..< dfj} = 0"
proof (rule sum.neutral, intro ballI, goal_cases)
case (1 k)
hence "b (int k - int off) = 0" unfolding b db
by (intro coeff_int_eq_0, auto)
thus ?case by simp
qed
also have "sum ?g_b {off .. off + db} = sum ?gb {0 .. db}"
using sum.atLeastAtMost_shift_bounds [of ?g_b 0 off db]
by (auto intro: sum.cong simp add: b ac_simps)
finally have id: "row ?G_F i ∙ col ?M j - ?H = ?Pj + sum ?gb {0 .. db} - ?H"
(is "_ = ?E")
by (simp add: ac_simps)
define E where "E = ?E"
let ?b = "coeff B"
have Bsum: "(∑k = 0..db. monom (?b k) k) = B" unfolding db
using atMost_atLeast0 poly_as_sum_of_monoms by auto
have "E = 0"
proof (cases "i = n - 1")
case i_n: False
hence id: "(i = n - 1) = False" by simp
with i have i: "i < n - 1" by auto
let ?ii = "int df - int i + int d"
have "?thesis = ([:f ?ii:] +
(∑k = 0..db.
[:g (int dg - int i + int (dfj - Suc k - off)):] * [:?b k:]) -
[:h ?ii:] = 0)" (is "_ = (?e = 0)") unfolding E_def id if_False by simp
also have "?e = [: f ?ii +
(∑k = 0..db.
g (int dg - int i + int (dfj - Suc k - off)) * ?b k) -
h ?ii:]" (is "_ = [: ?e :]")
proof (rule poly_eqI, goal_cases)
case (1 n)
show ?case unfolding coeff_diff coeff_add coeff_sum coeff_const
by (cases n, auto simp: ac_simps)
qed
also have "[: ?e :] = 0 ⟷ ?e = 0" by simp
also have "?e = (∑k = 0..db. g (int dg - int i + int (dfj - Suc k - off)) * ?b k)
- coeff_int (B * G) ?ii"
unfolding hfg by simp
also have "(B * G) = (∑k = 0..db. monom (?b k) k) * G" unfolding Bsum by simp
also have "… = (∑k = 0..db. monom (?b k) k * G)" by (rule sum_distrib_right)
also have "coeff_int … ?ii = (∑k = 0..db. g (?ii - k) * ?b k)"
unfolding coeff_int_sum coeff_int_monom_mult g by (simp add: ac_simps)
also have "… = (∑k = 0..db. g (int dg - int i + int (dfj - Suc k - off)) * ?b k)"
proof (rule sum.cong[OF refl], goal_cases)
case (1 k)
hence "k ≤ db" by simp
hence id: "int dg - int i + int (dfj - Suc k - off) = ?ii - k"
using False i j off dfg
unfolding dbfg d_def dfj_def n by linarith
show ?case unfolding id ..
qed
finally show ?thesis by simp
next
case True
let ?jj = "dgj - Suc d"
have zero: "int off - (dgj - Suc d) = 0" using dfg False j unfolding off dbfg dfj_def d_def dgj_def n
by linarith
from True have "E = monom 1 ?jj * F + (∑k = 0.. db.
monom 1 (k + off) * G * [: ?b k :]) - monom 1 ?jj * H"
(is "_ = ?A + ?sum - ?mon") unfolding id E_def by simp
also have "?mon = monom 1 ?jj * F + monom 1 ?jj * (B * G)"
unfolding FGH[symmetric] by (simp add: ring_distribs)
also have "?A + ?sum - … = ?sum - (monom 1 ?jj * G) * B" (is "_ = _ - ?GB * B") by simp
also have "?sum = (∑k = 0..db.
(monom 1 ?jj * G) * (monom 1 (k + off - ?jj) * [: ?b k :]))"
proof (rule sum.cong[OF refl], goal_cases)
case (1 k)
let ?one = "1 :: 'a"
have "int off ≥ int ?jj" using j False i True
unfolding off d_def dfj_def dgj_def dfbg n by linarith
hence "k + off = ?jj + (k + off - ?jj)" by linarith
hence id: "monom ?one (k + off) = monom (1 * 1) (?jj + (k + off - ?jj))" by simp
show ?case unfolding id[folded mult_monom] by (simp add: ac_simps)
qed
also have "… = (monom 1 ?jj * G) * (∑k = 0..db. monom 1 (k + off - ?jj) * [:?b k:])"
(is "_ = _ * ?sum")
unfolding sum_distrib_left ..
also have "… - (monom 1 ?jj * G) * B = (monom 1 ?jj * G) * (?sum - B)" by (simp add: ring_distribs)
also have "?sum = (∑k = 0..db. monom 1 k * [:?b k:])"
by (rule sum.cong[OF refl], insert zero, auto)
also have "… = (∑k = 0..db. monom (?b k) k)"
by (rule sum.cong[OF refl], rule poly_eqI, auto)
also have "… = B" unfolding Bsum ..
finally show ?thesis by simp
qed
from id[folded E_def, unfolded this]
show ?thesis using False unfolding d_def by simp
qed
also have "… = ?G_H $$ (i,j)" using i j by simp
finally show "(?G_F * ?M) $$ (i,j) = ?G_H $$ (i,j)" .
qed auto
finally show eq_18: "subresultant J F G = smult ?m1 (det ?G_H)" unfolding dfj_def dgj_def .
{
fix i j
assume ij: "i < j" and j: "j < n"
with dgh have "int dg - int i + int j > int dg" by auto
hence "g (int dg - int i + int j) = 0" unfolding g dg by (intro coeff_int_eq_0, auto)
} note g0 = this
{
assume *: "dh ≤ J" "J < dg"
have n_dfj: "n > dfj" using * unfolding n dfj_def by auto
note eq_18
also have "det ?G_H = prod_list (diag_mat ?G_H)"
proof (rule det_lower_triangular[of n])
fix i j
assume ij: "i < j" and j: "j < n"
from ij j have if_e: "i = n - 1 ⟷ False" by auto
have "?G_H $$ (i,j) = ?GH i j" using ij j by auto
also have "… = 0"
proof (cases "j < dfj")
case True
with True g0[OF ij j] show ?thesis unfolding if_e by simp
next
case False
have "h (int df - int i + int (j - dfj)) = 0" unfolding h
by (rule coeff_int_eq_0, insert False * ij j dfg, unfold dfj_def dh[symmetric], auto)
with False show ?thesis unfolding if_e by auto
qed
finally show "?G_H $$ (i,j) = 0" .
qed auto
also have "… = (∏i = 0..<n. ?GH i i)"
by (subst prod_list_diag_prod, simp)
also have "{0 ..< n} = {0 ..< dfj} ∪ {dfj ..< n}" unfolding n dfj_def by auto
also have "prod (λ i. ?GH i i) … = prod (λ i. ?GH i i) {0 ..< dfj} * prod (λ i. ?GH i i) {dfj ..< n}"
by (simp add: prod.union_disjoint)
also have "prod (λ i. ?GH i i) {0 ..< dfj} = prod (λ i. [: lead_coeff G :]) {0 ..< dfj}"
proof -
show ?thesis
by (rule prod.cong[OF refl], insert n_dfj, auto simp: g coeff_int_def dg)
qed
also have "… = [: (lead_coeff G)^dfj :]" by (simp add: poly_const_pow)
also have "{dfj ..< n} = {dfj ..< n-1} ∪ {n - 1}" using n_dfj by auto
also have "prod (λ i. ?GH i i) … = prod (λ i. ?GH i i) {dfj ..< n-1} * ?GH (n - 1) (n - 1)"
by (simp add: prod.union_disjoint)
also have "?GH (n - 1) (n - 1) = H"
proof -
have "dgj - 1 - (n - 1 - dfj) = 0" using n_dfj unfolding dgj_def dfj_def n by auto
with n_dfj show ?thesis by auto
qed
also have "prod (λ i. ?GH i i) {dfj ..< n-1} = prod (λ i. [:h (int df - dfj):]) {dfj ..< n-1}"
by (rule prod.cong[OF refl], auto intro!: arg_cong[of _ _ h])
also have "… = [: h (int df - dfj) ^ (n - 1 - dfj) :]"
unfolding prod_constant by (simp add: poly_const_pow)
also have "n - 1 - dfj = dg - J - 1" unfolding n dfj_def by simp
also have "int df - dfj = J" using * dfg unfolding dfj_def by auto
also have "h J = coeff H J" unfolding h coeff_int_def by simp
finally show "subresultant J F G = smult (?m1 * ?G * ?H) H" by (simp add: dfj_def ac_simps)
} note eq_19 = this
{
assume J: "J < dh"
define dhj where "dhj = dh - J"
have n_add: "n = (df - dh) + (dhj + dgj)" unfolding dhj_def dgj_def n using J dfg dgh by auto
let ?split = "split_block ?G_H (df - dh) (df - dh)"
have dim: "dim_row ?G_H = (df - dh) + (dhj + dgj)"
"dim_col ?G_H = (df - dh)