There is a simpler solution than Ferrari's method.

Assume the equation: ax^4 + bx^3 + cx^2 + dx + e = 0

Simplify by dividing by a (this should always be done):

x^4 + Bx^3 + Cx^2 + Dx + E = 0

where
B = b/a
C = c/a
D = d/a
E = e/a

Solve this cubic equation for z: z^3 + 2Cz^2 + (C^2 - 4E)z - D^2 = 0

Solve the equations:
p = sqrt(z)
q = (c + z - d/p)/2
r = -sqrt(z)
s = (c + z + d/p)/2

Solve these quadratic equations to obtain the four roots:

x^2 + px + q = 0
x^2 + rx + s = 0

You'll have to understand that Ferrari discovered the first method, and it was complicated.

This is simpler, but I can't guarantee that it will have better performance, since you'll have to calculate the cubic and quadratic solutions.

If you're interested, the theory behind this can be found here.

I've been dealing with a few numerical stability issues off and on. This may be a better solution because of that. I'm curious why the term B never shows up in the equations. It seems that B should matter!

Ups. Thanks for the warning, pTymN.

This is the correct method.

ax^4 + bx^3 + cx^2 + dx + e = 0 (=) Divide by a

x^4 + Bx^3 + Cx^2 + Dx + E = 0

B = b/a
C = c/a
D = d/a
E = e/a

(depressed quartic, like in Ferrari's method)
alpha = I = -3(B^2) / 8 + C
beta = J = (B^3)/8 - BC/2 + D
gamma = K = -3(B^4)/256 + C(B^2)/16 - BD/4 + E

Solve equation for z (one solution is enough):
z^3 + 2Iz^2 + (I^2 - 4K)z - J^2 = 0

p = sqrt(z)
r = -p
q = (I + z - J/p)/2
s = (I + z + J/p)/2

Solve (there are four u's, two for each equation):
u^2 + pu + q = 0
u^2 + ru + s = 0

Solution (x[]):
for each root u, x = u - B/4

[Edited by - Silver Phoenix on June 9, 2007 6:20:47 PM]

Quote:Original post by johnnyBravo
I just kinda learned how to do the cubic equation, but I don't understand how there are only 3 roots, because I have the two +- from the quadratic equation,

but from that I'd assume you would be able to use the ++, --, +- and -+

A polynomial of degree N will always have no more than N roots. So a cubic equation cannot have more than 3 roots.

You only need one root from the three that the cubic equation gives. The final solution will always be equal (not taking precision loss into account).

For each different cubic solution there will be different p, q, r and s.
These different solutions give different second degree polynomials, but there's an invariants:

(u^2 + pu + q)*(u^2 + ru + s) = 0 (=)

(=) u^4 + (p + r)u^3 + (pr + q + s)u^2 + (ps + qr)u + qs = 0

For simplicity we will just focus on the simpler monomials: (p + r)u^3 and qs

If we take any two solutions from the cubic, then if one solution decreases p, the other will have to increase r, because (p + r) is always constant.
Actually, by definition r = -p => p + r = 0 (depressed quartic).

The same goes for qs. If two solution give q with opposite signs , then s will also have opposite signs for the two solutions.

Also remember that we solved a system of equations with p, q, r and s on the left side.

Quote:Original post by johnnyBravo
btw x as in ax^3 + bx^2 + cx + d

x is the final solution. The solution for:

ax^4 + bx^3 + cx^2 + dx + e = 0

(or x^4 + Bx^3 + Cx^2 + Dx + E = 0, since they have the same roots)

Some notes:
There was a small but fatal mistake in the z equation. I already corrected it.

Compare edited with:

Quote:Original
z^3 + 2Iz^2 + (I^2 - 4K)z + J^2 = 0

It is possible to simplify the cubic if K or J is zero.

If z is zero => p = sqrt(z) = 0 => J/p = J/0 (divide by zero)

Slight nitpick, but a cubic polynomial can have less roots than 3 in the real numbers (3 in complex of course), and a polynomial of degree N can have MORE than N roots in a non-integral domain (such as the residue class ring Z(N) of integers modulo N, for non-prime N)

In Z(8),

x^2 - 1 = 0 has roots 1, 3, 5, and 7 (mod 8).

Of course, Z(8) has elements that are proper divisors of zero, i.e. a*b = 0 although a != 0 and b != 0. (namely, 2 and 4 (mod 8)).

bool solveQuadraticOther(double a, double b, double c, double &root){	if(a == 0.0 || abs(a/b) < 1.0e-4)	{		if(abs(b) < 1.0e-4) 			return false;		else		{			root = -c/b;			return true;		}	}	double discriminant = b * b - 4.0 * a * c;	if(discriminant >= 0.0)	{		discriminant = sqrt(discriminant);		root = (b - discriminant) * -0.5 / a;		return true;	}	return false;}//----------------------------------------------------------------------------bool solveQuartic(double a, double b, double c, double d, double e, double &root){    // I switched to this method, and it seems to be more numerically stable.    // http://www.gamedev.net/community/forums/topic.asp?topic_id=451048 	// When a or (a and b) are magnitudes of order smaller than C,D,E	// just ignore them entirely. This seems to happen because of numerical	// inaccuracies of the line-circle algorithm. I wanted a robust solver,	// so I put the fix here instead of there.	if(a == 0.0 || abs(a/b) < 1.0e-5 || abs(a/c) < 1.0e-5 || abs(a/d) < 1.0e-5)		return solveCubic(b, c, d, e, root);	double B = b/a, C = c/a, D = d/a, E = e/a;	double BB = B*B;	double I = -3.0*BB*0.125 + C;	double J = BB*B*0.125 - B*C*0.5 + D;	double K = -3*BB*BB/256.0 + C*BB/16.0 - B*D*0.25 + E;	double z;	bool foundRoot2 = false, foundRoot3 = false, foundRoot4 = false, foundRoot5 = false;	if(solveCubic(1.0, I+I, I*I - 4*K, -(J*J), z))	{		double value = z*z*z + z*z*(I+I) + z*(I*I - 4*K) - J*J;		double p = sqrt(z);		double r = -p;		double q = (I + z - J/p)*0.5;		double s = (I + z + J/p)*0.5;		bool foundRoot = false, foundARoot;		double aRoot;		foundRoot = solveQuadratic(1.0, p, q, root);		root -= B/4.0;		foundARoot = solveQuadratic(1.0, r, s, aRoot);		aRoot -= B/4.0;		if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 			|| root < 0.0)) || (!foundRoot && foundARoot)) 		{			root = aRoot;			foundRoot = true;		}		foundARoot = solveQuadraticOther(1.0, p, q, aRoot);		aRoot -= B/4.0;		if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 			|| root < 0.0)) || (!foundRoot && foundARoot)) 		{			root = aRoot;			foundRoot = true;		}		foundARoot = solveQuadraticOther(1.0, r, s, aRoot);		aRoot -= B/4.0;		if((foundRoot && foundARoot && ((aRoot < root && aRoot >= 0.0) 			|| root < 0.0)) || (!foundRoot && foundARoot)) 		{			root = aRoot;			foundRoot = true;		}		return foundRoot;	}	return false;}

This method is MUCH more stable than Ferrari's method. I highly encourage its use.

[Edited by - pTymN on June 13, 2007 11:14:04 PM]