100/cubic.ml

   1 (* ========================================================================= *)
   2 (* Cubic formula.                                                            *)
   3 (* ========================================================================= *)
   4
   5 needs "Complex/complex_transc.ml";;
   6
   7 prioritize_complex();;
   8
   9 (* ------------------------------------------------------------------------- *)
  10 (* Define cube roots (it doesn't matter which one we choose here)            *)
  11 (* ------------------------------------------------------------------------- *)
  12
  13 let ccbrt = new_definition
  14  `ccbrt(z) = if z = Cx(&0) then Cx(&0) else cexp(clog(z) / Cx(&3))`;;
  15
  16 let CCBRT = prove
  17  (`!z. ccbrt(z) pow 3 = z`,
  18   GEN_TAC THEN REWRITE_TAC[ccbrt] THEN COND_CASES_TAC THEN
  19   ASM_REWRITE_TAC[] THENL [CONV_TAC COMPLEX_RING; ALL_TAC] THEN
  20   REWRITE_TAC[COMPLEX_FIELD `z pow 3 = z * z * z:complex`; GSYM CEXP_ADD] THEN
  21   REWRITE_TAC[COMPLEX_FIELD `z / Cx(&3) + z / Cx(&3) + z / Cx(&3) = z`] THEN
  22   ASM_SIMP_TAC[CEXP_CLOG]);;
  23
  24 (* ------------------------------------------------------------------------- *)
  25 (* The reduction to a simpler form.                                          *)
  26 (* ------------------------------------------------------------------------- *)
  27
  28 let CUBIC_REDUCTION = COMPLEX_FIELD
  29   `~(a = Cx(&0)) /\
  30    x = y - b / (Cx(&3) * a) /\
  31    p = (Cx(&3) * a * c - b pow 2) / (Cx(&9) * a pow 2) /\
  32    q = (Cx(&9) * a * b * c - Cx(&2) * b pow 3 - Cx(&27) * a pow 2 * d) /
  33        (Cx(&54) * a pow 3)
  34    ==> (a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
  35         y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0))`;;
  36
  37 (* ------------------------------------------------------------------------- *)
  38 (* The solutions of the special form.                                        *)
  39 (* ------------------------------------------------------------------------- *)
  40
  41 let CUBIC_BASIC = COMPLEX_FIELD
  42  `!i t s.
  43     s pow 2 = q pow 2 + p pow 3 /\
  44     (s1 pow 3 = if p = Cx(&0) then Cx(&2) * q else q + s) /\
  45     s2 = --s1 * (Cx(&1) + i * t) / Cx(&2) /\
  46     s3 = --s1 * (Cx(&1) - i * t) / Cx(&2) /\
  47     i pow 2 + Cx(&1) = Cx(&0) /\
  48     t pow 2 = Cx(&3)
  49     ==> if p = Cx(&0) then
  50           (y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0) <=>
  51            y = s1 \/ y = s2 \/ y = s3)
  52         else
  53           ~(s1 = Cx(&0)) /\
  54           (y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0) <=>
  55                (y = s1 - p / s1 \/ y = s2 - p / s2 \/ y = s3 - p / s3))`;;
  56
  57 (* ------------------------------------------------------------------------- *)
  58 (* Explicit formula for the roots (doesn't matter which square/cube roots).  *)
  59 (* ------------------------------------------------------------------------- *)
  60
  61 let CUBIC = prove
  62  (`~(a = Cx(&0))
  63    ==> let p = (Cx(&3) * a * c - b pow 2) / (Cx(&9) * a pow 2)
  64        and q = (Cx(&9) * a * b * c - Cx(&2) * b pow 3 - Cx(&27) * a pow 2 * d) /
  65                (Cx(&54) * a pow 3) in
  66        let s = csqrt(q pow 2 + p pow 3) in
  67        let s1 = if p = Cx(&0) then ccbrt(Cx(&2) * q) else ccbrt(q + s) in
  68        let s2 = --s1 * (Cx(&1) + ii * csqrt(Cx(&3))) / Cx(&2)
  69        and s3 = --s1 * (Cx(&1) - ii * csqrt(Cx(&3))) / Cx(&2) in
  70        if p = Cx(&0) then
  71          a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
  72             x = s1 - b / (Cx(&3) * a) \/
  73             x = s2 - b / (Cx(&3) * a) \/
  74             x = s3 - b / (Cx(&3) * a)
  75        else
  76          ~(s1 = Cx(&0)) /\
  77          (a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
  78              x = s1 - p / s1 - b / (Cx(&3) * a) \/
  79              x = s2 - p / s2 - b / (Cx(&3) * a) \/
  80              x = s3 - p / s3 - b / (Cx(&3) * a))`,
  81   DISCH_TAC THEN REPEAT LET_TAC THEN
  82   ABBREV_TAC `y = x + b / (Cx(&3) * a)` THEN
  83   SUBGOAL_THEN
  84    `a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
  85     y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0)`
  86   SUBST1_TAC THENL
  87    [MATCH_MP_TAC CUBIC_REDUCTION THEN ASM_REWRITE_TAC[] THEN
  88     EXPAND_TAC "y" THEN CONV_TAC COMPLEX_RING;
  89     ALL_TAC] THEN
  90   ONCE_REWRITE_TAC[COMPLEX_RING `x = a - b <=> x + b = (a:complex)`] THEN
  91   ASM_REWRITE_TAC[] THEN MATCH_MP_TAC CUBIC_BASIC THEN
  92   MAP_EVERY EXISTS_TAC
  93    [`ii`; `csqrt(Cx(&3))`; `csqrt (q pow 2 + p pow 3)`] THEN
  94   ASM_REWRITE_TAC[] THEN REPEAT CONJ_TAC THENL
  95    [ASM_MESON_TAC[CSQRT];
  96     ASM_MESON_TAC[CCBRT];
  97     MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_RING;
  98     ASM_MESON_TAC[CSQRT]]);;