section ‹Algorithms to compute all complex and real roots of a quartic polynomial› theory Quartic_Polynomials imports Ferraris_Formula Cubic_Polynomials begin text ‹The complex case is straight-forward› definition solve_depressed_quartic_complex :: "complex ⇒ complex ⇒ complex ⇒ complex list" where "solve_depressed_quartic_complex p q r = remdups (if q = 0 then (concat (map (λ z. let y = csqrt z in [y,-y]) (croots2 [:r,p,1:]))) else let cubics = croots3 [: - (q^2), 2 * p^2 - 8 * r, 8 * p, 8:]; m = hd cubics; ― ‹select any root of the cubic polynomial› a = csqrt (2 * m); p2m = p / 2 + m; q2a = q / (2 * a); b1 = p2m - q2a; b2 = p2m + q2a in (croots2 [:b1,a,1:] @ croots2 [:b2,-a,1:]))" lemma solve_depressed_quartic_complex: "x ∈ set (solve_depressed_quartic_complex p q r) ⟷ (x^4 + p * x^2 + q * x + r = 0)" proof - note powers = field_simps power4_eq_xxxx power3_eq_cube power2_eq_square show ?thesis proof (cases "q = 0") case True have csqrt: "z = x^2 ⟷ (x = csqrt z ∨ x = - csqrt z)" for z by (metis power2_csqrt power2_eq_iff) have "(x ^ 4 + p * x⇧^{2}+ q * x + r = 0) ⟷ (x ^ 4 + p * x⇧^{2}+ r = 0)" unfolding True by simp also have "… ⟷ (∃z. z⇧^{2}+ p * z + r = 0 ∧ z = x⇧^{2})" unfolding biquadratic_solution by simp also have "… ⟷ (∃ z. poly [:r,p,1:] z = 0 ∧ z = x^2)" by (simp add: powers) also have "… ⟷ (∃ z ∈ set (croots2 [:r,p,1:]). z = x^2)" by (subst croots2[symmetric], auto) also have "… ⟷ (∃ z ∈ set (croots2 [:r,p,1:]). x = csqrt z ∨ x = - csqrt z)" unfolding csqrt .. also have "… ⟷ (x ∈ set (solve_depressed_quartic_complex p q r))" unfolding solve_depressed_quartic_complex_def id unfolding True Let_def by auto finally show ?thesis .. next case q0: False hence id: "(if q = 0 then x else y) = y" for x y :: "complex list" by auto note powers = field_simps power4_eq_xxxx power3_eq_cube power2_eq_square let ?poly = "[:- q⇧^{2}, 2 * p⇧^{2}- 8 * r, 8 * p, 8:]" from croots3[of ?poly] have croots: "set (croots3 ?poly) = {x. poly ?poly x = 0}" by auto from fundamental_theorem_of_algebra_alt[of ?poly] have "{x. poly ?poly x = 0} ≠ {}" by auto with croots have "croots3 ?poly ≠ []" by auto then obtain m rest where rts: "croots3 ?poly = m # rest" by (cases "croots3 ?poly", auto) hence hd: "hd (croots3 ?poly) = m" by auto from croots[unfolded rts] have "poly ?poly m = 0" by auto hence mrt: "8*m^3 + (8 * p) * m^2 + (2 * p^2 - 8 * r) * m - q^2 = 0" and m0: "m ≠ 0" using q0 by (auto simp: powers) define b1 where "b1 = p / 2 + m - q / (2 * csqrt (2 * m))" define b2 where "b2 = p / 2 + m + q / (2 * csqrt (2 * m))" have csqrt: "csqrt x * csqrt x = x" for x by (metis power2_csqrt power2_eq_square) show ?thesis unfolding solve_depressed_quartic_complex_def id Let_def set_remdups set_append hd unfolding b1_def[symmetric] b2_def[symmetric] apply (subst depressed_quartic_Ferrari[OF mrt q0 csqrt b1_def b2_def]) apply (subst (1 2) croots2[symmetric], auto) done qed qed text ‹The main difference in the real case is that a specific cubic root has to be used, namely a positive one. In the soundness proof we show that such a cubic root always exists.› definition solve_depressed_quartic_real :: "real ⇒ real ⇒ real ⇒ real list" where "solve_depressed_quartic_real p q r = remdups (if q = 0 then (concat (map (λ z. rroots2 [:-z,0,1:]) (rroots2 [:r,p,1:]))) else let cubics = rroots3 [: - (q^2), 2 * p^2 - 8 * r, 8 * p, 8:]; m = the (find (λ m. m > 0) cubics); ― ‹select any positive root of the cubic polynomial› a = sqrt (2 * m); p2m = p / 2 + m; q2a = q / (2 * a); b1 = p2m - q2a; b2 = p2m + q2a in (rroots2 [:b1,a,1:] @ rroots2 [:b2,-a,1:]))" lemma solve_depressed_quartic_real: "x ∈ set (solve_depressed_quartic_real p q r) ⟷ (x^4 + p * x^2 + q * x + r = 0)" proof - note powers = field_simps power4_eq_xxxx power3_eq_cube power2_eq_square show ?thesis proof (cases "q = 0") case True have sqrt: "z = x^2 ⟷ (x ∈ set (rroots2 [:-z,0,1:]))" for z by (subst rroots2[symmetric], auto simp: powers) have "(x ^ 4 + p * x⇧^{2}+ q * x + r = 0) ⟷ (x ^ 4 + p * x⇧^{2}+ r = 0)" unfolding True by simp also have "… ⟷ (∃z. z⇧^{2}+ p * z + r = 0 ∧ z = x⇧^{2})" unfolding biquadratic_solution by simp also have "… ⟷ (∃ z. poly [:r,p,1:] z = 0 ∧ z = x^2)" by (simp add: powers) also have "… ⟷ (∃ z ∈ set (rroots2 [:r,p,1:]). z = x^2)" by (subst rroots2[symmetric], auto) also have "… ⟷ (∃ z ∈ set (rroots2 [:r,p,1:]). x ∈ set (rroots2 [:-z,0,1:]))" unfolding sqrt .. also have "… ⟷ (x ∈ set (solve_depressed_quartic_real p q r))" unfolding solve_depressed_quartic_real_def id unfolding True Let_def by auto finally show ?thesis .. next case q0: False hence id: "(if q = 0 then x else y) = y" for x y :: "real list" by auto note powers = field_simps power4_eq_xxxx power3_eq_cube power2_eq_square let ?poly = "[:- q⇧^{2}, 2 * p⇧^{2}- 8 * r, 8 * p, 8:]" define cubics where "cubics = rroots3 ?poly" from rroots3[of ?poly, folded cubics_def] have rroots: "set cubics = {x. poly ?poly x = 0}" by auto from odd_degree_imp_real_root[of ?poly] have "{x. poly ?poly x = 0} ≠ {}" by auto with rroots have "cubics ≠ []" by auto have "∃ m. m ∈ set cubics ∧ m > 0" proof (rule ccontr) assume "¬ ?thesis" from this[unfolded rroots] have rt: "poly ?poly m = 0 ⟹ m ≤ 0" for m by auto have "poly ?poly 0 = - (q^2)" by simp also have "… < 0" using q0 by auto finally have lt: "poly ?poly 0 ≤ 0" by simp from poly_pinfty_gt_lc[of ?poly] obtain n0 where "⋀ n. n ≥ n0 ⟹ 8 ≤ poly ?poly n" by auto from this[of "max n0 0"] have "poly ?poly (max n0 0) ≥ 0" by auto from IVT[of "poly ?poly", OF lt this] obtain m where "m ≥ 0" and poly: "poly ?poly m = 0" by auto from rt[OF this(2)] this(1) have "m = 0" by auto thus False using poly q0 by simp qed hence "find (λ m. m > 0) cubics ≠ None" unfolding find_None_iff by auto then obtain m where find: "find (λ m. m > 0) cubics = Some m" by auto from find_Some_D[OF this] have m: "m ∈ set cubics" and m_0: "m > 0" by auto with rroots have "poly ?poly m = 0" by auto hence mrt: "8*m^3 + (8 * p) * m^2 + (2 * p^2 - 8 * r) * m - q^2 = 0" by (auto simp: powers) from m_0 have sqrt: "sqrt (2 * m) * sqrt (2 * m) = 2 * m" by simp define b1 where "b1 = p / 2 + m - q / (2 * sqrt (2 * m))" define b2 where "b2 = p / 2 + m + q / (2 * sqrt (2 * m))" show ?thesis unfolding solve_depressed_quartic_real_def id Let_def set_remdups set_append cubics_def[symmetric] find option.sel unfolding b1_def[symmetric] b2_def[symmetric] apply (subst depressed_quartic_Ferrari[OF mrt q0 sqrt b1_def b2_def]) apply (subst (1 2) rroots2[symmetric], auto) done qed qed text ‹Combining the various algorithms› lemma numeral_4_eq_4: "4 = Suc (Suc (Suc (Suc 0)))" by (simp add: eval_nat_numeral) lemma degree4_coeffs: "degree p = 4 ⟹ ∃ a b c d e. p = [: e, d, c, b, a :] ∧ a ≠ 0" using degree3_coeffs degree_pCons_eq_if nat.inject numeral_3_eq_3 numeral_4_eq_4 pCons_cases zero_neq_numeral by metis definition roots4_generic :: "('a :: field_char_0 ⇒ 'a ⇒ 'a ⇒ 'a list) ⇒ 'a poly ⇒ 'a list" where "roots4_generic depressed_solver p = (let cs = coeffs p; cs = coeffs p; a4 = cs ! 4; a3 = cs ! 3; a2 = cs ! 2; a1 = cs ! 1; a0 = cs ! 0; b = a3 / a4; c = a2 / a4; d = a1 / a4; e = a0 / a4; b2 = b * b; b3 = b2 * b; b4 = b3 * b; b4' = b / 4; p = c - 3/8 * b2; q = (b3 - 4*b*c + 8 * d) / 8; r = ( -3 * b4 + 256 * e - 64 * b * d + 16 * b2 * c) / 256; roots = depressed_solver p q r in map (λ y. y - b4') roots)" lemma roots4_generic: assumes deg: "degree p = 4" and solver: "⋀ p q r y. y ∈ set (depressed_solver p q r) ⟷ y^4 + p * y^2 + q * y + r = 0" shows "set (roots4_generic depressed_solver p) = {x. poly p x = 0}" proof - note powers = field_simps power4_eq_xxxx power3_eq_cube power2_eq_square from degree4_coeffs[OF deg] obtain a4 a3 a2 a1 a0 where p: "p = [:a0,a1,a2,a3,a4:]" and a4: "a4 ≠ 0" by auto have coeffs: "coeffs p ! 4 = a4" "coeffs p ! 3 = a3" "coeffs p ! 2 = a2" "coeffs p ! 1 = a1" "coeffs p ! 0 = a0" unfolding p using a4 by auto define b where "b = a3 / a4" define c where "c = a2 / a4" define d where "d = a1 / a4" define e where "e = a0 / a4" note def = roots4_generic_def[of depressed_solver p, unfolded Let_def coeffs, folded b_def c_def d_def e_def, folded power4_eq_xxxx, folded power3_eq_cube, folded power2_eq_square] let ?p = p { fix x define y where "y = x + b / 4" define p where "p = c - (3/8) * b^2" define q where "q = (b^3 - 4*b*c + 8 * d) / 8" define r where "r = ( -3 * b^4 + 256 * e - 64 * b * d + 16 * b^2 * c) / 256" note def = def[folded p_def q_def r_def] have xy: "x = y - b / 4" unfolding y_def by auto have "poly ?p x = 0 ⟷ a4 * x^4 + a3 * x^3 + a2 * x^2 + a1 * x + a0 = 0" unfolding p by (simp add: powers) also have "… ⟷ (y ^ 4 + p * y⇧^{2}+ q * y + r = 0)" unfolding to_depressed_quartic[OF a4 b_def c_def d_def e_def p_def q_def r_def xy] .. also have "… ⟷ y ∈ set (depressed_solver p q r)" unfolding solver .. also have "… ⟷ x ∈ set (roots4_generic depressed_solver ?p)" unfolding xy def by auto finally have "poly ?p x = 0 ⟷ x ∈ set (roots4_generic depressed_solver ?p)" by auto } thus ?thesis by simp qed definition croots4 :: "complex poly ⇒ complex list" where "croots4 = roots4_generic solve_depressed_quartic_complex" lemma croots4: assumes deg: "degree p = 4" shows "set (croots4 p) = { x. poly p x = 0}" unfolding croots4_def by (rule roots4_generic[OF deg solve_depressed_quartic_complex]) definition rroots4 :: "real poly ⇒ real list" where "rroots4 = roots4_generic solve_depressed_quartic_real" lemma rroots4: assumes deg: "degree p = 4" shows "set (rroots4 p) = { x. poly p x = 0}" unfolding rroots4_def by (rule roots4_generic[OF deg solve_depressed_quartic_real]) end