Theory Cubic_Polynomials

section ‹Algorithms to compute all complex and real roots of a cubic polynomial›

theory Cubic_Polynomials
  imports 
    Cardanos_Formula
    Complex_Roots
begin

text ‹The real case where a result is only delivered if the discriminant is negative›

definition solve_depressed_cubic_Cardano_real :: "real ⇒ real ⇒ real option" where
  "solve_depressed_cubic_Cardano_real e f = (
    if e = 0 then Some (root 3 (-f)) else
     let v = - (e ^ 3 / 27) in
     case rroots2 [:v,f,1:] of 
       [u,_] ⇒ let rt = root 3 u in Some (rt - e / (3 * rt))
     | _ ⇒ None)" 

lemma solve_depressed_cubic_Cardano_real: 
  assumes "solve_depressed_cubic_Cardano_real e f = Some y" 
  shows "{y. y^3 + e * y + f = 0} = {y}"
proof (cases "e = 0")
  case True
  have "{y. y^3 + e * y + f = 0} = {y. y^3 = -f}" unfolding True 
    by (auto simp add: field_simps)
  also have "… = {root 3 (-f)}" 
    using odd_real_root_unique[of 3 _ "-f"] odd_real_root_pow[of 3] by auto
  also have "root 3 (-f) = y" using assms unfolding True solve_depressed_cubic_Cardano_real_def
    by auto
  finally show ?thesis .
next
  case False
  define v where "v = - (e ^ 3 / 27)" 
  note * = assms[unfolded solve_depressed_cubic_Cardano_real_def Let_def, folded v_def]
  let ?rr = "rroots2 [:v,f,1:]" 
  from * False obtain u u' where rr: "?rr = [u,u']" 
    by (cases ?rr; cases "tl ?rr"; cases "tl (tl ?rr)"; auto split: if_splits)
  from *[unfolded rr list.simps] False 
  have y: "y = root 3 u - e / (3 * root 3 u)" by auto
  have "u ∈ set (rroots2 [:v,f,1:])" unfolding rr by auto
  also have "set (rroots2 [:v,f,1:]) = {u. poly [:v,f,1:] u = 0}" 
    by (subst rroots2, auto)
  finally have u: "u^2 + f * u + v = 0" by (simp add: field_simps power2_eq_square)
  note Cardano = solve_cubic_depressed_Cardano_real[OF False v_def u]
  have 2: "2 = Suc (Suc 0)" by simp
  from rr have 0: "f⇧2 - 4 * v ≠ 0" unfolding rroots2_def Let_def
    by (auto split: if_splits simp: 2)
  hence 0: "discriminant_cubic_depressed e f ≠ 0" 
    unfolding discriminant_cubic_depressed_def v_def by auto
  show ?thesis using Cardano(1) Cardano(2)[OF 0] unfolding y[symmetric] by blast
qed

text ‹The complex case›

definition solve_depressed_cubic_complex :: "complex ⇒ complex ⇒ complex list" where
  "solve_depressed_cubic_complex e f = (let
          ys = (if e = 0 then all_croots 3 (- f) else (let
       u = hd (croots2 [: - (e ^ 3 / 27) ,f,1:]); 
       zs = all_croots 3 u 
       in map (λ z. z - e / (3 * z)) zs))
      in remdups ys)" 

lemma solve_depressed_cubic_complex_code[code]: 
  "solve_depressed_cubic_complex e f = (let
          ys = (if e = 0 then all_croots 3 (- f) else (let
            f2 = f / 2;
            u = - f2 + csqrt (f2^2 + e ^ 3 / 27);
            zs = all_croots 3 u 
            in map (λ z. z - e / (3 * z)) zs))
      in remdups ys)" 
  unfolding solve_depressed_cubic_complex_def Let_def croots2_def 
  by (simp add: numeral_2_eq_2)


lemma solve_depressed_cubic_complex: "y ∈ set (solve_depressed_cubic_complex e f) 
  ⟷ (y^3 + e * y + f = 0)"
proof (cases "e = 0")
  case True
  thus ?thesis by (simp add: solve_depressed_cubic_complex_def Let_def all_croots eq_neg_iff_add_eq_0)
next
  case e0: False
  hence id: "(if e = 0 then x else y) = y" for x y :: "complex list" by simp
  define v where "v = - (e ^ 3 / 27)" 
  define p where "p = [:v, f, 1:]" 
  have p2: "degree p = 2" unfolding p_def by auto
  let ?u = "hd (croots2 p)" 
  define u where "u = ?u" 
  have "u ∈ set (croots2 p)" unfolding croots2_def Let_def u_def by auto
  with croots2[OF p2] have "poly p u = 0" by auto
  hence u: "u^2 + f * u + v = 0" unfolding p_def
    by (simp add: field_simps power2_eq_square)
  note cube_roots = all_croots[of 3, simplified]
  show ?thesis unfolding solve_depressed_cubic_complex_def Let_def set_remdups set_map id cube_roots
    unfolding v_def[symmetric] p_def[symmetric] set_concat set_map
      u_def[symmetric]
  proof - 
    have p: "{x. poly p x = 0} = {u. u^2 + f * u + v = 0}" unfolding p_def by (auto simp: field_simps power2_eq_square)
    have cube: "⋃ (set ` all_croots 3 ` {x. poly p x = 0}) = {z. ∃ u. u⇧2 + f * u + v = 0 ∧ z ^ 3 = u}" 
      unfolding p by (auto simp: cube_roots)
    show "(y ∈ (λz. z - e / (3 * z)) ` {y. y ^ 3 = u}) = (y ^ 3 + e * y + f = 0)"
      using solve_cubic_depressed_Cardano_complex[OF e0 v_def u] cube by blast
  qed
qed

text ‹For the general real case, we first try Cardano with negative discrimiant and only if it is not applicable,
   then we go for the calculation using complex numbers. Note that for for non-negative delta 
   no filter is required to identify the real roots from the list of complex roots, since in that case we 
   already know that all roots are real.›
definition solve_depressed_cubic_real :: "real ⇒ real ⇒ real list" where
  "solve_depressed_cubic_real e f = (case solve_depressed_cubic_Cardano_real e f 
      of Some y ⇒ [y] 
       | None ⇒ map Re (solve_depressed_cubic_complex (of_real e) (of_real f)))"

lemma solve_depressed_cubic_real_code[code]: "solve_depressed_cubic_real e f =
  (if e = 0 then [root 3 (-f)] else 
   let v = e ^ 3 / 27; 
       f2 = f / 2;
       f2v = f2^2 + v in
   if f2v > 0 then 
     let u = -f2 + sqrt f2v;
         rt = root 3 u
      in [rt - e / (3 * rt)]
  else 
  let ce3 = of_real e / 3; 
      u = - of_real f2 + csqrt (of_real f2v) in
   map Re (remdups (map (λrt. rt - ce3 / rt) (all_croots 3 u))))" 
proof -
  have id: "rroots2 [:v, f, 1:] = (let 
     f2 = f / 2;
     bac = f2⇧2 - v in 
     if bac = 0 then [- f2] else 
     if bac < 0 then [] else let e = sqrt bac in [- f2 + e, - f2 - e])" for v
    unfolding rroots2_def Let_def numeral_2_eq_2 by auto
  define foo :: "real ⇒ real ⇒ real option" where 
    "foo f2v f2 = (case (if f2v = 0 then [- f2] else []) of [] ⇒ None | _ ⇒ None)" 
    for f2v f2
  have "solve_depressed_cubic_real e f = (if e = 0 then [root 3 (-f)] else 
   let v = e ^ 3 / 27; 
       f2 = f / 2;
       f2v = f2⇧2 + v in
   if f2v > 0 then 
     let u = -f2 + sqrt f2v;
         rt = root 3 u
      in [rt - e / (3 * rt)]
  else 
  (case foo f2v f2 of
     None ⇒ let u = - cor f2 + csqrt (cor f2v) in
   map Re
    (remdups (map (λz. z - cor e / (3 * z)) (all_croots 3 u)))
   | Some y ⇒ []))" 
    unfolding solve_depressed_cubic_real_def solve_depressed_cubic_Cardano_real_def 
      solve_depressed_cubic_complex_code
      Let_def id foo_def
    by (auto split: if_splits)
  also have id: "foo f2v f2 = None" 
    for f2v f2 unfolding foo_def by auto
  ultimately show ?thesis by (auto simp: Let_def)
qed

lemma solve_depressed_cubic_real: "y ∈ set (solve_depressed_cubic_real e f) 
  ⟷ (y^3 + e * y + f = 0)" 
proof (cases "solve_depressed_cubic_Cardano_real e f")
  case (Some x)
  show ?thesis unfolding solve_depressed_cubic_real_def Some option.simps
    using solve_depressed_cubic_Cardano_real[OF Some] by auto
next
  case None
  from this[unfolded solve_depressed_cubic_Cardano_real_def Let_def rroots2_def]
  have disc: "0 ≤ discriminant_cubic_depressed e f" unfolding discriminant_cubic_depressed_def
    by (auto split: if_splits simp: numeral_2_eq_2)
  let ?c = "complex_of_real" 
  let ?y = "?c y" 
  let ?e = "?c e" 
  let ?f = "?c f" 
  have sub: "set (solve_depressed_cubic_complex ?e ?f) ⊆ ℝ" 
  proof 
    fix y
    assume y: "y ∈ set (solve_depressed_cubic_complex ?e ?f)" 
    show "y ∈ ℝ" 
      by (rule solve_cubic_depressed_Cardano_all_real_roots[OF disc y[unfolded solve_depressed_cubic_complex]])
  qed
  have "y^3 + e * y + f = 0 ⟷ (?c (y^3 + e * y + f) = ?c 0)" unfolding of_real_eq_iff by simp
  also have "… ⟷ ?y^3 + ?e * ?y + ?f = 0" by simp
  also have "… ⟷ ?y ∈ set (solve_depressed_cubic_complex ?e ?f)" 
    unfolding solve_depressed_cubic_complex ..
  also have "… ⟷ y ∈ Re ` set (solve_depressed_cubic_complex ?e ?f)" using sub by force
  finally show ?thesis unfolding solve_depressed_cubic_real_def None by auto
qed

text ‹Combining the various algorithms›

lemma degree3_coeffs: "degree p = 3 ⟹
  ∃ a b c d. p = [: d, c, b, a :] ∧ a ≠ 0"
  by (metis One_nat_def Suc_1 degree2_coeffs degree_pCons_eq_if nat.inject numeral_3_eq_3 pCons_cases zero_neq_numeral)

definition roots3_generic :: "('a :: field_char_0 ⇒ 'a ⇒ 'a list) ⇒ 'a poly ⇒ 'a list" where
  "roots3_generic depressed_solver p = (let 
     cs = coeffs p; 
     a = cs ! 3; b = cs ! 2; c = cs ! 1; d = cs ! 0;
     a3 = 3 * a;
     ba3 = b / a3;
     b2 = b * b;
     b3 = b2 * b;
     e = (c - b2 / a3) / a;
     f = (d + 2 * b3 / (27 * a^2) - b * c / a3) / a;
     roots = depressed_solver e f
     in map (λ y. y - ba3) roots)" 

lemma roots3_generic: assumes deg: "degree p = 3" 
  and solver: "⋀ e f y. y ∈ set (depressed_solver e f) ⟷ y^3 + e * y + f = 0" 
  shows "set (roots3_generic depressed_solver p) = {x. poly p x = 0}" 
proof -
  note powers = field_simps power3_eq_cube power2_eq_square
  from degree3_coeffs[OF deg] obtain a b c d where
    p: "p = [:d,c,b,a:]" and a: "a ≠ 0" by auto
  have coeffs: "coeffs p ! 3 = a" "coeffs p ! 2 = b" "coeffs p ! 1 = c" "coeffs p ! 0 = d" 
    unfolding p using a by auto
  define e where "e = (c - b^2 / (3 * a)) / a" 
  define f where "f = (d + 2 * b^3 / (27 * a^2) - b * c / (3 * a)) / a" 
  note def = roots3_generic_def[of depressed_solver p, unfolded Let_def coeffs,
      folded power3_eq_cube, folded power2_eq_square,  folded e_def f_def]
  {
    fix x :: 'a
    define y where "y = x + b / (3 * a)" 
    have xy: "x = y - b / (3 * a)" unfolding y_def by auto
    have "poly p x = 0 ⟷ a * x^3 + b * x^2 + c * x + d = 0" unfolding p
      by (simp add: powers)
    also have "… ⟷ (y ^ 3 + e * y + f = 0)" 
      unfolding to_depressed_cubic[OF a xy e_def f_def] ..
    also have "… ⟷ y ∈ set (depressed_solver e f)" 
      unfolding solver ..
    also have "… ⟷ x ∈ set (roots3_generic depressed_solver p)" unfolding xy def by auto
    finally have "poly p x = 0 ⟷ x ∈ set (roots3_generic depressed_solver p)" by auto
  }
  thus ?thesis by auto
qed

definition croots3 :: "complex poly ⇒ complex list" where
  "croots3 = roots3_generic solve_depressed_cubic_complex"

lemma croots3: assumes deg: "degree p = 3" 
  shows "set (croots3 p) = { x. poly p x = 0}" 
  unfolding croots3_def by (rule roots3_generic[OF deg solve_depressed_cubic_complex])

definition rroots3 :: "real poly ⇒ real list" where
  "rroots3 = roots3_generic solve_depressed_cubic_real"

lemma rroots3: assumes deg: "degree p = 3" 
  shows "set (rroots3 p) = { x. poly p x = 0}" 
  unfolding rroots3_def by (rule roots3_generic[OF deg solve_depressed_cubic_real])

end