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 
58 using namespace cxsc;
59 using namespace std;
60 
61 class 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 
74 int Deg ( const RPolynomial& );
75 istream& operator>> ( istream&, RPolynomial& );
76 ostream& operator<< ( ostream&, RPolynomial );
77 #endif
78 
79 
80 
81 

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 
55 using namespace cxsc;
56 using namespace std;
57 
58 int Deg ( const RPolynomial& p ) { return Ub(p.coeff); }
59 
60 RPolynomial::RPolynomial( int n )
61 {
62  Resize(coeff,0,n);
63  coeff = 0.0;
64 }
65 
66 RPolynomial::RPolynomial( const RPolynomial& p )
67 {
68  Resize(coeff,0,Deg(p));
69  coeff = p.coeff;
70 }
71 
72 istream& 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 
80 ostream& 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 

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 
56 using namespace cxsc;
57 using namespace std;
58 
59 extern char* RPolyEvalErrMsg ( int );
60 extern void RPolyEval ( RPolynomial, real, real&, interval&, int&, int& );
61 #endif
62 
63 
64 
65 

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 
60 using namespace cxsc;
61 using namespace std;
62 
63 const 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 //---------------------------------------------------------------------------
69 const 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 //----------------------------------------------------------------------------
76 char* 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 //----------------------------------------------------------------------------
115 void 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 

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 
23 using namespace cxsc;
24 using namespace std;
25 
26 
27 int 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) $.

cxsc::mid
cvector mid(const cimatrix_subv &mv)
Returns the middle of the matrix.
Definition: cimatrix.inl:739
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:54
cxsc::rmatrix
The Data Type rmatrix.
Definition: rmatrix.hpp:470
cxsc::rvector
The Data Type rvector.
Definition: rvector.hpp:57
cxsc::idotprecision
The Data Type idotprecision.
Definition: idot.hpp:47
cxsc::in
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
Definition: cinterval.cpp:654
cxsc::Col
cimatrix_subv Col(cimatrix &m, const int &i)
Returns one column of the matrix as a vector.
Definition: cimatrix.inl:242
cxsc::dotprecision
The Data Type dotprecision.
Definition: dot.hpp:111
cxsc::ivector
The Data Type ivector.
Definition: ivector.hpp:54
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::UlpAcc
int UlpAcc(const interval &x, int n)
Checks if the diameter of the interval is ulps.
Definition: interval.cpp:335
cxsc::Ub
int Ub(const cimatrix &rm, const int &i)
Returns the upper bound index.
Definition: cimatrix.inl:1163
cxsc::Resize
void Resize(cimatrix &A)
Resizes the matrix.
Definition: cimatrix.inl:1211
cxsc::real
The Scalar Type real.
Definition: real.hpp:113