10
$\begingroup$

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.

| cite | improve this question | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Godric Seer Sep 4 '12 at 11:16
  • 3
    $\begingroup$ 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. $\endgroup$ – Victor Liu Sep 4 '12 at 19:28
  • 2
    $\begingroup$ Have a look at the cubic/quartic solvers published in ACM TOMS (Algorithm 954). Code that makes it into that journal is usually of very high quality. The paper itself is behind a paywall but the code can be downloaded from this link. $\endgroup$ – GoHokies Jan 28 '17 at 8:40
  • $\begingroup$ ... (later edit) the ACM code is written in FORTRAN 90 but my first impression is that one could call it from C without putting in a lot of effort. $\endgroup$ – GoHokies Jan 28 '17 at 8:50
  • 1
    $\begingroup$ @GoHokies I think you should convert your comment into an answer because I think it is a good answer to this question. Especially since the linked paper manages to avoid the usual numerical instabilities, and that's absolutely not a trivial thing to do. $\endgroup$ – Kirill Jan 28 '17 at 21:29
20
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – NoChance Sep 4 '12 at 21:46
  • $\begingroup$ @EmmadKareem I've updated my answer above. $\endgroup$ – GertVdE Sep 5 '12 at 8:01
3
$\begingroup$

See these:

| cite | improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ 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) $\endgroup$ – GertVdE Sep 6 '12 at 11:20
1
$\begingroup$

Numerical recipes in c provides closed form expression for real roots of quadratic and cubic which presumably have decent precision. Since the algebraic solution of the quartic involves solving a cubic and then solving two quadratics maybe a closed form quartic w good precision isn't out of the question.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I just got the root of the cubic example cited to within 2e-16 (a tad over the precision of my floats) using the numerical recipes in c (press et al) cubic formulas. So there is reason to hope. $\endgroup$ – Nemocopperfield Jan 28 '17 at 1:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.