section ‹Cardano's formula for solving cubic equations› theory Cardanos_Formula imports Polynomial_Factorization.Explicit_Roots Polynomial_Interpolation.Ring_Hom_Poly Complex_Geometry.More_Complex Algebraic_Numbers.Complex_Roots_Real_Poly begin subsection ‹Translation to depressed case› text ‹Solving an arbitrary cubic equation can easily be turned into the depressed case, i.e., where there is no quadratic part.› lemma to_depressed_cubic: fixes a :: "'a :: field_char_0" assumes a: "a ≠ 0" and xy: "x = y - b / (3 * a)" and e: "e = (c - b^2 / (3 * a)) / a" and f: "f = (d + 2 * b^3 / (27 * a^2) - b * c / (3 * a)) / a" shows "(a * x ^ 3 + b * x⇧^{2}+ c * x + d = 0) ⟷ y^3 + e * y + f = 0" proof - let ?yexp = "y^3 + e * y + f" have "a * x^3 + b * x^2 + c * x + d = 0 ⟷ (a * x^3 + b * x^2 + c * x + d) / a = 0" using a by auto also have "(a * x^3 + b * x^2 + c * x + d) / a = ?yexp" unfolding xy e f power3_eq_cube power2_eq_square using a by (simp add: field_simps) finally show ?thesis . qed subsection ‹Solving the depressed case in arbitrary fields› lemma cubic_depressed: fixes e :: "'a :: field_char_0" assumes yz: "e ≠ 0 ⟹ z^2 - y * z - e / 3 = 0" and u: "e ≠ 0 ⟹ u = z^3" and v: "v = - (e ^ 3 / 27)" shows "y^3 + e * y + f = 0 ⟷ (if e = 0 then y^3 = -f else u⇧^{2}+ f * u + v = 0)" proof - let ?yexp = "y^3 + e * y + f" show ?thesis proof (cases "e = 0") case False note yz = yz[OF False] from yz have eyz: "e = 3 * (z^2 - y * z)" by auto from yz False have z0: "z ≠ 0" by auto have "?yexp = 0 ⟷ z^3 * ?yexp = 0" using z0 by simp also have "z^3 * ?yexp = z^6 + f * z^3 - e^3/27" unfolding eyz by algebra also have "… = u^2 + f * u + v" unfolding u[OF False] v by algebra finally show ?thesis using False by auto next case True show ?thesis unfolding True by (auto, algebra) qed qed subsection ‹Solving the depressed case for complex numbers› text ‹In the complex-numbers-case, the quadratic equation for u is always solvable, and the main challenge here is prove that it does not matter which solution of the quadratic equation is considered (this is the diff:False case in the proof below.)› lemma solve_cubic_depressed_Cardano_complex: fixes e :: complex assumes e0: "e ≠ 0" and v: "v = - (e ^ 3 / 27)" and u: "u^2 + f * u + v = 0" shows "y^3 + e * y + f = 0 ⟷ (∃ z. z^3 = u ∧ y = z - e / (3 * z))" proof - from v e0 have v0: "v ≠ 0" by auto from e0 have "(if e = 0 then x else y) = y" for x y :: bool by auto note main = cubic_depressed[OF _ _ v, unfolded this] show ?thesis (is "?l = ?r") proof assume ?r then obtain z where z: "z^3 = u" and y: "y = z - e / (3 * z)" by auto from u v0 have u0: "u ≠ 0" by auto from z u0 have z0: "z ≠ 0" by auto show ?l proof (subst main) show "u⇧^{2}+ f * u + v = 0" by fact show "u = z^3" unfolding z by simp show "z⇧^{2}- y * z - e / 3 = 0" unfolding y using z0 by (auto simp: field_simps power2_eq_square) qed next assume ?l let ?yexp = "y^3 + e * y + f" have y0: "?yexp = 0" using ‹?l› by auto define p where "p = [: -e/3, -y, 1:]" have deg: "degree p = 2" unfolding p_def by auto define z where "z = hd (croots2 p)" have "z ∈ set (croots2 p)" unfolding croots2_def Let_def z_def by auto with croots2[OF deg] have pz: "poly p z = 0" by auto from pz e0 have z0: "z ≠ 0" unfolding p_def by auto from pz have yz: "y * z = z * z - e / 3" unfolding p_def by (auto simp: field_simps) from arg_cong[OF this, of "λ x. x / z"] z0 have "y = z - e / (3 * z)" by (auto simp: field_simps) have "∃ u z. u⇧^{2}+ f * u + v = 0 ∧ z^3 = u ∧ y = z - e / (3 * z)" proof (intro exI conjI) show "y = z - e / (3 * z)" by fact from y0 have "0 = ?yexp * z^3" by auto also have "… = (y * z)^3 + e * (y * z) * z^2 + f * z^3" by algebra also have "… = (z^3)^2 + f * (z^3) + v" unfolding yz v by algebra finally show "(z^3)^2 + f * (z^3) + v = 0" by simp qed simp then obtain uu z where *: "uu⇧^{2}+ f * uu + v = 0" "z ^ 3 = uu" "y = z - e / (3 * z)" by blast show ?r proof (cases "uu = u") case True thus ?thesis using * by auto next case diff: False define p where "p = [:v,f,1:]" have p2: "degree p = 2" unfolding p_def by auto have poly: "poly p u = 0" "poly p uu = 0" using u *(1) unfolding p_def by (auto simp: field_simps power2_eq_square) have u0: "u ≠ 0" "uu ≠ 0" using poly v0 unfolding p_def by auto { from poly(1) have "[:-u,1:] dvd p" by (meson poly_eq_0_iff_dvd) then obtain q where pq: "p = q * [:-u,1:]" by auto from poly(2)[unfolded pq poly_mult] diff have "poly q uu = 0" by auto hence "[:-uu,1:] dvd q" by (meson poly_eq_0_iff_dvd) then obtain q' where qq': "q = q' * [:-uu,1:]" by auto with pq have pq: "p = q' * [:-uu,1:] * [:-u,1:]" by auto from pq[unfolded p_def] have q': "q' ≠ 0" by auto from arg_cong[OF pq, of degree, unfolded p2] have "2 = degree (q' * [:- uu, 1:] * [:- u, 1:])" . also have "… = degree q' + degree [:- uu, 1:] + degree [:- u, 1:]" apply (subst degree_mult_eq) subgoal using q' by (metis mult_eq_0_iff pCons_eq_0_iff zero_neq_one) subgoal by force by (subst degree_mult_eq[OF q'], auto) also have "… = degree q' + 2" by simp finally have dq: "degree q' = 0" by simp from dq obtain c where q': "q' = [: c:]" by (metis degree_eq_zeroE) from pq[unfolded q' p_def] have "c = 1" by auto with q' have "q' = 1" by simp with pq have "[: -u, 1:] * [: -uu, 1 :] = p" by simp } from this[unfolded p_def, simplified] have prod: "uu * u = v" by simp hence uu: "u = v / uu" using u0 by (simp add: field_simps) define zz where "zz = - e / (3 * z)" show ?r using *(2-) uu unfolding v using u0 by (intro exI[of _ zz], auto simp: zz_def field_simps) qed qed qed subsection ‹Solving the depressed case for real numbers› definition discriminant_cubic_depressed :: "'a :: comm_ring_1 ⇒ 'a ⇒ 'a" where "discriminant_cubic_depressed e f = - (4 * e^3 + 27 * f^2)" lemma discriminant_cubic_depressed: assumes "[:-x,1:] * [:-y,1:] * [:-z,1:] = [:f,e,0,1:]" shows "discriminant_cubic_depressed e f = (x-y)^2 * (x - z)^2 * (y - z)^2" proof - from assms have f: "f = - (z * (y * x))" and e: "e = y * x - z * (- y - x)" and z: "z = - y - x" by auto show ?thesis unfolding discriminant_cubic_depressed_def e f z by (simp add: power2_eq_square power3_eq_cube field_simps) qed text ‹If the discriminant is negative, then there is exactly one real root› lemma solve_cubic_depressed_Cardano_real: fixes e f v u :: real defines "y1 ≡ root 3 u - e / (3 * root 3 u)" and "Δ ≡ discriminant_cubic_depressed e f" assumes e0: "e ≠ 0" and v: "v = - (e ^ 3 / 27)" and u: "u⇧^{2}+ f * u + v = 0" (* this implies Δ ≤ 0 *) shows "y1^3 + e * y1 + f = 0" "Δ ≠ 0 ⟹ y^3 + e * y + f = 0 ⟹ y = y1" (* this is the case Δ < 0 *) proof - let ?c = complex_of_real let ?y = "?c y" let ?e = "?c e" let ?u = "?c u" let ?v = "?c v" let ?f = "?c f" { fix y :: real let ?y = "?c y" have "y^3 + e * y + f = 0 ⟷ ?c (y^3 + e * y + f) = ?c 0" using of_real_eq_iff by blast also have "… ⟷ ?y^3 + ?e * ?y + ?f = 0" by simp also have "… ⟷ (∃ z. z^3 = ?u ∧ ?y = z - ?e / (3 * z))" proof (rule solve_cubic_depressed_Cardano_complex) show "?e ≠ 0" using e0 by auto show "?v = - (?e ^ 3 / 27)" unfolding v by simp show "?u⇧^{2}+ ?f * ?u + ?v = 0" using arg_cong[OF u, of ?c] by simp qed finally have "y^3 + e * y + f = 0 ⟷ (∃ z. z^3 = ?u ∧ ?y = z - ?e / (3 * z))" . } note pre = this show y1: "y1^3 + e * y1 + f = 0" unfolding pre y1_def by (intro exI[of _ "?c (root 3 u)"], simp only: of_real_power[symmetric], simp del: of_real_power add: odd_real_root_pow) from u have "{z. poly [:v,f,1:] z = 0} ≠ {}" by (auto simp add: field_simps power2_eq_square) hence "set (rroots2 [:v,f,1:]) ≠ {}" by (subst rroots2[symmetric], auto) hence "rroots2 [:v,f,1:] ≠ []" by simp from this[unfolded rroots2_def Let_def, simplified] have "f^2 - 4 * v ≥ 0" by (auto split: if_splits simp: numeral_2_eq_2 field_simps power2_eq_square) hence delta_le_0: "Δ ≤ 0" unfolding Δ_def discriminant_cubic_depressed_def v by auto assume Delta_non_0: "Δ ≠ 0" with delta_le_0 have delta_neg: "Δ < 0" by simp let ?p = "[:f,e,0,1:]" have poly: "poly ?p y = 0 ⟷ y^3 + e * y + f = 0" for y by (simp add: field_simps power2_eq_square power3_eq_cube) from y1 have "poly ?p y1 = 0" unfolding poly . hence "[:-y1,1:] dvd ?p" using poly_eq_0_iff_dvd by blast then obtain q where pq: "?p = [:-y1,1:] * q" by blast { fix y2 assume "poly ?p y2 = 0" "y2 ≠ y1" from this[unfolded pq] poly_mult have "poly q y2 = 0" by auto from this[unfolded poly_eq_0_iff_dvd] obtain r where qr: "q = [:-y2,1:] * r" by blast { have r0: "r ≠ 0" using pq unfolding qr poly_mult by auto have "3 = degree ?p" by simp also have "… = 2 + degree r" unfolding pq qr apply (subst degree_mult_eq, force) subgoal using r0 pq qr by force by (subst degree_mult_eq[OF _ r0], auto) finally have "degree r = 1" by simp from degree1_coeffs[OF this] obtain yy a where r: "r = [:yy,a:]" by auto define y3 where "y3 = -yy" with r have r: "r = [:-y3,a:]" by auto from pq[unfolded qr r] have "a = 1" by auto with r have "∃ y3. r = [:-y3,1:]" by auto } then obtain y3 where r: "r = [:-y3,1:]" by auto have py: "?p = [:-y1,1:] * [:-y2,1:] * [:-y3,1:]" unfolding pq qr r by algebra from discriminant_cubic_depressed[OF this[symmetric], folded Δ_def] have delta: "Δ = (y1 - y2)⇧^{2}* (y1 - y3)⇧^{2}* (y2 - y3)⇧^{2}" . have d0: "Δ ≥ 0" unfolding delta by auto with delta_neg have False by auto } with y1 show "y^3 + e * y + f = 0 ⟹ y = y1" unfolding poly by auto qed text ‹If the discriminant is non-negative, then all roots are real› lemma solve_cubic_depressed_Cardano_all_real_roots: fixes e f v :: real and y :: complex defines "Δ ≡ discriminant_cubic_depressed e f" assumes Delta: "Δ ≥ 0" and rt: "y^3 + e * y + f = 0" shows "y ∈ ℝ" proof - note powers = field_simps power3_eq_cube power2_eq_square let ?c = complex_of_real let ?e = "?c e" let ?f = "?c f" let ?cp = "[:?f,?e,0,1:]" let ?p = "[:f,e,0,1:]" from odd_degree_imp_real_root[of ?p] obtain x1 where "poly ?p x1 = 0" by auto hence "[:-x1,1:] dvd ?p" using poly_eq_0_iff_dvd by blast then obtain q where pq: "?p = [:-x1,1:] * q" by auto from arg_cong[OF pq, of degree] have "3 = degree ([:-x1,1:] * q)" by simp also have "… = 1 + degree q" by (subst degree_mult_eq, insert pq, auto) finally have dq: "degree q = 2" by auto let ?cc = "map_poly ?c" let ?q = "?cc q" have cpq: "?cc ?p = ?cc [:-x1,1:] * ?q" unfolding pq hom_distribs by simp let ?x1 = "?c x1" have dq': "degree ?q = 2" using dq by simp have "¬ constant (poly ?q)" using dq by (simp add: constant_degree) from fundamental_theorem_of_algebra[OF this] obtain x2 where x2: "poly ?q x2 = 0" by blast have "x2 ∈ ℝ" proof (rule ccontr) assume x2r: "x2 ∉ ℝ" define x3 where "x3 = cnj x2" from x2r have x23: "x2 ≠ x3" unfolding x3_def using Reals_cnj_iff by force have x3: "poly ?q x3 = 0" unfolding x3_def by (rule complex_conjugate_root[OF _ x2], auto) from x2[unfolded poly_eq_0_iff_dvd] obtain r where qr: "?q = [:-x2,1:] * r" by auto from arg_cong[OF this[symmetric], of "λ x. poly x x3", unfolded poly_mult x3 mult_eq_0_iff] x23 have x3: "poly r x3 = 0" by auto from arg_cong[OF qr, of degree] have "2 = degree ([:-x2,1:] * r)" using dq' by simp also have "… = 1 + degree r" by (subst degree_mult_eq, insert pq qr, auto) finally have "degree r = 1" by simp from degree1_coeffs[OF this] obtain a b where r: "r = [:a,b:]" by auto from cpq[unfolded qr r] have b1: "b = 1" by simp with x3 r have "a + x3 = 0" by simp hence "a = - x3" by algebra with b1 r have r: "r = [:-x3,1:]" by auto have "?cc ?p = ?cc [:-x1,1:] * [:-x2,1:] * [:-x3,1:]" unfolding cpq qr r by algebra also have "?cc [:-x1,1:] = [:-?x1,1:]" by simp also have "?cc ?p = ?cp" by simp finally have id: "[:-?x1,1:] * [:-x2,1:] * [:-x3,1:] = ?cp" by simp define x23 where "x23 = - 4 * (Im x2)^2" define x12c where "x12c = ?x1 - x2" define x12 where "x12 = (Re x12c) ^ 2 + (Im x12c)^2" have x23_0: "x23 < 0" unfolding x23_def using x2r using complex_is_Real_iff by force have "Im x12c ≠ 0" unfolding x12c_def using x2r using complex_is_Real_iff by force hence "(Im x12c)^2 > 0" by simp hence x12: "x12 > 0" unfolding x12_def using sum_power2_gt_zero_iff by auto from discriminant_cubic_depressed[OF id] have "?c Δ = ((?x1 - x2)⇧^{2}* (?x1 - x3)⇧^{2}) * (x2 - x3)⇧^{2}" unfolding Δ_def discriminant_cubic_depressed_def by simp also have "(x2 - x3)^2 = ?c x23" unfolding x3_def x23_def by (simp add: complex_eq_iff power2_eq_square) also have "(?x1 - x2)⇧^{2}* (?x1 - x3)⇧^{2}= ((?x1 - x2) * (?x1 - x3))^2" by (simp add: power2_eq_square) also have "?x1 - x3 = cnj (?x1 - x2)" unfolding x3_def by simp also have "(?x1 - x2) = x12c" unfolding x12c_def .. also have "x12c * cnj x12c = ?c x12" by (simp add: x12_def complex_eq_iff power2_eq_square) finally have "?c Δ = ?c (x12^2 * x23)" by simp hence "Δ = x12^2 * x23" by (rule of_real_hom.injectivity) also have "… < 0" using x12 x23_0 by (meson mult_pos_neg zero_less_power) finally show False using Delta by simp qed with x2 obtain x2 where "poly ?q (?c x2) = 0" unfolding Reals_def by auto hence x2: "poly q x2 = 0" by simp from x2[unfolded poly_eq_0_iff_dvd] obtain r where qr: "q = [:-x2,1:] * r" by auto from arg_cong[OF qr, of degree] have "2 = degree ([:-x2,1:] * r)" using dq' by simp also have "… = 1 + degree r" by (subst degree_mult_eq, insert pq qr, auto) finally have "degree r = 1" by simp from degree1_coeffs[OF this] obtain a b where r: "r = [:a,b:]" by auto from pq[unfolded qr r] have b1: "b = 1" by simp define x3 where "x3 = -a" have r: "r = [:-x3,1:]" unfolding r b1 x3_def by simp let ?pp = "[:-x1,1:] * [:-x2,1:] * [:-x3,1:]" have id: "?p = ?pp" unfolding pq qr r by linarith have "True ⟷ y^3 + e * y + f = 0" using rt by auto also have "y^3 + e * y + f = poly (?cc ?p) y" by (simp add: powers) also have "… = poly (?cc ?pp) y" unfolding id by simp also have "?cc ?pp = [:-?c x1, 1:] * [:-?c x2, 1:] * [:- ?c x3, 1:]" by simp also have "poly … y = 0 ⟷ y = ?c x1 ∨ y = ?c x2 ∨ y = ?c x3" unfolding poly_mult mult_eq_0_iff by auto finally show "y ∈ ℝ" by auto qed end