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.
Algorithmic Description
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).
The RPolyEval Algorithm
Implementation and Examples
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 .
Module rpoly
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.
1/*
2** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
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().
1/*
2** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
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 .
1/*
2** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
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.
Polynomial p(t) in interval [1.9995 , 2.0005]
Polynomial q(t) in interval [0.99999 , 1.00001]
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().
34 cout << "Enter the degree of the polynomial (>=0): ";
35 cin >> n; cout << endl;
36 } while (n < 0);
37
38 RPolynomial p(n);
39
40 cout << SetPrecision(23,15) << Scientific; // Output format
41
42 cout << "Enter the coefficients of p in increasing order:" << endl;
43 cin >> p; cout << endl;
44 cout << "Enter the argument t for evaluation:" << endl;
45 cin >> t; cout << endl;
46
47 RPolyEval(p,t,y,yy,No,Err);
48
49if (!Err) {
50 cout << "Polynomial:" << endl << p <<endl;
51 cout << "Floating-point evaluation of p(t) using Horner's scheme:"
52 << endl << " " << y << endl << endl;
53 cout << "Verified inclusion of p(t):"
54 << endl << yy << endl << endl;
55 cout << "Number of iterations needed: " << No << endl;
56 }
57else
58 cout << RPolyEvalErrMsg(Err) << endl;
59
60return 0;
61}
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 .