| @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 @code{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 |
| |
| @deftypefun double gsl_poly_eval (const double @var{c}[], const int @var{len}, const double @var{x}) |
| This function evaluates the polynomial |
| @c{$c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1}$} |
| @math{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. |
| @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. |
| |
| @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{xa}, @var{ya}) 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}. |
| @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}. |
| @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 |
| |
| @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. |
| |
| @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}) |
| 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.out |
| @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 @asis |
| @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 Abramowitz and Stegun, |
| |
| @itemize @asis |
| @item |
| Abramowitz and Stegun, @cite{Handbook of Mathematical Functions}, |
| Sections 25.1.4 and 25.2.26. |
| @end itemize |