cubic-equation.egi

-- Cubic Formula (Cardano's formula)

declare symbol x, a, b, c, d, p, q

def cubicFormula : MathExpr -> MathExpr -> (MathExpr, MathExpr, MathExpr) := cF

def cF (f: MathExpr) (x: MathExpr) : (MathExpr, MathExpr, MathExpr) :=
  match coefficients f x as list mathExpr with
    | [$a_0, $a_1, $a_2, $a_3] -> cF' a_3 a_2 a_1 a_0

def cF' (a: MathExpr) (b: MathExpr) (c: MathExpr) (d: MathExpr) : (MathExpr, MathExpr, MathExpr) :=
  match (a, b, c, d) as (mathExpr, mathExpr, mathExpr, mathExpr) with
    | (#1, #0, $p, $q) ->
      let (s1, s2) := (2)#(rt 3 $1, rt 3 $2) (qF' 1 (27 * q) ((-27) * p ^ 3))
       in ((s1 + s2) / 3, (w ^ 2 * s1 + w * s2) / 3, (w * s1 + w ^ 2 * s2) / 3)
    | (#1, _, _, _) ->
      (3)#($1 - b / 3, $2 - b / 3, $3 - b / 3)
        (withSymbols [x, y]
          cF (substitute [(x, y - b / 3)] (x ^ 3 + b * x ^ 2 + c * x + d)) y)
    | (_, _, _, _) -> cF' 1 (b / a) (c / a) (d / a)

def w : MathExpr := ((-1) + i * sqrt 3) / 2

-- Solution for x^3 + p*x + q = 0 (depressed cubic)
(3)#$1 (cF (x ^ 3 + p * x + q) x)

-- Verify: (x-1)(x-2)(x-3) = x^3 - 6x^2 + 11x - 6
assertEqual "cubic (x-1)(x-2)(x-3)"
  (cF ((x - 1) * (x - 2) * (x - 3)) x)
  (2, (- sqrt 3 * rt 3 (3 * sqrt 3) + 6) / 3, (sqrt 3 * rt 3 (3 * sqrt 3) + 6) / 3)

-- General form
(3)#$1 (cF (a * x ^ 3 + b * x ^ 2 + c * x + d) x)