Roots of a quadratic equation

[ivan-pi]

Motivation

Equations such as $ax^2 + bx + c = 0$ are common across all fields of science and engineering. Often one would like to find the roots of such a quadratic equation. The solution of the quadratic is commonly given as

$$ x_{1,2} = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a} $$

The formula above suffers from issues such as cancellation errors in the determinant calculation and indeterminacy when $a = 0$. Writing a robust solver which takes into account such corner cases is not a trivial feat. More details can be found under the following pages:

• Numerically stable method for solving quadratic equations | Stack Exchange
• Numerically Stable Method for Solving Quadratic Equations (PDF, 28 KB)
• Numerically stable algorithm for solving the quadratic equation when 𝑎 is very small or 0 | Math StackExchange
• Quadratic equation - Avoiding Loss of significance | Wikipedia
• Nonweiler, T. R. (1968). Algorithm 326: Roots of low-order polynomial equations. Communications of the ACM, 11(4), 269-270. (PDF, 481 KB)
• On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic | W. Kahan, 2004 (see routine qdrtc)
• HOW DO YOU SOLVE A QUADRATIC EQUATION? BY GEORGE E. FORSYTHE, Technical Report No. CS40, 1966

Prior Art

MATLAB: roots
Boost: quadratic_roots
Julia: Polynomials.roots
NAG: quadratic_​real and quadratic_​complex

Additional Information

One possible interface might look as follows:

interface quadratic_roots
  subroutine squadratic_roots(a,b,c,x0,x1)
    import sp
    real(sp), intent(in) :: a, b, c
    real(sp), intent(out) :: x0, x1
  end subroutine
  subroutine dquadratic_roots(a,b,c,x0,x1)
    import dp
    real(dp), intent(in) :: a, b, c
    real(dp), intent(out) :: x0, x1
  end subroutine
  subroutine cquadratic_roots(a,b,c,x0,x1)
    import sp
    complex(sp), intent(in) :: a, b, c
    complex(sp), intent(out) :: x0, x1
  end subroutine
  subroutine zquadratic_roots(a,b,c,x0,x1)
    import dp
    complex(dp), intent(in) :: a, b, c
    complex(dp), intent(out) :: x0, x1
  end subroutine
end interface

If bind(C) were allowed, an implementation based upon Boost C++ libraries could look as follows:

#include <boost/math/tools/roots.hpp>

#define NAME(name,kind) kind ## name

#define QUADRATIC_ROOTS(kind,T)                                  \
void NAME(quadratic_roots,kind)                                  \
  (const T* a, const T* b, const T* c, T* x0, T* x1) {           \
  auto [r0, r1] = boost::math::tools::quadratic_roots(*a,*b,*c); \
  *x0 = r0;                                                      \
  *x1 = r1;                                                      \
};

#define FOR_EACH_DATA_TYPE(X) \
  X(s,float)                  \
  X(d,double)                 \
  X(c,std::complex<float>)    \
  X(z,std::complex<double>)         

extern "C" {
FOR_EACH_DATA_TYPE(QUADRATIC_ROOTS)
}

[arjenmarkus]

In the two cases of real arguments x0 and x1, how do you communicate that the determinant is negative, so that the roots are complex? What do you do with a coefficient "a" equal to zero? (Oh and just a nitpick: why x0 and x1 and not x1 and x2?)

[ivan-pi]

I arrived at the interface above as a wrapper of the routine quadratic_roots in the Boost C++ library, with the following behaviour:

To solve a quadratic ax2 + bx + c = 0, we may use

auto [x0, x1] = boost::math::tools::quadratic_roots(a, b, c);

If the roots are real, they are arranged so that x0 ≤ x1. If the roots are complex and the inputs are real, x0 and x1 are both std::numeric_limits::quiet_NaN(). In this case we must cast a, b and c to a complex type to extract the complex roots. If a, b and c are integral, then the roots are of type double.

If a == 0, Boost does the following:

• b == 0 and c /= 0: no roots, return x1 = x2 = NaN
• b == 0 and c == 0: the x-axis (any x would solve), return x1 = x2 = 0
• b /= 0: linear equation, return x1 = x2 = -c/b

The other possibility is to return an integer error flag like the NAG routine quadratic_​real does:

Subroutine c02ajf (a,b,c,zsm,zlg,ifail)
Integer, Intent (Inout)	         :: ifail
Real (Kind=nag_wp), Intent (In)	 :: a, b, c
Real (Kind=nag_wp), Intent (Out) :: zsm(2), zlg(2)

In IMSL on the other hand, they just provide a general routine for zeros of polynomials:

      INTEGER    NDEG
      PARAMETER  (NDEG=3)
!
      REAL       COEFF(NDEG+1)
      COMPLEX    ZERO(NDEG)
! ...                               Set values of COEFF
!
      CALL ZPORC (COEFF, ZERO)

A second alternative would be to return the solution as a derived type composed of two complex numbers, but there is no elegant way to do this without PDTs for kind genericism. A third approach might be to expect a Monic polynomial of the form x**2 + p*x + q = 0 such that two solutions are guaranteed to exist.

[Romendakil]

I like the general idea, however want to point out that especially for numerical testing purposes depending on a library like Boost could be problematic. I know of (C++) codes in my field that have eliminated the Boost dependency for that reason.

[ivan-pi]

We don't need to use Boost; I just used it as a stepping-stone to arrive at an interface quickly. We are free to tweak the procedure prototype and behaviour to our liking.

[ivan-pi]

Just wanted to mention two more quadratic equation solvers:

• SLICOT MC01VD
• NSWC Library of Subroutines dqdcrt

The interfaces are

      SUBROUTINE MC01VD( A, B, C, Z1RE, Z1IM, Z2RE, Z2IM, INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO
      DOUBLE PRECISION  A, B, C, Z1IM, Z1RE, Z2IM, Z2RE

   pure subroutine dqdcrt(a, zr, zi)
      real(wp), dimension(3), intent(in)    :: a    !! coefficients
      real(wp), dimension(2), intent(out)   :: zr   !! real components of roots
      real(wp), dimension(2), intent(out)   :: zi   !! imaginary components of roots

The main point of contention is how to retrieve the roots (i.e. as scalars, arrays, or complex numbers). A user can always write a small wrapper if he prefers one approach over the provided one.

[danieljprice]

possibly helpful, here's my implementation for a cubic solver which also handles quadratics: there are two routines, one that returns the real roots only (cubicsolve, returns x(3) and also an integer for the number of real roots) and one that returns all 3 roots as complex numbers:

https://github.com/danieljprice/phantom/blob/master/src/utils/cubicsolve.f90

Would suggest at least that the API should return an array of (complex or real) roots rather than x0,x1. Then can easily be generalised to cubic, quartic etc.