Theory Dichotomous_Lazard
section ‹Dichotomous Lazard›
text ‹This theory contains Lazard's optimization in the computation of
the subresultant PRS as described by Ducos \<^cite>‹‹Section 2› in "Ducos"›.›
theory Dichotomous_Lazard
imports
"HOL-Computational_Algebra.Polynomial_Factorial"
begin
lemma power_fract[simp]: "(Fract a b)^n = Fract (a^n) (b^n)"
by (induct n, auto simp: fract_collapse)
lemma range_to_fract_dvd_iff: assumes b: "b ≠ 0"
shows "Fract a b ∈ range to_fract ⟷ b dvd a"
proof
assume "b dvd a" then obtain c where a: "a = b * c" unfolding dvd_def by auto
have "Fract a b = Fract c 1" using b unfolding a by (simp add: eq_fract)
thus "Fract a b ∈ range to_fract" unfolding to_fract_def by auto
next
assume "Fract a b ∈ range to_fract"
then obtain c where "Fract a b = Fract c 1" unfolding to_fract_def by auto
hence "a = b * c" using b by (simp add: eq_fract)
thus "b dvd a" ..
qed
lemma Fract_cases_coprime [cases type: fract]:
fixes q :: "'a :: factorial_ring_gcd fract"
obtains (Fract) a b where "q = Fract a b" "b ≠ 0" "coprime a b"
proof -
obtain a b where q: "q = Fract a b" and b0: "b ≠ 0" by (cases q, auto)
define g where g: "g = gcd a b"
define A where A: "A = a div g"
define B where B: "B = b div g"
have a: "a = A * g" unfolding A g by simp
have b: "b = B * g" unfolding B g by simp
from b0 b have 0: "B ≠ 0" by auto
have q: "q = Fract A B" unfolding q a b
by (subst eq_fract, auto simp: b0 0 g)
have cop: "coprime A B" unfolding A B g using b0
by (simp add: div_gcd_coprime)
assume "⋀a b. q = Fract a b ⟹ b ≠ 0 ⟹ coprime a b ⟹ thesis"
from this[OF q 0 cop] show ?thesis .
qed
lemma to_fract_power_le: fixes a :: "'a :: factorial_ring_gcd fract"
assumes no_fract: "a * b ^ e ∈ range to_fract"
and a: "a ∈ range to_fract"
and le: "f ≤ e"
shows "a * b ^ f ∈ range to_fract"
proof -
obtain bn bd where b: "b = Fract bn bd" and bd: "bd ≠ 0" and copb: "coprime bn bd" by (cases b, auto)
obtain an where a: "a = Fract an 1" using a unfolding to_fract_def by auto
have id: "a * b ^ e = Fract (an * bn^e) (bd^e)" unfolding a b power_fract mult_fract by simp
have 0: "bd^e ≠ 0" for e using bd by auto
from no_fract[unfolded id range_to_fract_dvd_iff[OF 0]] have dvd: "bd ^ e dvd an * bn ^ e" .
from copb have copb: "coprime (bd ^ e) (bn ^ e)" for e
by (simp add: ac_simps)
from dvd copb [of e] bd have "bd ^ e dvd an"
by (simp add: coprime_dvd_mult_left_iff)
hence "bd ^ f dvd an" using le by (rule power_le_dvd)
hence dvd: "bd ^ f dvd an * bn ^ f" by simp
from le obtain g where e: "e = f + g" using le_Suc_ex by blast
have id': "a * b ^ f = Fract (an * bn^f) (bd^f)" unfolding a b power_fract mult_fract by simp
show ?thesis unfolding id' range_to_fract_dvd_iff[OF 0] by (rule dvd)
qed
lemma div_divide_to_fract: assumes "x ∈ range to_fract"
and "x = (y :: 'a :: idom_divide fract) / z"
and "x' = y' div z'"
and "y = to_fract y'" "z = to_fract z'"
shows "x = to_fract x'"
proof (cases "z' = 0")
case True
thus ?thesis using assms by auto
next
case False
from assms obtain r where "to_fract y' / to_fract z' = to_fract r" by auto
thus ?thesis using False assms
by (simp add: eq_fract(1) to_fract_def)
qed
declare Euclidean_Rings.divmod_nat_def [termination_simp]
fun dichotomous_Lazard :: "'a :: idom_divide ⇒ 'a ⇒ nat ⇒ 'a" where
"dichotomous_Lazard x y n = (if n ≤ 1 then if n = 1 then x else 1 else
let (d,r) = Euclidean_Rings.divmod_nat n 2;
rec = dichotomous_Lazard x y d;
recsq = rec * rec div y in
if r = 0 then recsq else recsq * x div y)"
lemma dichotomous_Lazard_main: fixes x :: "'a :: idom_divide"
assumes "⋀ i. i ≤ n ⟹ (to_fract x)^i / (to_fract y)^(i - 1) ∈ range to_fract"
shows "to_fract (dichotomous_Lazard x y n) = (to_fract x)^n / (to_fract y)^(n-1)"
using assms
proof (induct x y n rule: dichotomous_Lazard.induct)
case (1 x y n)
let ?f = to_fract
consider (0) "n = 0" | (1) "n = 1" | (n) "¬ n ≤ 1" by linarith
thus ?case
proof cases
case n
obtain d r where n2: "Euclidean_Rings.divmod_nat n 2 = (d,r)" by force
from Euclidean_Rings.divmod_nat_def[of n 2] n2 have dr: "d = n div 2" "r = n mod 2" by auto
hence r: "r = 0 ∨ r = 1" by auto
define rec where "rec = dichotomous_Lazard x y d"
let ?sq = "rec * rec div y"
have res: "dichotomous_Lazard x y n = (if r = 0 then ?sq else ?sq * x div y)"
unfolding dichotomous_Lazard.simps[of x y n] n2 Let_def rec_def using n by auto
have ndr: "n = d + d + r" unfolding dr by presburger
from ndr r n have d0: "d ≠ 0" by auto
have IH: "?f rec = ?f x ^ d / ?f y ^ (d - 1)"
using 1(1)[OF n refl n2[symmetric] 1(2), folded rec_def] ndr by auto
have "?f (rec * rec) = ?f x ^ d / ?f y ^ (d - 1) * ?f x ^ d / ?f y ^ (d - 1)" using IH by simp
also have "… = ?f x ^ (d + d) / ?f y ^ (d - 1 + (d - 1))" unfolding power_add by simp
also have "d - 1 + (d - 1) = d + d - 2" using d0 by simp
finally have id: "?f (rec * rec) = ?f x ^ (d + d) / ?f y ^ (d + d - 2)" .
let ?dd = "(?f x ^ (d + d) / ?f y ^ (d + d - 2)) / ?f y"
let ?d = "?f x ^ (d + d) / ?f y ^ (d + d - 1)"
have dd: "?dd = ?d" using d0 by (cases d, auto)
have sq: "?f ?sq = ?d" unfolding dd[symmetric]
proof (rule sym, rule div_divide_to_fract[OF _ refl refl id[symmetric] refl], unfold dd)
show "?d ∈ range ?f" by (rule 1(2), insert ndr, auto)
qed
show ?thesis
proof (cases "r = 0")
case True
with res sq show ?thesis unfolding ndr by auto
next
case False
with r have r: "r = 1" by auto
let ?sq' = "?sq * x div y"
from False res have res: "dichotomous_Lazard x y n = ?sq'" by simp
from sq have id: "?f (?sq * x) = ?f x ^ (d + d + r) / ?f y ^ (d + d - 1)"
unfolding r by simp
let ?dd = "(?f x ^ (d + d + r) / ?f y ^ (d + d - 1)) / ?f y"
let ?d = "?f x ^ (d + d + r) / ?f y ^ (d + d + r - 1)"
have dd: "?dd = ?d" using d0 unfolding r by (cases d, auto)
have sq': "?f ?sq' = ?d" unfolding dd[symmetric]
proof (rule sym, rule div_divide_to_fract[OF _ refl refl id[symmetric] refl], unfold dd)
show "?d ∈ range ?f" by (rule 1(2), unfold ndr, auto)
qed
show ?thesis unfolding res sq' unfolding ndr by simp
qed
qed auto
qed
lemma dichotomous_Lazard: fixes x :: "'a :: factorial_ring_gcd"
assumes "(to_fract x)^n / (to_fract y)^(n-1) ∈ range to_fract"
shows "to_fract (dichotomous_Lazard x y n) = (to_fract x)^n / (to_fract y)^(n-1)"
proof (rule dichotomous_Lazard_main)
fix i
assume i: "i ≤ n"
show "to_fract x ^ i / to_fract y ^ (i - 1) ∈ range to_fract"
proof (cases i)
case (Suc j)
have id: "to_fract x ^ i / to_fract y ^ (i - 1) = to_fract x * (to_fract x / to_fract y) ^ j"
unfolding Suc by (simp add: power_divide)
from Suc i have "n ≠ 0" and j: "j ≤ n - 1" by auto
hence idd: "to_fract x * (to_fract x / to_fract y) ^ (n - 1) = (to_fract x)^n / (to_fract y)^(n-1)"
by (cases n, auto simp: power_divide)
show ?thesis unfolding id
by (rule to_fract_power_le[OF _ _ j], unfold idd, insert assms, auto)
next
case 0
have "1 = to_fract 1" by simp
hence "1 ∈ range to_fract" by blast
thus ?thesis using 0 by auto
qed
qed
declare dichotomous_Lazard.simps[simp del]
end