@cindex polynomials, roots of 

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 @file{gsl_poly.h}.

@menu
* Polynomial Evaluation::       
* Divided Difference Representation of Polynomials::  
* Quadratic Equations::         
* Cubic Equations::             
* General Polynomial Equations::  
* Roots of Polynomials Examples::  
* Roots of Polynomials References and Further Reading::  
@end menu

@node Polynomial Evaluation
@section Polynomial Evaluation
@cindex polynomial evaluation
@cindex evaluation of polynomials

The functions described here evaluate the polynomial 
@c{$P(x) = c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1}$}
@math{P(x) = c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^@{len-1@}} using
Horner's method for stability. @inlinefns{}

@deftypefun double gsl_poly_eval (const double @var{c}[], const int @var{len}, const double @var{x})
This function evaluates a polynomial with real coefficients for the real variable @var{x}.
@end deftypefun

@deftypefun gsl_complex gsl_poly_complex_eval (const double @var{c}[], const int @var{len}, const gsl_complex @var{z})
This function evaluates a polynomial with real coefficients for the complex variable @var{z}.
@end deftypefun

@deftypefun gsl_complex gsl_complex_poly_complex_eval (const gsl_complex @var{c}[], const int @var{len}, const gsl_complex @var{z})
This function evaluates a polynomial with complex coefficients for the complex variable @var{z}.
@end deftypefun

@deftypefun int gsl_poly_eval_derivs (const double @var{c}[], const size_t @var{lenc}, const double @var{x}, double @var{res}[], const size_t @var{lenres})
This function evaluates a polynomial and its derivatives storing the
results in the array @var{res} of size @var{lenres}.  The output array
contains the values of @math{d^k P/d x^k} for the specified value of
@var{x} starting with @math{k = 0}.
@end deftypefun

@node Divided Difference Representation of Polynomials
@section Divided Difference Representation of Polynomials
@cindex divided differences, polynomials
@cindex evaluation of polynomials, in divided difference form

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 and 25.2.26, and
Burden and Faires, chapter 3, and discussed briefly below.

@noindent
Given a function @math{f(x)}, an @math{n}th degree interpolating polynomial @math{P_{n}(x)}
can be constructed which agrees with @math{f} at @math{n+1} distinct points
@math{x_0,x_1,...,x_{n}}. This polynomial can be written in a
form known as Newton's divided-difference representation:
@tex
\beforedisplay
$$
P_{n}(x) = f(x_0) + \sum_{k=1}^n [x_0,x_1,...,x_k] (x-x_0)(x-x_1) \cdots (x-x_{k-1})
$$
\afterdisplay
@end tex
@ifinfo

@example
P_n(x) = f(x_0) + \sum_(k=1)^n [x_0,x_1,...,x_k] (x-x_0)(x-x_1)...(x-x_(k-1))
@end example

@end ifinfo
where the divided differences @math{[x_0,x_1,...,x_k]} are defined in section 25.1.4 of
Abramowitz and Stegun. Additionally, it is possible to construct an interpolating
polynomial of degree @math{2n+1} which also matches the first derivatives of @math{f}
at the points @math{x_0,x_1,...,x_n}. This is called the Hermite interpolating
polynomial and is defined as
@tex
\beforedisplay
$$
H_{2n+1}(x) = f(z_0) + \sum_{k=1}^{2n+1} [z_0,z_1,...,z_k] (x-z_0)(x-z_1) \cdots (x-z_{k-1})
$$
\afterdisplay
@end tex
@ifinfo

@example
H_(2n+1)(x) = f(z_0) + \sum_(k=1)^(2n+1) [z_0,z_1,...,z_k] (x-z_0)(x-z_1)...(x-z_(k-1))
@end example

@end ifinfo
where the elements of @math{z = \{x_0,x_0,x_1,x_1,...,x_n,x_n\}} are defined by
@math{z_{2k} = z_{2k+1} = x_k}. The divided-differences @math{[z_0,z_1,...,z_k]}
are discussed in Burden and Faires, section 3.4.

@deftypefun int gsl_poly_dd_init (double @var{dd}[], const double @var{xa}[], const double @var{ya}[], size_t @var{size})
This function computes a divided-difference representation of the
interpolating polynomial for the points (@var{x}, @var{y}) stored in
the arrays @var{xa} and @var{ya} of length @var{size}.  On output the
divided-differences of (@var{xa},@var{ya}) are stored in the array
@var{dd}, also of length @var{size}. Using the notation above,
@math{dd[k] = [x_0,x_1,...,x_k]}.
@end deftypefun

@deftypefun double gsl_poly_dd_eval (const double @var{dd}[], const double @var{xa}[], const size_t @var{size}, const double @var{x})
This function evaluates the polynomial stored in divided-difference form
in the arrays @var{dd} and @var{xa} of length @var{size} at the point
@var{x}.  @inlinefn{}
@end deftypefun

@deftypefun int gsl_poly_dd_taylor (double @var{c}[], double @var{xp}, const double @var{dd}[], const double @var{xa}[], size_t @var{size}, double @var{w}[])
This function converts the divided-difference representation of a
polynomial to a Taylor expansion.  The divided-difference representation
is supplied in the arrays @var{dd} and @var{xa} of length @var{size}.
On output the Taylor coefficients of the polynomial expanded about the
point @var{xp} are stored in the array @var{c} also of length
@var{size}.  A workspace of length @var{size} must be provided in the
array @var{w}.
@end deftypefun

@deftypefun int gsl_poly_dd_hermite_init (double @var{dd}[], double @var{za}[], const double @var{xa}[], const double @var{ya}[], const double @var{dya}[], const size_t @var{size})
This function computes a divided-difference representation of the
interpolating Hermite polynomial for the points (@var{x}, @var{y}) stored in
the arrays @var{xa} and @var{ya} of length @var{size}. Hermite interpolation
constructs polynomials which also match first derivatives @math{dy/dx} which are
provided in the array @var{dya} also of length @var{size}. The first derivatives can be
incorported into the usual divided-difference algorithm by forming a new
dataset @math{z = \{x_0,x_0,x_1,x_1,...\}}, which is stored in the array
@var{za} of length 2*@var{size} on output. On output the
divided-differences of the Hermite representation are stored in the array
@var{dd}, also of length 2*@var{size}. Using the notation above,
@math{dd[k] = [z_0,z_1,...,z_k]}. The resulting Hermite polynomial
can be evaluated by calling @code{gsl_poly_dd_eval} and using
@var{za} for the input argument @var{xa}.
@end deftypefun

@node  Quadratic Equations
@section Quadratic Equations
@cindex quadratic equation, solving

@deftypefun int gsl_poly_solve_quadratic (double @var{a}, double @var{b}, double @var{c}, double * @var{x0}, double * @var{x1})
This function finds the real roots of the quadratic equation,
@tex
\beforedisplay
$$
a x^2 + b x + c = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
a x^2 + b x + c = 0
@end example

@end ifinfo
@noindent
The number of real roots (either zero, one or two) is returned, and
their locations are stored in @var{x0} and @var{x1}.  If no real roots
are found then @var{x0} and @var{x1} are not modified.  If one real root
is found (i.e. if @math{a=0}) then it is stored in @var{x0}.  When two
real roots are found they are stored in @var{x0} and @var{x1} in
ascending order.  The case of coincident roots is not considered
special.  For example @math{(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
@math{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.

@end deftypefun

@deftypefun int gsl_poly_complex_solve_quadratic (double @var{a}, double @var{b}, double @var{c}, gsl_complex * @var{z0}, gsl_complex * @var{z1})

This function finds the complex roots of the quadratic equation,
@tex
\beforedisplay
$$
a z^2 + b z + c = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
a z^2 + b z + c = 0
@end example

@end ifinfo
@noindent
The number of complex roots is returned (either one or two) and the
locations of the roots are stored in @var{z0} and @var{z1}.  The roots
are returned in ascending order, sorted first by their real components
and then by their imaginary components.  If only one real root is found
(i.e. if @math{a=0}) then it is stored in @var{z0}.

@end deftypefun


@node Cubic Equations
@section Cubic Equations
@cindex cubic equation, solving

@deftypefun int gsl_poly_solve_cubic (double @var{a}, double @var{b}, double @var{c}, double * @var{x0}, double * @var{x1}, double * @var{x2})

This function finds the real roots of the cubic equation,
@tex
\beforedisplay
$$
x^3 + a x^2 + b x + c = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
x^3 + a x^2 + b x + c = 0
@end example

@end ifinfo
@noindent
with a leading coefficient of unity.  The number of real roots (either
one or three) is returned, and their locations are stored in @var{x0},
@var{x1} and @var{x2}.  If one real root is found then only @var{x0}
is modified.  When three real roots are found they are stored in
@var{x0}, @var{x1} and @var{x2} in ascending order.  The case of
coincident roots is not considered special.  For example, the equation
@math{(x-1)^3=0} will have three roots with exactly equal values.  As
in the quadratic case, finite precision may cause equal or
closely-spaced real roots to move off the real axis into the complex
plane, leading to a discrete change in the number of real roots.
@end deftypefun

@deftypefun int gsl_poly_complex_solve_cubic (double @var{a}, double @var{b}, double @var{c}, gsl_complex * @var{z0}, gsl_complex * @var{z1}, gsl_complex * @var{z2})

This function finds the complex roots of the cubic equation,
@tex
\beforedisplay
$$
z^3 + a z^2 + b z + c = 0
$$
\afterdisplay
@end tex
@ifinfo

@example
z^3 + a z^2 + b z + c = 0
@end example

@end ifinfo
@noindent
The number of complex roots is returned (always three) and the locations
of the roots are stored in @var{z0}, @var{z1} and @var{z2}.  The roots
are returned in ascending order, sorted first by their real components
and then by their imaginary components.

@end deftypefun


@node General Polynomial Equations
@section General Polynomial Equations
@cindex general polynomial equations, solving

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.

@deftypefun {gsl_poly_complex_workspace *} gsl_poly_complex_workspace_alloc (size_t @var{n})
@tindex gsl_poly_complex_workspace
This function allocates space for a @code{gsl_poly_complex_workspace}
struct and a workspace suitable for solving a polynomial with @var{n}
coefficients using the routine @code{gsl_poly_complex_solve}.

The function returns a pointer to the newly allocated
@code{gsl_poly_complex_workspace} if no errors were detected, and a null
pointer in the case of error.
@end deftypefun

@deftypefun void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * @var{w})
This function frees all the memory associated with the workspace
@var{w}.
@end deftypefun

@deftypefun int gsl_poly_complex_solve (const double * @var{a}, size_t @var{n}, gsl_poly_complex_workspace * @var{w}, gsl_complex_packed_ptr @var{z})
This function computes the roots of the general polynomial 
@c{$P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1}$} 
@math{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 @var{n}
specifies the length of the coefficient array.  The coefficient of the
highest order term must be non-zero.  The function requires a workspace
@var{w} of the appropriate size.  The @math{n-1} roots are returned in
the packed complex array @var{z} of length @math{2(n-1)}, alternating
real and imaginary parts.

The function returns @code{GSL_SUCCESS} if all the roots are found. If
the QR reduction does not converge, the error handler is invoked with
an error code of @code{GSL_EFAILED}.  Note that due to finite precision,
roots of higher multiplicity are returned as a cluster of simple roots
with reduced accuracy.  The solution of polynomials with higher-order
roots requires specialized algorithms that take the multiplicity
structure into account (see e.g. Z. Zeng, Algorithm 835, ACM
Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp
218--236).
@end deftypefun

@node Roots of Polynomials Examples
@section Examples

To demonstrate the use of the general polynomial solver we will take the
polynomial @math{P(x) = x^5 - 1} which has the following roots,
@tex
\beforedisplay
$$
1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5}
$$
\afterdisplay
@end tex
@ifinfo

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

@end ifinfo
@noindent
The following program will find these roots.

@example
@verbatiminclude examples/polyroots.c
@end example

@noindent
The output of the program is,

@example
$ ./a.out 
@verbatiminclude examples/polyroots.txt
@end example

@noindent
which agrees with the analytic result, @math{z_n = \exp(2 \pi n i/5)}.

@node Roots of Polynomials References and Further Reading
@section References and Further Reading

The balanced-QR method and its error analysis are described in the
following papers,

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

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

@item
A. Edelman and H. Murakami, ``Polynomial roots from companion matrix
eigenvalues'', @cite{Mathematics of Computation}, Vol.@: 64, No.@: 210
(1995), 763--776.
@end itemize

@noindent
The formulas for divided differences are given in the following texts,

@itemize @w{}
@item
Abramowitz and Stegun, @cite{Handbook of Mathematical Functions},
Sections 25.1.4 and 25.2.26.

@item
R. L. Burden and J. D. Faires, @cite{Numerical Analysis}, 9th edition,
ISBN 0-538-73351-9, 2011.
@end itemize