Theory Dichotomous_Lazard

theory Dichotomous_Lazard
imports Polynomial_Factorial
section ‹Dichotomous Lazard›
  
text ‹This theory contains Lazard's optimization in the computation of
  the subresultant PRS as described by Ducos \cite[Section 2]{Ducos}.›
theory Dichotomous_Lazard
imports 
  "HOL-Computational_Algebra.Polynomial_Factorial" (* to_fract *)
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 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) = Divides.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: "Divides.divmod_nat n 2 = (d,r)" by force
    from 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