Theory Explicit_Roots

theory Explicit_Roots
imports Missing_Polynomial Sqrt_Babylonian
(*  
    Author:      René Thiemann 
                 Akihisa Yamada
    License:     BSD
*)
section ‹Explicit Formulas for Roots›

text ‹We provide algorithms which use the explicit formulas to 
  compute the roots of polynomials of degree up to 2. We also define the formula for 
  polynomials of degree 3, but did not (yet) prove anything about it.›  

theory Explicit_Roots
imports   
  Polynomial_Interpolation.Missing_Polynomial
  Sqrt_Babylonian.Sqrt_Babylonian
begin

lemma roots0: assumes p: "p ≠ 0" and p0: "degree p = 0" 
  shows "{x. poly p x = 0} = {}"
  using degree0_coeffs[OF p0] p by auto

definition roots1 :: "'a :: field poly ⇒ 'a" where
  "roots1 p = (- coeff p 0 / coeff p 1)"

lemma roots1: fixes p :: "'a :: field poly"
  assumes p1: "degree p = 1" 
  shows "{x. poly p x = 0} = {roots1 p}"
  using degree1_coeffs[OF p1] unfolding roots1_def 
  by (auto simp: add_eq_0_iff nonzero_neg_divide_eq_eq2)

lemma roots2: fixes p :: "'a :: field_char_0 poly"
  assumes p2: "p = [: c, b, a :]" and a: "a ≠ 0"
  shows "{x. poly p x = 0} = { - ( b / (2 * a)) + e | e. e^2 = ( b / (2 * a))^2 - c/a}" (is "?l = ?r")
proof -
  define b2a where "b2a = b / (2 * a)"
  {
    fix x
    have "(x ∈ ?l) = (x * x * a + x * b + c = 0)" unfolding p2 by (simp add: field_simps)
    also have "… = ((x * x + 2 * x * b2a) + c/a = 0)" using a by (auto simp: b2a_def field_simps)
    also have "x * x + 2 * x * b2a = (x * x + 2 * x * b2a + b2a^2) - b2a^2" by simp
    also have "… = (x + b2a) ^ 2 - b2a ^ 2" 
      by (simp add: field_simps power2_eq_square)
    also have "(… + c / a = 0) = ((x + b2a) ^ 2 = b2a^2 - c/a)" by algebra
    also have "… = (x ∈ ?r)" unfolding b2a_def[symmetric] by (auto simp: field_simps)
    finally have "(x ∈ ?l) = (x ∈ ?r)" .
  }
  thus ?thesis by auto
qed

definition croots2 :: "complex poly ⇒ complex list" where
  "croots2 p = (let a = coeff p 2; b = coeff p 1; c = coeff p 0; b2a = b / (2 * a);
    bac = b2a^2 - c/a;
    e = csqrt bac 
    in
     remdups [- b2a + e, - b2a - e])"

definition complex_rat :: "complex ⇒ bool" where
  "complex_rat x = (Re x ∈ ℚ ∧ Im x ∈ ℚ)"

lemma croots2: assumes "degree p = 2"
  shows "{x. poly p x = 0} = set (croots2 p)"
proof -
  from degree2_coeffs[OF assms] obtain a b c 
  where p: "p = [:c, b, a:]" and a: "a ≠ 0" by auto
  note main = roots2[OF p a]
  have 2: "2 = Suc (Suc 0)" by simp
  have coeff: "coeff p 2 = a" "coeff p 1 = b" "coeff p 0 = c" unfolding p by (auto simp: 2)
  let ?b2a = "b / (2 * a)"
  define b2a where "b2a = ?b2a"
  let ?bac = "b2a^2 - c/a"
  define bac where "bac = ?bac"
  have roots: "set (croots2 p) = {- b2a + csqrt bac, - b2a - csqrt bac}"
    unfolding croots2_def Let_def coeff b2a_def[symmetric] bac_def[symmetric]
    by (auto split: if_splits)
  show ?thesis unfolding roots main b2a_def[symmetric] bac_def[symmetric]
    using power2_eq_iff by fastforce
qed

definition rroots2 :: "real poly ⇒ real list" where
  "rroots2 p = (let a = coeff p 2; b = coeff p 1; c = coeff p 0; b2a = b / (2 * a);
    bac = b2a^2 - c/a
  in if bac = 0 then [- b2a] else if bac < 0 then []
    else let e = sqrt bac
    in
     [- b2a + e, - b2a - e])"

definition rat_roots2 :: "rat poly ⇒ rat list" where
  "rat_roots2 p = (let a = coeff p 2; b = coeff p 1; c = coeff p 0; b2a = b / (2 * a);
    bac = b2a^2 - c/a
  in map (λ e. - b2a + e) (sqrt_rat bac))"

lemma rroots2: assumes "degree p = 2"
  shows "{x. poly p x = 0} = set (rroots2 p)"
proof -
  from degree2_coeffs[OF assms] obtain a b c 
  where p: "p = [:c, b, a:]" and a: "a ≠ 0" by auto
  note main = roots2[OF p a]
  have 2: "2 = Suc (Suc 0)" by simp
  have coeff: "coeff p 2 = a" "coeff p 1 = b" "coeff p 0 = c" unfolding p by (auto simp: 2)
  let ?b2a = "b / (2 * a)"
  define b2a where "b2a = ?b2a"
  let ?bac = "b2a^2 - c/a"
  define bac where "bac = ?bac"
  have roots: "set (rroots2 p) = (if bac < 0 then {} else {- b2a + sqrt bac, - b2a - sqrt bac})"
    unfolding rroots2_def Let_def coeff b2a_def[symmetric] bac_def[symmetric]
    by (auto split: if_splits)
  show ?thesis unfolding roots main b2a_def[symmetric] bac_def[symmetric]
    by auto
qed

lemma rat_roots2: assumes "degree p = 2"
  shows "{x. poly p x = 0} = set (rat_roots2 p)"
proof -
  from degree2_coeffs[OF assms] obtain a b c 
  where p: "p = [:c, b, a:]" and a: "a ≠ 0" by auto
  note main = roots2[OF p a]
  have 2: "2 = Suc (Suc 0)" by simp
  have coeff: "coeff p 2 = a" "coeff p 1 = b" "coeff p 0 = c" unfolding p by (auto simp: 2)
  let ?b2a = "b / (2 * a)"
  define b2a where "b2a = ?b2a"
  let ?bac = "b2a^2 - c/a"
  define bac where "bac = ?bac"
  have roots: "(rat_roots2 p) = (map (λ e. -b2a + e) (sqrt_rat bac))"
    unfolding rat_roots2_def Let_def coeff b2a_def[symmetric] bac_def[symmetric] by auto
  show ?thesis unfolding roots main b2a_def[symmetric] bac_def[symmetric]
    by (auto simp: power2_eq_square)
qed


text ‹Determinining roots of complex polynomials of degree up to 2.›

definition croots :: "complex poly ⇒ complex list" where
  "croots p = (if p = 0 ∨ degree p > 2 then [] 
    else (if degree p = 0 then [] else if degree p = 1 then [roots1 p]
    else croots2 p))"

lemma croots: assumes "p ≠ 0" "degree p ≤ 2"
  shows "set (croots p) = {x. poly p x = 0}"
  using assms unfolding croots_def
  using roots0[of p] roots1[of p] croots2[of p]
  by (auto split: if_splits)

text ‹Determinining roots of real polynomials of degree up to 2.›

definition rroots :: "real poly ⇒ real list" where
  "rroots p = (if p = 0 ∨ degree p > 2 then [] 
    else (if degree p = 0 then [] else if degree p = 1 then [roots1 p]
    else rroots2 p))"

lemma rroots: assumes "p ≠ 0" "degree p ≤ 2"
  shows "set (rroots p) = {x. poly p x = 0}"
  using assms unfolding rroots_def
  using roots0[of p] roots1[of p] rroots2[of p]
  by (auto split: if_splits)

text ‹Although there is a closed form for cubic roots, 
  which is specified below, we did not yet integrate it into the
 @{const croots} and @{const rroots} algorithms.
 One obstracle is that for complex numbers, the cubic root is not
 even defined. Therefore, we also did not proof soundness of the croots3 algorithm.›

context
  fixes croot :: "nat ⇒ complex ⇒ complex"
begin

definition croots3 :: "complex poly ⇒ complex × complex × complex × complex" where
  "croots3 p = (let a = coeff p 3; b = coeff p 2; c = coeff p 1; d = coeff p 0; 
    Δ0 = b^2 - 3 * a * c;
    Δ1 = 2 * b^3 - 9 * a * b * c + 27 * a^2 * d;
    C = croot 3 ((Δ1 + csqrt ( Δ1^2 - 4 * Δ0^3)) / 2);
    u1 = 1;
    u2 = (-1 + 𝗂 * csqrt 3) / 2;
    u3 = (-1 - 𝗂 * csqrt 3) / 2;
    xk = (λ u. (-1 / (3 * a)) * (b + u * C + Δ0 / (u * C)))
    in
     (xk u1, xk u2, xk u3, a))"
end
end