| home | send comment | send link | add bookmark |

Exact solutions of cubic polynomial equations

by Stephen R. Schmitt


Algorithm

An exact solution of the cubic polynomial equation:

x3 + a·x2 + b·x + c = 0
was first published by Gerolamo Cardano (1501-1576) in his treatise, Ars Magna. He did not discoverer of the solution; a professor of mathematics at the University of Bologna named Scipione del Ferro (ca. 1465-1526) is credited as the first to find an exact solution. In the years since, several improvements to the original solution have been discovered.

Zeno source code

Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.

const pi : real := 2*arccos( 0 )

program

    var a, b, c, z1, z2, z3, im : real

    put "solve: x^3 + a*x^2 + b*x + c = 0"
    put "a b c"...
    get a, b, c
    
    im := cubic( z1, z2, z3, a, b, c )

    if im ~= 0.0 then
        put "{ ", z1:6:4, ", ", 
                  z2:6:4, " + ", im:6:4, "j", ", ", 
                  z3:6:4, " - ", im:6:4, "j", " }"
    else
        put "{ ", z1:6:4, ", ", z2:6:4, ", ", z3:6:4, " }"
    end if

end program

% compute real or complex roots of cubic polynomial
function cubic( var z1, z2, z3 : real, a2, a1, a0 : real ) : real

    var Q, R, D, S, T : real
    var im, th : real

    Q := (3*a1 - a2^2)/9
    R := (9*a1*a2 - 27*a0 - 2*a2^3)/54
    D := Q^3 + R^2                          % polynomial discriminant

    if (D >= 0) then                        % complex or duplicate roots

        S := sgn(R + sqrt(D))*rabs(R + sqrt(D))^(1/3)
        T := sgn(R - sqrt(D))*rabs(R - sqrt(D))^(1/3)

        z1 := -a2/3 + (S + T)               % real root
        z2 := -a2/3 - (S + T)/2             % real part of complex root
        z3 := -a2/3 - (S + T)/2             % real part of complex root
        im := rabs(sqrt(3)*(S - T)/2)       % complex part of root pair

    else                                    % distinct real roots

        th := arccos(R/sqrt( -Q^3))
        
        z1 := 2*sqrt(-Q)*cos(th/3) - a2/3
        z2 := 2*sqrt(-Q)*cos((th + 2*pi)/3) - a2/3
        z3 := 2*sqrt(-Q)*cos((th + 4*pi)/3) - a2/3
        im := 0

    end if

    return im                               % imaginary part

end function

% sign of number
function sgn( x : real ) : int
    if x < 0.0 then 
        return -1
    end if
    return 1
end function

% absolute value of real number
function rabs( x : real ) : real
    if x < 0.0 then x := -x end if
    return x
end function

Sample output

solve: x^3 + a*x^2 + b*x + c = 0
a b c? -2 1 -2
{ 2.0000, 0.0000 + 1.0000j, 0.0000 - 1.0000j }

Reference

Wolfram Research : Cubic Equation

AbCd Classics - free on-line books


Copyright © 2005, Stephen R. Schmitt