C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
Evaluation of Polynomials

Overview

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 $ p : R \to R $ defined as

\[
p(t) = \sum \limits_{i=0}^n p_i t^i \;,\; p_i , t \in R \;,\; i = 0 , ... , n \;,\; p_n \not= 0
\]

We assume the coefficients $ p_i $ 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 $ p(t) = \sum \limits_{i=0}^n p_i t^i $ with maximum accuracy. Except for the special cases $ n = 0 $ and $ n = 1 $, which can be calculated directly, an iterative solution method is used. We first compute a floating-point approximation of $ p(t) $. 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. $ rr[i] $ denotes the interval $ [r]_i $.

Module rpoly

The header file of module rpoly supplies the definition of the class RPolynomial representing a real polynomial $ p(t) = \sum \limits_{i=0}^n p_i t^i \;,\; p_i , t \in R $. 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)
3**
4** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5** Universitaet Karlsruhe, Germany
6** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7** Universitaet Wuppertal, Germany
8**
9** This library is free software; you can redistribute it and/or
10** modify it under the terms of the GNU Library General Public
11** License as published by the Free Software Foundation; either
12** version 2 of the License, or (at your option) any later version.
13**
14** This library is distributed in the hope that it will be useful,
15** but WITHOUT ANY WARRANTY; without even the implied warranty of
16** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17** Library General Public License for more details.
18**
19** You should have received a copy of the GNU Library General Public
20** License along with this library; if not, write to the Free
21** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22*/
23
24/* CVS $Id: rpoly.hpp,v 1.15 2014/01/30 17:49:27 cxsc Exp $ */
25
26//============================================================================
27//
28// Program/Module
29// from
30// C++ TOOLBOX FOR VERIFIED COMPUTING I
31// Basic Numerical Problems
32//
33// Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz
34//
35// For details on theory, algorithms, and programs, see the book
36//
37// R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for
38// Verified Computing I - Basic Numerical Problems. Springer-Verlag,
39// Heidelberg, New York, 1995.
40//
41//============================================================================
42//----------------------------------------------------------------------------
43// File: rpoly (header)
44// Purpose: Definition of the class for the representation of a real
45// polynomial by its coefficients.
46// Class RPolynomial:
47// Deg() : to get the degree of the polynomial
48// RPolynomial(): constructors
49// operator [] : component access
50// operator >> : input operator for data of type RPolynomial
51// operator << : output operator for data of type RPolynomial
52//----------------------------------------------------------------------------
53#ifndef __RPOLY_HPP
54#define __RPOLY_HPP
55
56#include <rvector.hpp> // Real vector arithmetic
57
58using namespace cxsc;
59using namespace std;
60
61class RPolynomial {
62 private:
63 rvector coeff;
64 public:
65 RPolynomial ( int );
66 RPolynomial ( const RPolynomial& );
67 real& operator[] ( int i ) { return coeff[i]; }
68
69 friend int Deg ( const RPolynomial& );
70 friend istream& operator>> ( istream&, RPolynomial& );
71 friend ostream& operator<< ( ostream&, RPolynomial );
72};
73
74int Deg ( const RPolynomial& );
75istream& operator>> ( istream&, RPolynomial& );
76ostream& operator<< ( ostream&, RPolynomial );
77#endif
78
79
80
81
The Scalar Type real.
Definition real.hpp:114
The Data Type rvector.
Definition rvector.hpp:58
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29

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.

1/*
2** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3**
4** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5** Universitaet Karlsruhe, Germany
6** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7** Universitaet Wuppertal, Germany
8**
9** This library is free software; you can redistribute it and/or
10** modify it under the terms of the GNU Library General Public
11** License as published by the Free Software Foundation; either
12** version 2 of the License, or (at your option) any later version.
13**
14** This library is distributed in the hope that it will be useful,
15** but WITHOUT ANY WARRANTY; without even the implied warranty of
16** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17** Library General Public License for more details.
18**
19** You should have received a copy of the GNU Library General Public
20** License along with this library; if not, write to the Free
21** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22*/
23
24/* CVS $Id: rpoly.cpp,v 1.15 2014/01/30 17:49:27 cxsc Exp $ */
25
26//============================================================================
27//
28// Program/Module
29// from
30// C++ TOOLBOX FOR VERIFIED COMPUTING I
31// Basic Numerical Problems
32//
33// Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz
34//
35// For details on theory, algorithms, and programs, see the book
36//
37// R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for
38// Verified Computing I - Basic Numerical Problems. Springer-Verlag,
39// Heidelberg, New York, 1995.
40//
41//============================================================================
42//----------------------------------------------------------------------------
43// File: rpoly (implementation)
44// Purpose: Definition of the class for the representation of a real
45// polynomial by its coefficients.
46// Class RPolynomial:
47// Deg() : to get the degree of the polynomial
48// RPolynomial(): constructors
49// operator [] : component access
50// operator >> : input operator for data of type RPolynomial
51// operator << : output operator for data of type RPolynomial
52//----------------------------------------------------------------------------
53#include <rpoly.hpp>
54
55using namespace cxsc;
56using namespace std;
57
58int Deg ( const RPolynomial& p ) { return Ub(p.coeff); }
59
60RPolynomial::RPolynomial( int n )
61{
62 Resize(coeff,0,n);
63 coeff = 0.0;
64}
65
66RPolynomial::RPolynomial( const RPolynomial& p )
67{
68 Resize(coeff,0,Deg(p));
69 coeff = p.coeff;
70}
71
72istream& operator>> ( istream& in, RPolynomial& p )
73{
74 cout << " x^0 * "; in >> p[0];
75 for (int i = 1; i <= Deg(p); i++)
76 { cout << "+ x^" << i << " * "; in >> p[i]; }
77 return in;
78}
79
80ostream& operator<< ( ostream& out, RPolynomial p )
81{
82 int PolyIsZero, n = Deg(p);
83
84 PolyIsZero = 1;
85 for (int i = 0; i <= n; i++) {
86 if (p[i] == 0.0) continue;
87 if (PolyIsZero)
88 out << " ";
89 else
90 out << "+ ";
91 out << p[i] << " * x^" << i << endl;
92 PolyIsZero = 0;
93 }
94 if (PolyIsZero) out << " 0 (= zero polynomial)" << endl;
95 return out;
96}
97
98
99
100
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
void Resize(cimatrix &A) noexcept
Resizes the matrix.

Module rpeval

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)
3**
4** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5** Universitaet Karlsruhe, Germany
6** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7** Universitaet Wuppertal, Germany
8**
9** This library is free software; you can redistribute it and/or
10** modify it under the terms of the GNU Library General Public
11** License as published by the Free Software Foundation; either
12** version 2 of the License, or (at your option) any later version.
13**
14** This library is distributed in the hope that it will be useful,
15** but WITHOUT ANY WARRANTY; without even the implied warranty of
16** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17** Library General Public License for more details.
18**
19** You should have received a copy of the GNU Library General Public
20** License along with this library; if not, write to the Free
21** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22*/
23
24/* CVS $Id: rpeval.hpp,v 1.14 2014/01/30 17:49:27 cxsc Exp $ */
25
26//============================================================================
27//
28// Program/Module
29// from
30// C++ TOOLBOX FOR VERIFIED COMPUTING I
31// Basic Numerical Problems
32//
33// Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz
34//
35// For details on theory, algorithms, and programs, see the book
36//
37// R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for
38// Verified Computing I - Basic Numerical Problems. Springer-Verlag,
39// Heidelberg, New York, 1995.
40//
41//============================================================================
42//----------------------------------------------------------------------------
43// File: rpeval (header)
44// Purpose: Evaluation of a real polynomial p with maximum accuracy.
45// Global functions:
46// RPolyEval() : computes an enclosure for the exact value p(t) with
47// maximum accuracy
48// RPolyEvalErrMsg(): delivers an error message text
49//----------------------------------------------------------------------------
50#ifndef __RPEVAL_HPP
51#define __RPEVAL_HPP
52
53#include <interval.hpp> // Interval arithmetic
54#include <rpoly.hpp> // Real polynomial type
55
56using namespace cxsc;
57using namespace std;
58
59extern char* RPolyEvalErrMsg ( int );
60extern void RPolyEval ( RPolynomial, real, real&, interval&, int&, int& );
61#endif
62
63
64
65
The Scalar Type interval.
Definition interval.hpp:55

The function RPolyEval() is an implementation of the algorithm. It computes an approximation and an enclosure of the value of the real polynomial $ p(t) = \sum \limits_{i=0}^n p_i t^i $ with $ p_i \in R $ at a point $ t \in R $.

1/*
2** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3**
4** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5** Universitaet Karlsruhe, Germany
6** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7** Universitaet Wuppertal, Germany
8**
9** This library is free software; you can redistribute it and/or
10** modify it under the terms of the GNU Library General Public
11** License as published by the Free Software Foundation; either
12** version 2 of the License, or (at your option) any later version.
13**
14** This library is distributed in the hope that it will be useful,
15** but WITHOUT ANY WARRANTY; without even the implied warranty of
16** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17** Library General Public License for more details.
18**
19** You should have received a copy of the GNU Library General Public
20** License along with this library; if not, write to the Free
21** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22*/
23
24/* CVS $Id: rpeval.cpp,v 1.14 2014/01/30 17:49:27 cxsc Exp $ */
25
26//============================================================================
27//
28// Program/Module
29// from
30// C++ TOOLBOX FOR VERIFIED COMPUTING I
31// Basic Numerical Problems
32//
33// Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz
34//
35// For details on theory, algorithms, and programs, see the book
36//
37// R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for
38// Verified Computing I - Basic Numerical Problems. Springer-Verlag,
39// Heidelberg, New York, 1995.
40//
41//============================================================================
42//----------------------------------------------------------------------------
43// File: rpeval (implementation)
44// Purpose: Evaluation of a real polynomial p with maximum accuracy.
45// Method: Transformation of Horner's scheme for evaluating p(t) to a linear
46// system of equations A(t)*x = p, where A(t) is a bidiagonal, Toeplitz
47// matrix and x = (x_n,...,x_0). By iterative refinement the floating-
48// point approximation is improved and the exact solution x is enclosed in
49// an interval vector [x].
50// Global functions:
51// RPolyEval() : computes an enclosure for the exact value p(t) with
52// maximum accuracy
53// RPolyEvalErrMsg(): delivers an error message text
54//----------------------------------------------------------------------------
55#include <string.h> // String handling
56#include <imatrix.hpp> // Interval matrix/vector arithmetic
57#include <mvi_util.hpp> // Interval matrix/vector utility functions
58#include <rpeval.hpp>
59
60using namespace cxsc;
61using namespace std;
62
63const int kmax = 10; // Maximum number of iteration steps
64
65//---------------------------------------------------------------------------
66// Error codes used in this module. In the comments below p[0],..., p[n] are
67// the coefficients of a polynomial p.
68//---------------------------------------------------------------------------
69const int
70 NoError = 0, // No error occurred
71 ItFailed = 1; // Maximum number of iterations exceeded
72
73//----------------------------------------------------------------------------
74// Error messages depending on the error code.
75//----------------------------------------------------------------------------
76char* RPolyEvalErrMsg ( int Err )
77{
78 static char Msg[80] = "";
79
80 if (Err != NoError) {
81 char Hlp[60];
82
83 if (Err == ItFailed)
84 sprintf(Hlp,"Maximum number of iterations (=%d) exceeded",kmax);
85 else
86 strcpy(Hlp,"Code not defined");
87 sprintf(Msg,"Error: %s!",Hlp);
88 }
89 return(Msg);
90}
91
92//----------------------------------------------------------------------------
93// Purpose: Determination of p(t) (a polynomial p with argument t) with
94// maximum accuracy.
95// Parameters:
96// In : 'p' : represents a real polynomial by its coefficients
97// 't' : specifies the point of evaluation
98// Out: 'z' : floating-point approximation computed by Horner's scheme
99// 'zz' : enclosure of p(t)
100// 'k' : number of iterations needed
101// 'Err': error flag
102// Description: The polynomial 'p' is defined by its coefficients and, 't'
103// denotes the evaluation point. Horner's scheme for evaluating p is equi-
104// valent to computing the solution of a linear system of equations
105// A(t)*x = p, where A(t) is a bidiagonal, Toeplitz matrix and
106// x =(x_n,...,x_0). The component x_0 of x is equal to p(t). The solution
107// x is enclosed in an interval vector [x] by iterative refinement. The
108// first element of [x] is an enclosure of p(t). It is returned in the
109// variable 'zz'. A floating-point approximation computed by Horner's
110// scheme is returned in 'z'. The number of iterations is returned in 'k'
111// and the state of success is returned in 'Err'.
112// Remark: The polynomial's coefficients and the argument are assumed to be
113// exact floating-point numbers!
114//----------------------------------------------------------------------------
115void RPolyEval ( RPolynomial p, real t, real& z,
116 interval& zz, int& k, int& Err )
117{
118 int i, j, n = Deg(p); // The j-th column of the matrix
119 ivector rr(0,n), yy(0,n); // 'y' is used to store the j-th
120 rvector x(0,n); // correction for the solution of
121 rmatrix y(0,n,0,kmax); // A(t)*y = p.
122 dotprecision Accu; //-------------------------------
123 idotprecision IAccu;
124
125 k = 0; // Initialization
126 Err = NoError; //---------------
127
128 if (n == 0) // Constant polynomial
129 { z = p[0]; zz = p[0]; } //--------------------
130
131 else if (n == 1) { // Linear polynomial
132 Accu = p[0]; //------------------
133 accumulate(Accu,t,p[1]);
134 rnd(Accu,z);
135 rnd(Accu,zz);
136 }
137 else { // n > 1
138 x = 0.0; y = 0.0; // Initialization x = 0 and y = 0
139 //--------------------------------
140
141 x[n] = p[n]; // Computation of a first approxi-
142 for (i = n-1; i >= 0; i--) // mation using Horner's scheme
143 x[i] = x[i+1]*t + p[i]; //--------------------------------
144 z = x[0];
145
146 y[Col(0)] = x; // Iterative refinement for the
147 do { // solution of A*x = p
148 //--------------------------------
149 if (k > 0) // If a residual was computed, its middle
150 y[Col(k)] = mid(yy); // is stored as the next correction of 'y'
151
152 yy[n] = 0.0; // Computation of the residual [r]
153 for (i = n-1; i >= 0; i--) { // and evaluation of the interval
154 Accu = p[i]; // system A*[y] = [r]
155 for(j = 0; j <= k; j++) Accu -= y[i][j];
156 for(j = 0; j <= k; j++) accumulate(Accu,t,y[i+1][j]);
157 rnd(Accu,rr[i]);
158 yy[i] = yy[i+1]*t + rr[i];
159 }
160
161 IAccu = yy[0]; // Determination of a new
162 for (j = 0; j <= k; j++) IAccu += y[0][j]; // enclosure [z] of p(t)
163 rnd(IAccu,zz);
164
165 k++;
166 } while ( !(UlpAcc(zz,1) || (k > kmax)) );
167
168 if ( !UlpAcc(zz,1) ) Err = ItFailed;
169 }
170}
171
172
173
174
The Data Type dotprecision.
Definition dot.hpp:112
The Data Type idotprecision.
Definition idot.hpp:48
The Data Type ivector.
Definition ivector.hpp:55
The Data Type rmatrix.
Definition rmatrix.hpp:471
cimatrix_subv Col(cimatrix &m, const int &i) noexcept
Returns one column of the matrix as a vector.
Definition cimatrix.inl:242
int UlpAcc(const interval &x, int n)
Checks if the diameter of the interval is ulps.
Definition interval.cpp:335
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739

Examples

We consider the polynomial $ p(t) = t^4 - 8t^3 + 24t^2 - 32t + 16 $ to be evaluated in the neighborhood of the real value $ t=2.0001 $ and the polynomial $ q(t) = -t^3 + 3t^2 - 3t + 1 $ to be evaluated in the neighborhood of the real value $ t = 1.000005 $. To make sure that the arguments are representable on the computer, we use the machine numbers $ t = (2.0001) $ and $ t = (1.000005) $, 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().

1//============================================================================
2//
3// Program/Module
4// from
5// C++ TOOLBOX FOR VERIFIED COMPUTING I
6// Basic Numerical Problems
7//
8// Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz
9//
10// For details on theory, algorithms, and programs, see the book
11//
12// R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for
13// Verified Computing I - Basic Numerical Problems. Springer-Verlag,
14// Heidelberg, New York, 1995.
15//
16//============================================================================
17//----------------------------------------------------------------------------
18// Example: Evaluation of a real polynomial with maximum accuracy
19//----------------------------------------------------------------------------
20#include <rpeval.hpp> // Polynomial evaluation
21
22
23using namespace cxsc;
24using namespace std;
25
26
27int main ( )
28{
29 int Err, No, n;
30 real t, y;
31 interval yy;
32
33 do {
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
49 if (!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 }
57 else
58 cout << RPolyEvalErrMsg(Err) << endl;
59
60 return 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 $ p(t) $ and $ q(t)$. 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 $ p( 2.0001) $ and $ q( 1.000005) $. 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 $ t = (2.0001) $ and $ t = (1.000005) $.