src/HOL/ex/Cubic_Quartic.thy
author wenzelm
Wed Jun 22 10:09:20 2016 +0200 (2016-06-22)
changeset 63343 fb5d8a50c641
parent 63054 1b237d147cc4
child 63589 58aab4745e85
permissions -rw-r--r--
bundle lifting_syntax;
kleing@59190
     1
(*  Title:      HOL/ex/Cubic_Quartic.thy
kleing@59190
     2
    Author:     Amine Chaieb
kleing@59190
     3
*)
kleing@59190
     4
wenzelm@63054
     5
section \<open>The Cubic and Quartic Root Formulas\<close>
kleing@59190
     6
kleing@59190
     7
theory Cubic_Quartic
kleing@59190
     8
imports Complex_Main
kleing@59190
     9
begin
kleing@59190
    10
wenzelm@63054
    11
section \<open>The Cubic Formula\<close>
kleing@59190
    12
kleing@59190
    13
definition "ccbrt z = (SOME (w::complex). w^3 = z)"
kleing@59190
    14
kleing@59190
    15
lemma ccbrt: "(ccbrt z) ^ 3 = z"
wenzelm@63054
    16
proof -
wenzelm@63054
    17
  from rcis_Ex obtain r a where ra: "z = rcis r a"
wenzelm@63054
    18
    by blast
kleing@59190
    19
  let ?r' = "if r < 0 then - root 3 (-r) else root 3 r"
kleing@59190
    20
  let ?a' = "a/3"
wenzelm@63054
    21
  have "rcis ?r' ?a' ^ 3 = rcis r a"
wenzelm@63054
    22
    by (cases "r < 0") (simp_all add: DeMoivre2)
wenzelm@63054
    23
  then have *: "\<exists>w. w^3 = z"
wenzelm@63054
    24
    unfolding ra by blast
wenzelm@63054
    25
  from someI_ex [OF *] show ?thesis
wenzelm@63054
    26
    unfolding ccbrt_def by blast
kleing@59190
    27
qed
kleing@59190
    28
wenzelm@63054
    29
wenzelm@63054
    30
text \<open>The reduction to a simpler form:\<close>
kleing@59190
    31
kleing@59190
    32
lemma cubic_reduction:
kleing@59190
    33
  fixes a :: complex
wenzelm@63054
    34
  assumes
wenzelm@63054
    35
    "a \<noteq> 0 \<and> x = y - b / (3 * a) \<and>  p = (3* a * c - b^2) / (9 * a^2) \<and>
wenzelm@63054
    36
      q = (9 * a * b * c - 2 * b^3 - 27 * a^2 * d) / (54 * a^3)"
kleing@59190
    37
  shows "a * x^3 + b * x^2 + c * x + d = 0 \<longleftrightarrow> y^3 + 3 * p * y - 2 * q = 0"
wenzelm@63054
    38
proof -
wenzelm@63054
    39
  from assms have "3 * a \<noteq> 0" "9 * a^2 \<noteq> 0" "54 * a^3 \<noteq> 0" by auto
wenzelm@63054
    40
  then have *:
wenzelm@63054
    41
      "x = y - b / (3 * a) \<longleftrightarrow> (3*a) * x = (3*a) * y - b"
wenzelm@63054
    42
      "p = (3* a * c - b^2) / (9 * a^2) \<longleftrightarrow> (9 * a^2) * p = (3* a * c - b^2)"
wenzelm@63054
    43
      "q = (9 * a * b * c - 2 * b^3 - 27 * a^2 * d) / (54 * a^3) \<longleftrightarrow>
wenzelm@63054
    44
        (54 * a^3) * q = (9 * a * b * c - 2 * b^3 - 27 * a^2 * d)"
kleing@59190
    45
    by (simp_all add: field_simps)
wenzelm@63054
    46
  from assms [unfolded *] show ?thesis
wenzelm@63054
    47
    by algebra
kleing@59190
    48
qed
kleing@59190
    49
kleing@59190
    50
wenzelm@63054
    51
text \<open>The solutions of the special form:\<close>
wenzelm@63054
    52
wenzelm@63054
    53
lemma cubic_basic:
kleing@59190
    54
  fixes s :: complex
wenzelm@63054
    55
  assumes
wenzelm@63054
    56
    "s^2 = q^2 + p^3 \<and>
wenzelm@63054
    57
      s1^3 = (if p = 0 then 2 * q else q + s) \<and>
wenzelm@63054
    58
      s2 = -s1 * (1 + i * t) / 2 \<and>
wenzelm@63054
    59
      s3 = -s1 * (1 - i * t) / 2 \<and>
wenzelm@63054
    60
      i^2 + 1 = 0 \<and>
wenzelm@63054
    61
      t^2 = 3"
wenzelm@63054
    62
  shows
kleing@59190
    63
    "if p = 0
kleing@59190
    64
     then y^3 + 3 * p * y - 2 * q = 0 \<longleftrightarrow> y = s1 \<or> y = s2 \<or> y = s3
kleing@59190
    65
     else s1 \<noteq> 0 \<and>
wenzelm@63054
    66
      (y^3 + 3 * p * y - 2 * q = 0 \<longleftrightarrow> (y = s1 - p / s1 \<or> y = s2 - p / s2 \<or> y = s3 - p / s3))"
wenzelm@63054
    67
proof (cases "p = 0")
wenzelm@63054
    68
  case True
wenzelm@63054
    69
  with assms show ?thesis
wenzelm@63054
    70
    by (simp add: field_simps) algebra
wenzelm@63054
    71
next
wenzelm@63054
    72
  case False
wenzelm@63054
    73
  with assms have *: "s1 \<noteq> 0" by (simp add: field_simps) algebra
wenzelm@63054
    74
  with assms False have "s2 \<noteq> 0" "s3 \<noteq> 0"
wenzelm@63054
    75
    by (simp_all add: field_simps) algebra+
wenzelm@63054
    76
  with * have **:
wenzelm@63054
    77
      "y = s1 - p / s1 \<longleftrightarrow> s1 * y = s1^2 - p"
wenzelm@63054
    78
      "y = s2 - p / s2 \<longleftrightarrow> s2 * y = s2^2 - p"
wenzelm@63054
    79
      "y = s3 - p / s3 \<longleftrightarrow> s3 * y = s3^2 - p"
wenzelm@63054
    80
    by (simp_all add: field_simps power2_eq_square)
wenzelm@63054
    81
  from assms False show ?thesis
wenzelm@63054
    82
    unfolding ** by (simp add: field_simps) algebra
kleing@59190
    83
qed
kleing@59190
    84
wenzelm@63054
    85
wenzelm@63054
    86
text \<open>Explicit formula for the roots:\<close>
kleing@59190
    87
kleing@59190
    88
lemma cubic:
kleing@59190
    89
  assumes a0: "a \<noteq> 0"
wenzelm@63054
    90
  shows
wenzelm@63054
    91
    "let
wenzelm@63054
    92
      p = (3 * a * c - b^2) / (9 * a^2);
kleing@59190
    93
      q = (9 * a * b * c - 2 * b^3 - 27 * a^2 * d) / (54 * a^3);
kleing@59190
    94
      s = csqrt(q^2 + p^3);
kleing@59190
    95
      s1 = (if p = 0 then ccbrt(2 * q) else ccbrt(q + s));
kleing@59190
    96
      s2 = -s1 * (1 + ii * csqrt 3) / 2;
kleing@59190
    97
      s3 = -s1 * (1 - ii * csqrt 3) / 2
wenzelm@63054
    98
     in
wenzelm@63054
    99
      if p = 0 then
wenzelm@63054
   100
        a * x^3 + b * x^2 + c * x + d = 0 \<longleftrightarrow>
wenzelm@63054
   101
          x = s1 - b / (3 * a) \<or>
wenzelm@63054
   102
          x = s2 - b / (3 * a) \<or>
wenzelm@63054
   103
          x = s3 - b / (3 * a)
kleing@59190
   104
      else
kleing@59190
   105
        s1 \<noteq> 0 \<and>
wenzelm@63054
   106
          (a * x^3 + b * x^2 + c * x + d = 0 \<longleftrightarrow>
kleing@59190
   107
            x = s1 - p / s1 - b / (3 * a) \<or>
kleing@59190
   108
            x = s2 - p / s2 - b / (3 * a) \<or>
kleing@59190
   109
            x = s3 - p / s3 - b / (3 * a))"
wenzelm@63054
   110
proof -
kleing@59190
   111
  let ?p = "(3 * a * c - b^2) / (9 * a^2)"
kleing@59190
   112
  let ?q = "(9 * a * b * c - 2 * b^3 - 27 * a^2 * d) / (54 * a^3)"
wenzelm@63054
   113
  let ?s = "csqrt (?q^2 + ?p^3)"
kleing@59190
   114
  let ?s1 = "if ?p = 0 then ccbrt(2 * ?q) else ccbrt(?q + ?s)"
kleing@59190
   115
  let ?s2 = "- ?s1 * (1 + ii * csqrt 3) / 2"
kleing@59190
   116
  let ?s3 = "- ?s1 * (1 - ii * csqrt 3) / 2"
kleing@59190
   117
  let ?y = "x + b / (3 * a)"
wenzelm@63054
   118
  from a0 have zero: "9 * a^2 \<noteq> 0" "a^3 * 54 \<noteq> 0" "(a * 3) \<noteq> 0"
wenzelm@63054
   119
    by auto
wenzelm@63054
   120
  have eq: "a * x^3 + b * x^2 + c * x + d = 0 \<longleftrightarrow> ?y^3 + 3 * ?p * ?y - 2 * ?q = 0"
kleing@59190
   121
    by (rule cubic_reduction) (auto simp add: field_simps zero a0)
kleing@59190
   122
  have "csqrt 3^2 = 3" by (rule power2_csqrt)
wenzelm@63054
   123
  then have th0:
wenzelm@63054
   124
    "?s^2 = ?q^2 + ?p ^ 3 \<and> ?s1^ 3 = (if ?p = 0 then 2 * ?q else ?q + ?s) \<and>
wenzelm@63054
   125
      ?s2 = - ?s1 * (1 + ii * csqrt 3) / 2 \<and>
wenzelm@63054
   126
      ?s3 = - ?s1 * (1 - ii * csqrt 3) / 2 \<and>
wenzelm@63054
   127
      ii^2 + 1 = 0 \<and> csqrt 3^2 = 3"
nipkow@62361
   128
    using zero by (simp add: field_simps ccbrt)
kleing@59190
   129
  from cubic_basic[OF th0, of ?y]
wenzelm@63054
   130
  show ?thesis
kleing@59190
   131
    apply (simp only: Let_def eq)
nipkow@62361
   132
    using zero apply (simp add: field_simps ccbrt)
kleing@59190
   133
    using zero
wenzelm@63054
   134
    apply (cases "a * (c * 3) = b^2")
wenzelm@63054
   135
    apply (simp_all add: field_simps)
kleing@59190
   136
    done
kleing@59190
   137
qed
kleing@59190
   138
kleing@59190
   139
wenzelm@63054
   140
section \<open>The Quartic Formula\<close>
kleing@59190
   141
kleing@59190
   142
lemma quartic:
wenzelm@63054
   143
  "(y::real)^3 - b * y^2 + (a * c - 4 * d) * y - a^2 * d + 4 * b * d - c^2 = 0 \<and>
wenzelm@63054
   144
    R^2 = a^2 / 4 - b + y \<and>
wenzelm@63054
   145
    s^2 = y^2 - 4 * d \<and>
wenzelm@63054
   146
    (D^2 = (if R = 0 then 3 * a^2 / 4 - 2 * b + 2 * s
wenzelm@63054
   147
                     else 3 * a^2 / 4 - R^2 - 2 * b + (4 * a * b - 8 * c - a^3) / (4 * R))) \<and>
wenzelm@63054
   148
    (E^2 = (if R = 0 then 3 * a^2 / 4 - 2 * b - 2 * s
wenzelm@63054
   149
                     else 3 * a^2 / 4 - R^2 - 2 * b - (4 * a * b - 8 * c - a^3) / (4 * R)))
kleing@59190
   150
  \<Longrightarrow> x^4 + a * x^3 + b * x^2 + c * x + d = 0 \<longleftrightarrow>
kleing@59190
   151
      x = -a / 4 + R / 2 + D / 2 \<or>
kleing@59190
   152
      x = -a / 4 + R / 2 - D / 2 \<or>
kleing@59190
   153
      x = -a / 4 - R / 2 + E / 2 \<or>
kleing@59190
   154
      x = -a / 4 - R / 2 - E / 2"
wenzelm@63054
   155
  apply (cases "R = 0")
wenzelm@63054
   156
   apply (simp_all add: field_simps divide_minus_left[symmetric] del: divide_minus_left)
wenzelm@63054
   157
   apply algebra
wenzelm@63054
   158
  apply algebra
wenzelm@63054
   159
  done
kleing@59190
   160
nipkow@62390
   161
end