C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
|
On this page, we consider the evaluation of a polynomial function of a single variable. We usually compute the value of an arithmetic function by replacing each arithmetic operation by its corresponding floating-point machine operation. Roundoff errors and cancellations sometimes cause the calculated result to be drastically wrong. For similar reasons, a naive interval evaluation of a polynomial may lead to intervals so large as to be practically useless. Roundoff and cancellation errors are especially dangerous if we are evaluating a function close to a root, as we will see when we compute verified enclosures of zeroes of polynomials.
We present an algorithm to evaluate a real polynomial defined as
We assume the coefficients to be representable in the floating-point number system of the host computer. The algorithm achieves maximum accuracy, even in the neighborhood of a root where cancellation dooms an ordinary floating-point evaluation.
We present the algorithm RPolyEval for the evaluation of a real polynomial with maximum accuracy. Except for the special cases
and
, which can be calculated directly, an iterative solution method is used. We first compute a floating-point approximation of
. We then carry out a residual iteration by solving a linear system of equations. The new solution interval determined in the next step is checked for being of maximum accuracy, i.e. for being exact to one unit in the last place of the mantissa (1 ulp).
We list the C++ program code for the evaluation of a real polynomial with maximum accuracy. Interval variables are named with double characters, e.g. denotes the interval
.
The header file of module rpoly supplies the definition of the class RPolynomial representing a real polynomial . Since arithmetic operations for real polynomials are not required by the algorithm, they are not provided by module rpoly.
The input operator >> echos the indices on the standard output device to give an information on the index of the coefficient actually to be entered.
The header file of module rpeval supplies the interface of the function RPolyEval(), for evaluating a real polynomial with maximum accuracy and the interface of the function RPolyEvalErrMsg() which can be used to get an error message for the error code returned by RPolyEval().
The function RPolyEval() is an implementation of the algorithm. It computes an approximation and an enclosure of the value of the real polynomial with
at a point
.
We consider the polynomial to be evaluated in the neighborhood of the real value
and the polynomial
to be evaluated in the neighborhood of the real value
. To make sure that the arguments are representable on the computer, we use the machine numbers
and
, respectively.
To illustrate the difficulties that may occur with the calculation of a polynomial value when using floating-point arithmetic, these two polynomials have been evaluated in floating-point arithmetic for 100 values in the neighborhood of their roots. The corresponding plots are given in the above two figures. These two examples are solved reliably using the following sample program to call RPolyEval().
Our implementation of the algorithm produces the output listed below. In both cases the floating-point approximation, naively using Horner's nested multiplication form, does not lie within the verified enclosure of and
. An even worse side effect of rounding errors occuring during the evaluation using Horner's scheme is that the standard algorithm returns not only wrong digits but also an incorrect sign for the values of
and
. To avoid an incorrect interpretation of the resulting interval, we stress that this interval is numerically proved to be an enclosure of the exact value of the given polynomials for the machine numbers
and
.