htswanalysis/MACS/lib/gsl/gsl-1.11/doc/poly.texi

   1 @cindex polynomials, roots of
   2
   3 This chapter describes functions for evaluating and solving polynomials.
   4 There are routines for finding real and complex roots of quadratic and
   5 cubic equations using analytic methods.  An iterative polynomial solver
   6 is also available for finding the roots of general polynomials with real
   7 coefficients (of any order).  The functions are declared in the header
   8 file @code{gsl_poly.h}.
   9
  10 @menu
  11 * Polynomial Evaluation::
  12 * Divided Difference Representation of Polynomials::
  13 * Quadratic Equations::
  14 * Cubic Equations::
  15 * General Polynomial Equations::
  16 * Roots of Polynomials Examples::
  17 * Roots of Polynomials References and Further Reading::
  18 @end menu
  19
  20 @node Polynomial Evaluation
  21 @section Polynomial Evaluation
  22 @cindex polynomial evaluation
  23 @cindex evaluation of polynomials
  24
  25 The functions described here evaluate the polynomial
  26 @c{$c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1}$}
  27 @math{c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^@{len-1@}} using
  28 Horner's method for stability. @inlinefns{}
  29
  30 @deftypefun double gsl_poly_eval (const double @var{c}[], const int @var{len}, const double @var{x})
  31 This function evaluates a polynomial with real coefficients for the real variable @var{x}.
  32 @end deftypefun
  33
  34 @deftypefun gsl_complex gsl_poly_complex_eval (const double @var{c}[], const int @var{len}, const gsl_complex @var{z})
  35 This function evaluates a polynomial with real coefficients for the complex variable @var{z}.
  36 @end deftypefun
  37
  38 @deftypefun gsl_complex gsl_complex_poly_complex_eval (const gsl_complex @var{c}[], const int @var{len}, const gsl_complex @var{z})
  39 This function evaluates a polynomial with complex coefficients for the complex variable @var{z}.
  40 @end deftypefun
  41
  42 @node Divided Difference Representation of Polynomials
  43 @section Divided Difference Representation of Polynomials
  44 @cindex divided differences, polynomials
  45 @cindex evaluation of polynomials, in divided difference form
  46
  47 The functions described here manipulate polynomials stored in Newton's
  48 divided-difference representation.  The use of divided-differences is
  49 described in Abramowitz & Stegun sections 25.1.4 and 25.2.26.
  50
  51 @deftypefun int gsl_poly_dd_init (double @var{dd}[], const double @var{xa}[], const double @var{ya}[], size_t @var{size})
  52 This function computes a divided-difference representation of the
  53 interpolating polynomial for the points (@var{xa}, @var{ya}) stored in
  54 the arrays @var{xa} and @var{ya} of length @var{size}.  On output the
  55 divided-differences of (@var{xa},@var{ya}) are stored in the array
  56 @var{dd}, also of length @var{size}.
  57 @end deftypefun
  58
  59 @deftypefun double gsl_poly_dd_eval (const double @var{dd}[], const double @var{xa}[], const size_t @var{size}, const double @var{x})
  60 This function evaluates the polynomial stored in divided-difference form
  61 in the arrays @var{dd} and @var{xa} of length @var{size} at the point
  62 @var{x}.  @inlinefn{}
  63 @end deftypefun
  64
  65 @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}[])
  66 This function converts the divided-difference representation of a
  67 polynomial to a Taylor expansion.  The divided-difference representation
  68 is supplied in the arrays @var{dd} and @var{xa} of length @var{size}.
  69 On output the Taylor coefficients of the polynomial expanded about the
  70 point @var{xp} are stored in the array @var{c} also of length
  71 @var{size}.  A workspace of length @var{size} must be provided in the
  72 array @var{w}.
  73 @end deftypefun
  74
  75 @node  Quadratic Equations
  76 @section Quadratic Equations
  77 @cindex quadratic equation, solving
  78
  79 @deftypefun int gsl_poly_solve_quadratic (double @var{a}, double @var{b}, double @var{c}, double * @var{x0}, double * @var{x1})
  80 This function finds the real roots of the quadratic equation,
  81 @tex
  82 \beforedisplay
  83 $$
  84 a x^2 + b x + c = 0
  85 $$
  86 \afterdisplay
  87 @end tex
  88 @ifinfo
  89
  90 @example
  91 a x^2 + b x + c = 0
  92 @end example
  93
  94 @end ifinfo
  95 @noindent
  96 The number of real roots (either zero, one or two) is returned, and
  97 their locations are stored in @var{x0} and @var{x1}.  If no real roots
  98 are found then @var{x0} and @var{x1} are not modified.  If one real root
  99 is found (i.e. if @math{a=0}) then it is stored in @var{x0}.  When two
 100 real roots are found they are stored in @var{x0} and @var{x1} in
 101 ascending order.  The case of coincident roots is not considered
 102 special.  For example @math{(x-1)^2=0} will have two roots, which happen
 103 to have exactly equal values.
 104
 105 The number of roots found depends on the sign of the discriminant
 106 @math{b^2 - 4 a c}.  This will be subject to rounding and cancellation
 107 errors when computed in double precision, and will also be subject to
 108 errors if the coefficients of the polynomial are inexact.  These errors
 109 may cause a discrete change in the number of roots.  However, for
 110 polynomials with small integer coefficients the discriminant can always
 111 be computed exactly.
 112
 113 @end deftypefun
 114
 115 @deftypefun int gsl_poly_complex_solve_quadratic (double @var{a}, double @var{b}, double @var{c}, gsl_complex * @var{z0}, gsl_complex * @var{z1})
 116
 117 This function finds the complex roots of the quadratic equation,
 118 @tex
 119 \beforedisplay
 120 $$
 121 a z^2 + b z + c = 0
 122 $$
 123 \afterdisplay
 124 @end tex
 125 @ifinfo
 126
 127 @example
 128 a z^2 + b z + c = 0
 129 @end example
 130
 131 @end ifinfo
 132 @noindent
 133 The number of complex roots is returned (either one or two) and the
 134 locations of the roots are stored in @var{z0} and @var{z1}.  The roots
 135 are returned in ascending order, sorted first by their real components
 136 and then by their imaginary components.  If only one real root is found
 137 (i.e. if @math{a=0}) then it is stored in @var{z0}.
 138
 139 @end deftypefun
 140
 141
 142 @node Cubic Equations
 143 @section Cubic Equations
 144 @cindex cubic equation, solving
 145
 146 @deftypefun int gsl_poly_solve_cubic (double @var{a}, double @var{b}, double @var{c}, double * @var{x0}, double * @var{x1}, double * @var{x2})
 147
 148 This function finds the real roots of the cubic equation,
 149 @tex
 150 \beforedisplay
 151 $$
 152 x^3 + a x^2 + b x + c = 0
 153 $$
 154 \afterdisplay
 155 @end tex
 156 @ifinfo
 157
 158 @example
 159 x^3 + a x^2 + b x + c = 0
 160 @end example
 161
 162 @end ifinfo
 163 @noindent
 164 with a leading coefficient of unity.  The number of real roots (either
 165 one or three) is returned, and their locations are stored in @var{x0},
 166 @var{x1} and @var{x2}.  If one real root is found then only @var{x0}
 167 is modified.  When three real roots are found they are stored in
 168 @var{x0}, @var{x1} and @var{x2} in ascending order.  The case of
 169 coincident roots is not considered special.  For example, the equation
 170 @math{(x-1)^3=0} will have three roots with exactly equal values.  As
 171 in the quadratic case, finite precision may cause equal or
 172 closely-spaced real roots to move off the real axis into the complex
 173 plane, leading to a discrete change in the number of real roots.
 174 @end deftypefun
 175
 176 @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})
 177
 178 This function finds the complex roots of the cubic equation,
 179 @tex
 180 \beforedisplay
 181 $$
 182 z^3 + a z^2 + b z + c = 0
 183 $$
 184 \afterdisplay
 185 @end tex
 186 @ifinfo
 187
 188 @example
 189 z^3 + a z^2 + b z + c = 0
 190 @end example
 191
 192 @end ifinfo
 193 @noindent
 194 The number of complex roots is returned (always three) and the locations
 195 of the roots are stored in @var{z0}, @var{z1} and @var{z2}.  The roots
 196 are returned in ascending order, sorted first by their real components
 197 and then by their imaginary components.
 198
 199 @end deftypefun
 200
 201
 202 @node General Polynomial Equations
 203 @section General Polynomial Equations
 204 @cindex general polynomial equations, solving
 205
 206 The roots of polynomial equations cannot be found analytically beyond
 207 the special cases of the quadratic, cubic and quartic equation.  The
 208 algorithm described in this section uses an iterative method to find the
 209 approximate locations of roots of higher order polynomials.
 210
 211 @deftypefun {gsl_poly_complex_workspace *} gsl_poly_complex_workspace_alloc (size_t @var{n})
 212 This function allocates space for a @code{gsl_poly_complex_workspace}
 213 struct and a workspace suitable for solving a polynomial with @var{n}
 214 coefficients using the routine @code{gsl_poly_complex_solve}.
 215
 216 The function returns a pointer to the newly allocated
 217 @code{gsl_poly_complex_workspace} if no errors were detected, and a null
 218 pointer in the case of error.
 219 @end deftypefun
 220
 221 @deftypefun void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * @var{w})
 222 This function frees all the memory associated with the workspace
 223 @var{w}.
 224 @end deftypefun
 225
 226 @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})
 227 This function computes the roots of the general polynomial
 228 @c{$P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1}$}
 229 @math{P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_@{n-1@} x^@{n-1@}} using
 230 balanced-QR reduction of the companion matrix.  The parameter @var{n}
 231 specifies the length of the coefficient array.  The coefficient of the
 232 highest order term must be non-zero.  The function requires a workspace
 233 @var{w} of the appropriate size.  The @math{n-1} roots are returned in
 234 the packed complex array @var{z} of length @math{2(n-1)}, alternating
 235 real and imaginary parts.
 236
 237 The function returns @code{GSL_SUCCESS} if all the roots are found. If
 238 the QR reduction does not converge, the error handler is invoked with
 239 an error code of @code{GSL_EFAILED}.  Note that due to finite precision,
 240 roots of higher multiplicity are returned as a cluster of simple roots
 241 with reduced accuracy.  The solution of polynomials with higher-order
 242 roots requires specialized algorithms that take the multiplicity
 243 structure into account (see e.g. Z. Zeng, Algorithm 835, ACM
 244 Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp
 245 218--236).
 246 @end deftypefun
 247
 248 @node Roots of Polynomials Examples
 249 @section Examples
 250
 251 To demonstrate the use of the general polynomial solver we will take the
 252 polynomial @math{P(x) = x^5 - 1} which has the following roots,
 253 @tex
 254 \beforedisplay
 255 $$
 256 1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5}
 257 $$
 258 \afterdisplay
 259 @end tex
 260 @ifinfo
 261
 262 @example
 263 1, e^@{2\pi i /5@}, e^@{4\pi i /5@}, e^@{6\pi i /5@}, e^@{8\pi i /5@}
 264 @end example
 265
 266 @end ifinfo
 267 @noindent
 268 The following program will find these roots.
 269
 270 @example
 271 @verbatiminclude examples/polyroots.c
 272 @end example
 273
 274 @noindent
 275 The output of the program is,
 276
 277 @example
 278 $ ./a.out
 279 @verbatiminclude examples/polyroots.out
 280 @end example
 281
 282 @noindent
 283 which agrees with the analytic result, @math{z_n = \exp(2 \pi n i/5)}.
 284
 285 @node Roots of Polynomials References and Further Reading
 286 @section References and Further Reading
 287
 288 The balanced-QR method and its error analysis are described in the
 289 following papers,
 290
 291 @itemize @asis
 292 @item
 293 R.S. Martin, G. Peters and J.H. Wilkinson, ``The QR Algorithm for Real
 294 Hessenberg Matrices'', @cite{Numerische Mathematik}, 14 (1970), 219--231.
 295
 296 @item
 297 B.N. Parlett and C. Reinsch, ``Balancing a Matrix for Calculation of
 298 Eigenvalues and Eigenvectors'', @cite{Numerische Mathematik}, 13 (1969),
 299 293--304.
 300
 301 @item
 302 A. Edelman and H. Murakami, ``Polynomial roots from companion matrix
 303 eigenvalues'', @cite{Mathematics of Computation}, Vol.@: 64, No.@: 210
 304 (1995), 763--776.
 305 @end itemize
 306
 307 @noindent
 308 The formulas for divided differences are given in Abramowitz and Stegun,
 309
 310 @itemize @asis
 311 @item
 312 Abramowitz and Stegun, @cite{Handbook of Mathematical Functions},
 313 Sections 25.1.4 and 25.2.26.
 314 @end itemize