cubic solver

[Albert Chan]

We can use same trick (post #8) to depress quartic, a*x^4 + b*x^3 + c*x^2 + d*x + e = 0

{a, b, c, d, e} = {(A = 4a), 4b, 4c, 4d, 4e} = {1, 4b, 4Ac, 4A²d, 4A³e} / A

Synthetic Division to force coefficient 4b to 0

Code:

-b>   1   4b    4Ac          4A²d               4A³e
-b>   1   3b    4Ac-3b²      4A²d-b*(4Ac-3b²)   4A³e-b*(4A²d-b*(4Ac-3b²)) = E
-b>   1   2b    4Ac-5b²      4A²d-b*(8Ac-8b²) = D
-b>   1    b    4Ac-6b² = C
      1    0

{a, b, c, d, e} = ({1, 0, C, D, E} - b) / A

y^4 + C*y^2 + D*y + E = 0
(y^2+s)^2 = (2s-C)*y^2 - D*y + (s^2-E)      // Add an offset s, making LHS perfect square

Let k² = 2s-C ≠ 0, and force RHS as perfect square too.

{(k² = 2s-C), -D, (s²-E)} = {1, -D, (D²/4 = (s²-E)*k²)} / k²

D² = 4*(s²-E)*k² = ((k²+C)² - 4E) * k²

(k²)^3 + 2*C*(k²)^2 + (C²-4E)*(k²) - D² = 0

The reason we solve for k² instead of s is to easily check k²≠0
k² roots all zero ⇒ {C,D,E} all zeroes ⇒ y roots all zero too.

(y²+s)² = k² * (y-(D/2)/k²)²      // let t = D/2/k
(y²+s) = ±(k*y - t)

(y^2 - ky + s+t) * (y^2 + ky + s-t) = 0

(s±t), one of them may cause huge cancellation, we only use the good one.
The other can be derived by product of all roots = (s+t) * (s-t) = E

Last step is simply map y back to x: (x roots) = ((y roots) - b) / A

Code:

function quartic(a,b,c,d,e) -- a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
    a = 4*a
    c, d, e = 4*a*c, 4*a^2*d, 4*a^3*e
    e = e - b*(d - b*(c - 3*b*b))
    d = d - b*(c+c - 8*b*b)
    c = c - 6*b*b           -- y coefs = 1, 0, c, d, e
    local k2 = {cubic(1, 2*c, c*c-4*e, -d*d)}
    table.sort(k2, function(a,b) return I.abs(a) > I.abs(b) end)
    if I.abs(k2[1])==0 then c=I(-b/a); return c,c,c,c end
    local k = I.sqrt(k2[1]) -- (y^2+s)^2 = (k*y-t)^2
    c, d = k2[1]+c, d/k     -- 2*s, 2*t
    c, d = (c+d)/2, (c-d)/2 -- quadratic constant term
    if I.abs(c) < I.abs(d) then c=d; k=-k end
    local y1, y2 = quadratic(1, -k,  c)
    local y3, y4 = quadratic(1, k, e/c)
    return (y1-b)/a, (y2-b)/a, (y3-b)/a, (y4-b)/a
end

lua> quartic(1,2,3,4,5)
(0.28781547955764775+1.4160930801719083*I)
(-1.287815479557648+0.8578967583284901*I)
(0.2878154795576481-1.4160930801719085*I)
(-1.287815479557648-0.8578967583284899*I)