quartic-equation.egi

-- Quartic Formula (Ferrari's method)

declare symbol x, y

def quarticFormula : MathExpr -> MathExpr -> (MathExpr, MathExpr, MathExpr, MathExpr) := qtF

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

def qtF' (a: MathExpr) (b: MathExpr) (c: MathExpr) (d: MathExpr) (e: MathExpr) : (MathExpr, MathExpr, MathExpr, MathExpr) :=
  match (a, b, c, d, e) as
    (mathExpr, mathExpr, mathExpr, mathExpr, mathExpr) with
    | (#1, #0, $p, #0, $q) ->
      let (s1, s2) := qF' 1 p q
          (r1, r2) := qF' 1 0 (- s1)
          (r3, r4) := qF' 1 0 (- s2)
       in (r1, r2, r3, r4)
    | (#1, #0, $p, $q, $r) ->
      let u := (3)#$1
                 (withSymbols [u]
                   cF (u * (p + u) ^ 2 + (-4) * r * u + (- (q ^ 2))) u)
          (r1, r2) := qF (y ^ 2 + (p + u) / 2 + sqrt u * (y - q / (2 * u))) y
          (r3, r4) := qF
                        (y ^ 2 + (p + u) / 2 + (- sqrt u) * (y - q / (2 * u)))
                        y
       in (r1, r2, r3, r4)
    | (#1, _, _, _, _) ->
      (4)#($1 - b / 4, $2 - b / 4, $3 - b / 4, $4 - b / 4)
        (withSymbols [x, y]
          qtF
            (substitute
               [(x, y - b / 4)]
               (x ^ 4 + b * x ^ 3 + c * x ^ 2 + d * x + e))
            y)
    | (_, _, _, _, _) -> qtF' 1 (b / a) (c / a) (d / a) (e / a)

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

-- Verify: (x-1)(x-2)(x-3)(x-4) should give roots 1, 2, 3, 4
assertEqual "quartic (x-1)(x-2)(x-3)(x-4)"
  (qtF ((x - 1) * (x - 2) * (x - 3) * (x - 4)) x)
  (4, 1, 3, 2)