Roots of a quadratic equation

ivan-pi commented 2 years ago

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:

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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 commented 4 months ago

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.

fortran-lang / stdlib