# Solution of quartic equation

Is there a open C-implementation for the solution of quartic equations:

$$ax⁴+bx³+cx²+dx+e=0$$

I am thinking of an implementation of Ferrari's solution. On Wikipedia I read that the solution is computational stable only for some of the possible sign combinations of the coefficients. But maybe I am lucky ... I got a pragmatic solution by solving analytically using a computer algebra system and exporting to C. But if there is a tested implementation I would prefer to use this. I search a fast method and prefer not to use a general root finder.

I need only real solutions.

-
Do you need all the (real) solutions simultaneously? As GertVdE says below, if you have stability issues with a closed form solution, there isn't really a good reason to not use some root-finding algorithm. – Godric Seer Sep 4 '12 at 11:16
I find it funny that this was tagged nonlinear-algebra, since you could simply calculate the eigenvalues of the companion matrix, which is already in Hessenberg form and applying QR sweeps would be pretty simple. – Victor Liu Sep 4 '12 at 19:28

I would strongly advice against using closed form solutions since they tend to be numerically very unstable. You need to take extreme care in the way and order of your evaluations of the discriminant and other parameters.

The classical example is the one for the quadratic equation $ax^2+bx+c=0$. Calculating the roots as $$x_{1,2} = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$$ will get you into trouble for polynomials where $b\gg 4ac$ since then you get cancellation in the numerator. You need to calculate $$x_1 = \frac{-(b+sign(b)\sqrt{b^2-4ac})}{2a}; x_2 = \frac{c}{a}\frac{1}{x_1}$$.

Higham in his masterpiece "Accuracy and Stability of Numerical Algorithms" (2nd ed, SIAM) uses a direct search method to find coefficients of a cubic polynomial for which the classical analytical cubic solution gives very inaccurate results. The example he gives is $[a,b,c] = [1.732, 1, 1.2704]$. For this polynomial the roots are well-separated and hence the problem is not ill-conditioned. However, if he calculates the roots using the analytical approach, and evaluates the polynomial in these roots, he obtains a residue of $\mathcal{O}(10^{-2})$ while using a stable standard method (the companion matrix method), the residue is of order $\mathcal{O}(10^{-15})$. He proposes a slight modification to the algorithm, but even then, he can find a set of coefficients leading to residues of $\mathcal{O}(10^{-11})$ which is definitely not good. See p480-481 of the above mentioned book.

In your case, I would apply Bairstow's method. It uses an iterative combination of Newton iteration on quadratic forms (and then the roots of the quadratic are solved) and deflation. It is easily implemented and there are even some implementations available on the web.

-
Could you please explain what do you mean by "I would strongly advice against using closed form solutions since they tend to be numerically very unstable. ". Does this apply to 4th degree polynomials only or is this a general rule? – Emmad Kareem Sep 4 '12 at 21:46
@EmmadKareem I've updated my answer above. – GertVdE Sep 5 '12 at 8:01
@+1, very good, thank you. – Emmad Kareem Sep 5 '12 at 21:11

See these:

-
 Using this code on the polynomial with coefficients given in my answer, I find the following: $x_1 = -1.602644912244132e+00$, which has a relative error of $\mathcal{O}(10^{-8})$ compared to the real root (calculated using Octave's roots command which uses the companion matrix method). It has a residue of $\mathcal{O}(10^{-7})$ while the companion matrix method root has a residue of $\mathcal{O}(10^{-15})$. Up to you whether this is good enough (for computer graphics it might be, for some other applications, it won't be) – GertVdE Sep 6 '12 at 11:20