This chapter describes functions for evaluating and solving polynomials.
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 polynomials with real
coefficients (of any order). The functions are declared in the header
file gsl_poly.h.
This function evaluates the polynomial
c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1} using
Horner's method for stability. The function is inlined when possible.
The functions described here manipulate polynomials stored in Newton's
divided-difference representation. The use of divided-differences is
described in Abramowitz & Stegun sections 25.1.4, 25.2.26.
This function computes a divided-difference representation of the
interpolating polynomial for the points (xa, ya) stored in
the arrays xa and ya of length size. On output the
divided-differences of (xa,ya) are stored in the array
dd, also of length size.
This function converts the divided-difference representation of a
polynomial to a Taylor expansion. The divided-difference representation
is supplied in the arrays dd and xa of length size.
On output the Taylor coefficients of the polynomial expanded about the
point xp are stored in the array c also of length
size. A workspace of length size must be provided in the
array w.
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.
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.
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.
This function allocates space for a gsl_poly_complex_workspace
struct and a workspace 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 a null
pointer in the case of error.
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} 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 function returns GSL_SUCCESS if all the roots are found and
GSL_EFAILED if the QR reduction does not converge.