This chapter describes functions for solving polynomial equations. There are routines for finding real and complex roots of quadratic and cubic equations using analytic methods. An iterative polynomial solver is also available for finding the roots of general real polynomials (of any order). The functions are defined in the header file gsl_poly.h.

Quadratic equations
===================

 -- Function: int gsl_poly_solve_quadratic (double A, double B, double C, double *X0, double *X1)
     This function finds the real roots of the quadratic equation,

          a x^2 + b x + c = 0

     The number of real roots (either zero or two) is returned, and their locations are stored in X0 and X1. If no real roots are found then X0 and X1 are not modified. When two real roots are found they are stored in X0 and X1 in ascending order. The case of coincident roots is not considered special. For example (x-1)^2=0 will have two roots, which happen to have exactly equal values.

     The number of roots found depends on the sign of the discriminant b^2 - 4 a c.  This will be subject to rounding and cancellation errors when computed in double precision, and will also be subject to errors if the coefficients of the polynomial are inexact. These errors may cause a discrete change in the number of roots. However, for polynomials with small integer coefficients the discriminant can always be computed exactly.

 -- Function: int gsl_poly_complex_solve_quadratic (double A, double B, double C, gsl_complex *Z0, gsl_complex *Z1)

     This function finds the complex roots of the quadratic equation,

          a z^2 + b z + c = 0

     The number of complex roots is returned (always two) and the locations of the roots are stored in Z0 and Z1. The roots are returned in ascending order, sorted first by their real components and then by their imaginary components.

Cubic equations
===============

 -- Function: int gsl_poly_solve_cubic (double A, double B, double C, double *X0, double *X1, double *X2)

     This function finds the real roots of the cubic equation,

          x^3 + a x^2 + b x + c = 0

     with a leading coefficient of unity. The number of real roots (either one or three) is returned, and their locations are stored in X0, X1 and X2. If one real root is found then only X0 is modified. When three real roots are found they are stored in X0, X1 and X2 in ascending order. The case of coincident roots is not considered special. For example, the equation (x-1)^3=0 will have three roots with exactly equal values.

 -- Function: int gsl_poly_complex_solve_cubic (double A, double B, double C, gsl_complex *Z0, gsl_complex *Z1, gsl_complex *Z2)

     This function finds the complex roots of the cubic equation,

          z^3 + a z^2 + b z + c = 0

     The number of complex roots is returned (always three) and the locations of the roots are stored in Z0, Z1 and Z2. The roots are returned in ascending order, sorted first by their real components and then by their imaginary components.

General polynomial equations
============================

The roots of polynomial equations cannot be found analytically beyond the special cases of the quadratic, cubic and quartic equation. The algorithm described in this section uses an iterative method to find the approximate locations of roots of higher order polynomials.

 -- Function: gsl_poly_complex_workspace * gsl_poly_complex_workspace_alloc (size_t N)
     This function allocates space for a gsl_poly_complex_workspace struct and a scratch area suitable for solving a polynomial with N coefficients using the routine gsl_poly_complex_solve.

     The function returns a pointer to the newly allocated gsl_poly_complex_workspace if no errors were detected, and 0 in the case of error.

 -- Function: void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * W)
     This function frees all the memory associated with the workspace W.

 -- Function: int gsl_poly_complex_solve (const double * A, size_t N, gsl_poly_complex_workspace * W, gsl_complex_packed_ptr Z)
     This function computes the roots of the general polynomial $P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_n-1 x^n-1$ P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using balanced-QR reduction of the companion matrix.  The parameter N specifies the length of the coefficient array. The coefficient of the highest order term must be non-zero. The function requires a workspace W of the appropriate size. The n-1 roots are returned in the packed complex array Z of length 2(n-1), alternating real and imaginary parts.

     The functions returns GSL_SUCCESS if all the roots are found and GSL_EFAILED if the QR reduction does not converge.

Examples
========

To demonstrate the use of the general polynomial solver we will take the polynomial P(x) = x^5 - 1 which has the following roots,

     1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5}

The following program will find these roots.

     #include <stdio.h>
     #include <gsl_poly.h>

     int
     main ()
     {
       int i;
       double a[6] = { -1, 0, 0, 0, 0, 1 };  /*  P(x) =  x^5 - 1  */
       double z[10];

       gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (6) ;

       gsl_poly_complex_solve (a, 6, w, z) ;

       gsl_poly_complex_workspace_free (w) ;

       for (i = 0; i < 5 ; i++)
         {
           printf("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]) ;
         }
     }

The results of running the program are,

     bash$ ./a.out
     z0 = -0.809016994374947451 +0.587785252292473137
     z1 = -0.809016994374947451 -0.587785252292473137
     z2 = +0.309016994374947451 +0.951056516295153642
     z3 = +0.309016994374947451 -0.951056516295153642
     z4 = +1.000000000000000000 +0.000000000000000000

References and Further Reading
==============================

The balanced-QR method and its error analysis is described in the following papers.

  R.S. Martin, G. Peters and J.H. Wilkinson, ``The QR Algorithm for Real Hessenberg Matrices'', Numerische Mathematik, 14 (1970), 219--231.

  B.N. Parlett and C. Reinsch, ``Balancing a Matrix for Calculation of Eigenvalues and Eigenvectors'', Numerische Mathematik, 13 (1969), 293--304.

  A. Edelman and H. Murakami, ``Polynomial roots from companion matrix eigenvalues'', Mathematics of Computation, Vol. 64 No. 210 (1995), 763--776.