(* ========================================================================= *)
(* Cubic formula.                                                            *)
(* ========================================================================= *)

needs "Complex/complex_transc.ml";;

prioritize_complex();;

(* ------------------------------------------------------------------------- *)
(* Define cube roots (it doesn't matter which one we choose here)            *)
(* ------------------------------------------------------------------------- *)

let ccbrt = new_definition
 `ccbrt(z) = if z = Cx(&0) then Cx(&0) else cexp(clog(z) / Cx(&3))`;;
let CCBRT = 
prove
 (`!z. ccbrt(z) pow 3 = z`,

  GEN_TAC THEN REWRITE_TAC[ccbrt] THEN COND_CASES_TAC THEN
  ASM_REWRITE_TAC[] THENL [CONV_TAC COMPLEX_RING; ALL_TAC] THEN
  REWRITE_TAC[COMPLEX_FIELD `z pow 3 = z * z * z:complex`; GSYM CEXP_ADD] THEN
  REWRITE_TAC[COMPLEX_FIELD `z / Cx(&3) + z / Cx(&3) + z / Cx(&3) = z`] THEN
  ASM_SIMP_TAC[CEXP_CLOG]);;
(* ------------------------------------------------------------------------- *)
(* The reduction to a simpler form.                                          *)
(* ------------------------------------------------------------------------- *)

let CUBIC_REDUCTION = COMPLEX_FIELD
  `~(a = Cx(&0)) /\
   x = y - b / (Cx(&3) * a) /\
   p = (Cx(&3) * a * c - b pow 2) / (Cx(&9) * a pow 2) /\
   q = (Cx(&9) * a * b * c - Cx(&2) * b pow 3 - Cx(&27) * a pow 2 * d) /
       (Cx(&54) * a pow 3)
   ==> (a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
        y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0))`;;

(* ------------------------------------------------------------------------- *)
(* The solutions of the special form.                                        *)
(* ------------------------------------------------------------------------- *)

let CUBIC_BASIC = COMPLEX_FIELD
 `!i t s.
    s pow 2 = q pow 2 + p pow 3 /\
    (s1 pow 3 = if p = Cx(&0) then Cx(&2) * q else q + s) /\
    s2 = --s1 * (Cx(&1) + i * t) / Cx(&2) /\
    s3 = --s1 * (Cx(&1) - i * t) / Cx(&2) /\
    i pow 2 + Cx(&1) = Cx(&0) /\
    t pow 2 = Cx(&3)
    ==> if p = Cx(&0) then 
          (y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0) <=>
           y = s1 \/ y = s2 \/ y = s3)
        else
          ~(s1 = Cx(&0)) /\
          (y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0) <=>
               (y = s1 - p / s1 \/ y = s2 - p / s2 \/ y = s3 - p / s3))`;;

(* ------------------------------------------------------------------------- *)
(* Explicit formula for the roots (doesn't matter which square/cube roots).  *)
(* ------------------------------------------------------------------------- *)

let CUBIC = 
prove
 (`~(a = Cx(&0))
   ==> let p = (Cx(&3) * a * c - b pow 2) / (Cx(&9) * a pow 2)
       and q = (Cx(&9) * a * b * c - Cx(&2) * b pow 3 - Cx(&27) * a pow 2 * d) /
               (Cx(&54) * a pow 3) in
       let s = csqrt(q pow 2 + p pow 3) in
       let s1 = if p = Cx(&0) then ccbrt(Cx(&2) * q) else ccbrt(q + s) in
       let s2 = --s1 * (Cx(&1) + ii * csqrt(Cx(&3))) / Cx(&2)
       and s3 = --s1 * (Cx(&1) - ii * csqrt(Cx(&3))) / Cx(&2) in
       if p = Cx(&0) then 
         a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
            x = s1 - b / (Cx(&3) * a) \/
            x = s2 - b / (Cx(&3) * a) \/               
            x = s3 - b / (Cx(&3) * a)
       else
         ~(s1 = Cx(&0)) /\                          
         (a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
             x = s1 - p / s1 - b / (Cx(&3) * a) \/
             x = s2 - p / s2 - b / (Cx(&3) * a) \/
             x = s3 - p / s3 - b / (Cx(&3) * a))`,

  DISCH_TAC THEN REPEAT LET_TAC THEN
  ABBREV_TAC `y = x + b / (Cx(&3) * a)` THEN
  SUBGOAL_THEN 
   `a * x pow 3 + b * x pow 2 + c * x + d = Cx(&0) <=>
    y pow 3 + Cx(&3) * p * y - Cx(&2) * q = Cx(&0)`
  SUBST1_TAC THENL
   [MATCH_MP_TAC CUBIC_REDUCTION THEN ASM_REWRITE_TAC[] THEN
    EXPAND_TAC "y" THEN CONV_TAC COMPLEX_RING;
    ALL_TAC] THEN
  ONCE_REWRITE_TAC[COMPLEX_RING `x = a - b <=> x + b = (a:complex)`] THEN
  ASM_REWRITE_TAC[] THEN MATCH_MP_TAC CUBIC_BASIC THEN
  MAP_EVERY EXISTS_TAC 
   [`ii`; `csqrt(Cx(&3))`; `csqrt (q pow 2 + p pow 3)`] THEN
  ASM_REWRITE_TAC[] THEN REPEAT CONJ_TAC THENL
   [ASM_MESON_TAC[CSQRT];
    ASM_MESON_TAC[CCBRT];
    MP_TAC COMPLEX_POW_II_2 THEN CONV_TAC COMPLEX_RING;
    ASM_MESON_TAC[CSQRT]]);;