Roots of Polynomials

Classes are provided for solving quadratic, cubic, and quartic equations as well as general polynomials. There is a standard nomenclature: classes which handle polynomials with real coefficients and real roots end with the suffix _real (o2scl::quadratic_real, o2scl::cubic_real and o2scl::quartic_real), classes which handle real coefficients and complex roots end with the suffix _real_coeff (o2scl::quadratic_real_coeff, o2scl::cubic_real_coeff, o2scl::quartic_real_coeff, and o2scl::poly_real_coeff), and classes which handle complex polynomials with complex coefficients end with the suffix _complex (o2scl::quadratic_complex, o2scl::cubic_complex, o2scl::quartic_complex, and o2scl::poly_complex). As a reminder, complex roots may not occur in conjugate pairs if the coefficients are not real. Most of these routines will return an error if the leading coefficient is zero.

In the public interfaces to the polynomial solvers, the complex type std::complex<double> is used.

For quadratics, o2scl::quadratic_real_coeff_gsl is the best if the coefficients are real, while if the coefficients are complex, use o2scl::quadratic_complex_std. For cubics with real coefficients, o2scl::cubic_real_coeff_cern is the best, while if the coefficients are complex, use o2scl::cubic_complex_std.

For a quartic polynomial with real coefficients, o2scl::quartic_real_coeff_cern is the best, unless the coefficients of odd powers happen to be small, in which case, o2scl::quartic_real_gsl2 tends to work better. For quartics, generic polynomial solvers such as o2scl::poly_real_coeff_gsl can provide more accurate (but slower) results. If the coefficients are complex, then you can use o2scl::quartic_complex_simple.

Polynomial solver example

This example shows how to find the roots of the second-, third-, fourth-, and fifth-order Chebyshev polynomials

\begin{eqnarray*} &2x^2-1& \nonumber \\ &4x^3-3 x& \nonumber \\ &8x^4-8x^2+1& \nonumber \\ &16x^5-20x^3+5x& \nonumber \end{eqnarray*}

For the Chebyshev polynomial of order $ n $, the roots are given by

\[ \cos \left[ \frac{\pi(k-1/2)}{n}\right] \]

for $ k = 1,\ldots,n $ These roots are used in o2scl::cheb_approx to approximate functions using Chebyshev polynomials .

/* Example: ex_poly.cpp
Demonstrate the solution of the Chebyshev polynomials
#include <boost/numeric/ublas/vector.hpp>
#include <o2scl/poly.h>
// For pi
#include <o2scl/constants.h>
#include <o2scl/vector.h>
#include <o2scl/test_mgr.h>
using namespace std;
using namespace o2scl;
using namespace o2scl_const;
int main(void) {
// Quadratic solver
// Cubic solver
// Quartic solver
// Generic polynomial solver
// Storage for the roots
ubvector v(5);
double d;
std::complex<double> ca[5];
// The second order polynomial
cout << "Second order roots: " << endl;
// Sort the roots and compare with the exact results
for(size_t i=0;i<2;i++) {
double exact=cos(pi*(((double)(2-i))-0.5)/2.0);
cout << v[i] << " " << exact << endl;
t.test_abs(v[i],exact,1.0e-14,"2nd order");
cout << endl;
// The third order polynomial
cout << "Third order roots: " << endl;
// Sort the roots and compare with the exact results
for(size_t i=0;i<3;i++) {
double exact=cos(pi*(((double)(3-i))-0.5)/3.0);
cout << v[i] << " " << exact << endl;
if (i==1) {
t.test_abs(v[i],exact,1.0e-14,"3rd order");
} else {
t.test_abs(v[i],exact,1.0e-14,"3rd order");
cout << endl;
// The fourth order polynomial
cout << "Fourth order roots: " << endl;
// Sort the roots and compare with the exact results
for(size_t i=0;i<4;i++) v[i]=ca[i].real();
for(size_t i=0;i<4;i++) {
double exact=cos(pi*(((double)(4-i))-0.5)/4.0);
cout << v[i] << " " << exact << endl;
t.test_abs(v[i],exact,1.0e-14,"4th order");
cout << endl;
// The fifth order polynomial
cout << "Fifth order roots: " << endl;
double co[6]={16.0,0.0,-20.0,0.0,5.0,0.0};
// Sort the roots and compare with the exact results
for(size_t i=0;i<5;i++) v[i]=ca[i].real();
for(size_t i=0;i<5;i++) {
double exact=cos(pi*(((double)(5-i))-0.5)/5.0);
cout << v[i] << " " << exact << endl;
t.test_abs(v[i],exact,1.0e-14,"5th order");
cout << endl;;
return 0;
// End of example

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).