SetHelp ("QuadraticFormula", "equation_solving",
         "Find roots of a quadratic polynomial (given as vector of coefficients)")
function QuadraticFormula(p) = (
	if not IsPoly(p) or elements(p) != 3 or p@(3) == 0 then
		(error("QuadraticFormula: argument 1 must be a quadratic polynomial");bailout);
	[ (- p@(2) + sqrt(p@(2)^2 - 4*p@(1)*p@(3))) / (2 * p@(3)) 
	 (- p@(2) - sqrt(p@(2)^2 - 4*p@(1)*p@(3))) / (2 * p@(3))]
)
protect ("QuadraticFormula")

# see http://en.wikipedia.org/wiki/Cubic_equation
SetHelp ("CubicFormula", "equation_solving",
         "Find roots of a cubic polynomial (given as vector of coefficients)")
function CubicFormula(pol) = (
	if not IsPoly(pol) or elements(pol) != 4 or pol@(4) == 0 then
		(error("CubicFormula: argument 1 must be a cubic polynomial");bailout);
	c := pol@(1) / pol@(4);
	b := pol@(2) / pol@(4);
	a := pol@(3) / pol@(4);

    p := b - (a^2)/3;
    q := c + (2*a^3 - 9*a*b)/27;
    
    if p == 0 then
    	if q == 0 then
    		return [-a/3,-a/3,-a/3]
    	else
    		u := q^(1/3)
    else if q == 0 then
    	u := sqrt(p/3)
    else
    	u := ( q/2 + sqrt(((q^2) / 4) + (p^3)/27) )^(1/3);

    omega1 := -(1/2)+1i*sqrt(3)/2;
    omega2 := -(1/2)-1i*sqrt(3)/2;
    
    # If real coefficients and the discriminant is negative
    # or if q == 0, then we have all real roots 
    if (IsReal(a) and IsReal(b) and IsReal(c)) then (
	DELTA := - a^2*b^2 + 4*b^3 + 4*a^3*c + 27*c^2 - 18*a*b*c;
    	# If the discriminant is negative
    	# or if q == 0, then we have all real roots 
	if DELTA <= 0 or q == 0 then
		[RealPart (p/(3*u) - u - a/3)
		 RealPart (p/(3*u*omega1) - u*omega1 - a/3)
		 RealPart (p/(3*u*omega2) - u*omega2 - a/3)]
	# if the discriminant is positive then we have one real root
	else #if DELTA > 0 then
		[RealPart (p/(3*u) - u - a/3)
		 p/(3*u*omega1) - u*omega1 - a/3
		 p/(3*u*omega2) - u*omega2 - a/3]
    ) else
    	[p/(3*u) - u - a/3
     	 p/(3*u*omega1) - u*omega1 - a/3
     	 p/(3*u*omega2) - u*omega2 - a/3]
)
protect ("CubicFormula")

## see planetmath
SetHelp ("QuarticFormula", "equation_solving",
         "Find roots of a quartic polynomial (given as vector of coefficients)")
function QuarticFormula(pol) = (
	if not IsPoly(pol) or elements(pol) != 5 or pol@(5) == 0 then
		(error("QuarticFormula: argument 1 must be a quartic polynomial");bailout);
	d := pol@(1) / pol@(5);
	c := pol@(2) / pol@(5);
	b := pol@(3) / pol@(5);
	a := pol@(4) / pol@(5);
	
	# resolvent cubic
	t := CubicFormula([c^2+a^2*d-a*b*c,b^2+a*c-4*d,-2*b,1]);

	# compute r1 * r2 and r3 * r4
	tt := t@(2)+t@(3)-t@(1);
	A := sqrt(tt^2 - 16*d);
	r1r2 := ((tt) + A)/4;
	r3r4 := ((tt) - A)/4;

	# note that it is hard to tell apart r1 + r2 and r3 + r4,
	# hence we'll try both below
	A := sqrt(a^2 - 4*t@(1));
	r1pr2 := (-a + A)/2;
	r3pr4 := (-a - A)/2;

	# compare the minus of the correct symetric function of the roots
	# to coeff c and see which is better and use that
	if (|r1r2*r3pr4 + r1pr2*r3r4 + c| > |r1r2*r1pr2 + r3pr4*r3r4 + c|) then (
		r1pr2 := (-a - A)/2;
		r3pr4 := (-a + A)/2
	);
    	
	A1 := sqrt(r1pr2^2 -4*r1r2);
	A2 := sqrt(r3pr4^2 -4*r3r4);
	[(r1pr2 + A1)/2
	 (r1pr2 - A1)/2
	 (r3pr4 + A2)/2
	 (r3pr4 - A2)/2]
)
protect ("QuarticFormula")

SetHelp ("PolynomialRoots", "equation_solving",
         "Find roots of a polynomial (given as vector of coefficients)")
function PolynomialRoots(p) = (
	if not IsPoly(p) then
		(error("PolynomialRoots: argument 1 must be a polynomial");bailout);
	p = TrimPoly (p);
	if elements(p) < 2 or elements(p) > 5 then
		(error("PolynomialRoots: Solving for polynomials only of degrees 1 through 4 are implemented");bailout);
	
	if elements(p) == 2 then
		[-p@(1)/p@(2)]
	else if elements(p) == 3 then
		QuadraticFormula (p)
	else if elements(p) == 4 then
		CubicFormula (p)
	else # if elements(p) == 5 then
		QuarticFormula (p)
)
protect ("PolynomialRoots")