defined as
We present an algorithm to evaluate a real polynomial 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 \]](form_86.png)
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.
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
![$ rr[i] $](form_92.png)
denotes the interval ![$ [r]_i $](form_93.png)
.
. Since arithmetic operations for real polynomials are not required by the algorithm, they are not provided by module rpoly.
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2011 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: rpoly.hpp,v 1.10 2011/06/07 15:17:33 cxsc Exp $ */ 00025 00026 //============================================================================ 00027 // 00028 // Program/Module 00029 // from 00030 // C++ TOOLBOX FOR VERIFIED COMPUTING I 00031 // Basic Numerical Problems 00032 // 00033 // Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz 00034 // 00035 // This program/module is free software for non-commercial use. For details 00036 // on theory, algorithms, and programs, see the book 00037 // 00038 // R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for 00039 // Verified Computing I - Basic Numerical Problems. Springer-Verlag, 00040 // Heidelberg, New York, 1995. 00041 // 00042 // This program/module is distributed WITHOUT ANY WARRANTY. For details, 00043 // see the "Disclaimer / Legal Matters" of the book (page iv). 00044 // 00045 //============================================================================ 00046 //---------------------------------------------------------------------------- 00047 // File: rpoly (header) 00048 // Purpose: Definition of the class for the representation of a real 00049 // polynomial by its coefficients. 00050 // Class RPolynomial: 00051 // Deg() : to get the degree of the polynomial 00052 // RPolynomial(): constructors 00053 // operator [] : component access 00054 // operator >> : input operator for data of type RPolynomial 00055 // operator << : output operator for data of type RPolynomial 00056 //---------------------------------------------------------------------------- 00057 #ifndef __RPOLY_HPP 00058 #define __RPOLY_HPP 00059 00060 #include <rvector.hpp> // Real vector arithmetic 00061 00062 using namespace cxsc; 00063 using namespace std; 00064 00065 class RPolynomial { 00066 private: 00067 rvector coeff; 00068 public: 00069 RPolynomial ( int ); 00070 RPolynomial ( const RPolynomial& ); 00071 real& operator[] ( int i ) { return coeff[i]; } 00072 00073 friend int Deg ( const RPolynomial& ); 00074 friend istream& operator>> ( istream&, RPolynomial& ); 00075 friend ostream& operator<< ( ostream&, RPolynomial ); 00076 }; 00077 00078 int Deg ( const RPolynomial& ); 00079 istream& operator>> ( istream&, RPolynomial& ); 00080 ostream& operator<< ( ostream&, RPolynomial ); 00081 #endif 00082 00083 00084 00085
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.
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2011 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: rpoly.cpp,v 1.10 2011/06/07 15:17:33 cxsc Exp $ */ 00025 00026 //============================================================================ 00027 // 00028 // Program/Module 00029 // from 00030 // C++ TOOLBOX FOR VERIFIED COMPUTING I 00031 // Basic Numerical Problems 00032 // 00033 // Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz 00034 // 00035 // This program/module is free software for non-commercial use. For details 00036 // on theory, algorithms, and programs, see the book 00037 // 00038 // R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for 00039 // Verified Computing I - Basic Numerical Problems. Springer-Verlag, 00040 // Heidelberg, New York, 1995. 00041 // 00042 // This program/module is distributed WITHOUT ANY WARRANTY. For details, 00043 // see the "Disclaimer / Legal Matters" of the book (page iv). 00044 // 00045 //============================================================================ 00046 //---------------------------------------------------------------------------- 00047 // File: rpoly (implementation) 00048 // Purpose: Definition of the class for the representation of a real 00049 // polynomial by its coefficients. 00050 // Class RPolynomial: 00051 // Deg() : to get the degree of the polynomial 00052 // RPolynomial(): constructors 00053 // operator [] : component access 00054 // operator >> : input operator for data of type RPolynomial 00055 // operator << : output operator for data of type RPolynomial 00056 //---------------------------------------------------------------------------- 00057 #include <rpoly.hpp> 00058 00059 using namespace cxsc; 00060 using namespace std; 00061 00062 int Deg ( const RPolynomial& p ) { return Ub(p.coeff); } 00063 00064 RPolynomial::RPolynomial( int n ) 00065 { 00066 Resize(coeff,0,n); 00067 coeff = 0.0; 00068 } 00069 00070 RPolynomial::RPolynomial( const RPolynomial& p ) 00071 { 00072 Resize(coeff,0,Deg(p)); 00073 coeff = p.coeff; 00074 } 00075 00076 istream& operator>> ( istream& in, RPolynomial& p ) 00077 { 00078 cout << " x^0 * "; in >> p[0]; 00079 for (int i = 1; i <= Deg(p); i++) 00080 { cout << "+ x^" << i << " * "; in >> p[i]; } 00081 return in; 00082 } 00083 00084 ostream& operator<< ( ostream& out, RPolynomial p ) 00085 { 00086 int PolyIsZero, n = Deg(p); 00087 00088 PolyIsZero = 1; 00089 for (int i = 0; i <= n; i++) { 00090 if (p[i] == 0.0) continue; 00091 if (PolyIsZero) 00092 out << " "; 00093 else 00094 out << "+ "; 00095 out << p[i] << " * x^" << i << endl; 00096 PolyIsZero = 0; 00097 } 00098 if (PolyIsZero) out << " 0 (= zero polynomial)" << endl; 00099 return out; 00100 } 00101 00102 00103 00104
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2011 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: rpeval.hpp,v 1.9 2011/06/07 15:17:33 cxsc Exp $ */ 00025 00026 //============================================================================ 00027 // 00028 // Program/Module 00029 // from 00030 // C++ TOOLBOX FOR VERIFIED COMPUTING I 00031 // Basic Numerical Problems 00032 // 00033 // Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz 00034 // 00035 // This program/module is free software for non-commercial use. For details 00036 // on theory, algorithms, and programs, see the book 00037 // 00038 // R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for 00039 // Verified Computing I - Basic Numerical Problems. Springer-Verlag, 00040 // Heidelberg, New York, 1995. 00041 // 00042 // This program/module is distributed WITHOUT ANY WARRANTY. For details, 00043 // see the "Disclaimer / Legal Matters" of the book (page iv). 00044 // 00045 //============================================================================ 00046 //---------------------------------------------------------------------------- 00047 // File: rpeval (header) 00048 // Purpose: Evaluation of a real polynomial p with maximum accuracy. 00049 // Global functions: 00050 // RPolyEval() : computes an enclosure for the exact value p(t) with 00051 // maximum accuracy 00052 // RPolyEvalErrMsg(): delivers an error message text 00053 //---------------------------------------------------------------------------- 00054 #ifndef __RPEVAL_HPP 00055 #define __RPEVAL_HPP 00056 00057 #include <interval.hpp> // Interval arithmetic 00058 #include <rpoly.hpp> // Real polynomial type 00059 00060 using namespace cxsc; 00061 using namespace std; 00062 00063 extern char* RPolyEvalErrMsg ( int ); 00064 extern void RPolyEval ( RPolynomial, real, real&, interval&, int&, int& ); 00065 #endif 00066 00067 00068 00069
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 
.
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2011 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: rpeval.cpp,v 1.9 2011/06/07 15:17:33 cxsc Exp $ */ 00025 00026 //============================================================================ 00027 // 00028 // Program/Module 00029 // from 00030 // C++ TOOLBOX FOR VERIFIED COMPUTING I 00031 // Basic Numerical Problems 00032 // 00033 // Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz 00034 // 00035 // This program/module is free software for non-commercial use. For details 00036 // on theory, algorithms, and programs, see the book 00037 // 00038 // R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for 00039 // Verified Computing I - Basic Numerical Problems. Springer-Verlag, 00040 // Heidelberg, New York, 1995. 00041 // 00042 // This program/module is distributed WITHOUT ANY WARRANTY. For details, 00043 // see the "Disclaimer / Legal Matters" of the book (page iv). 00044 // 00045 //============================================================================ 00046 //---------------------------------------------------------------------------- 00047 // File: rpeval (implementation) 00048 // Purpose: Evaluation of a real polynomial p with maximum accuracy. 00049 // Method: Transformation of Horner's scheme for evaluating p(t) to a linear 00050 // system of equations A(t)*x = p, where A(t) is a bidiagonal, Toeplitz 00051 // matrix and x = (x_n,...,x_0). By iterative refinement the floating- 00052 // point approximation is improved and the exact solution x is enclosed in 00053 // an interval vector [x]. 00054 // Global functions: 00055 // RPolyEval() : computes an enclosure for the exact value p(t) with 00056 // maximum accuracy 00057 // RPolyEvalErrMsg(): delivers an error message text 00058 //---------------------------------------------------------------------------- 00059 #include <string.h> // String handling 00060 #include <imatrix.hpp> // Interval matrix/vector arithmetic 00061 #include <mvi_util.hpp> // Interval matrix/vector utility functions 00062 #include <rpeval.hpp> 00063 00064 using namespace cxsc; 00065 using namespace std; 00066 00067 const int kmax = 10; // Maximum number of iteration steps 00068 00069 //--------------------------------------------------------------------------- 00070 // Error codes used in this module. In the comments below p[0],..., p[n] are 00071 // the coefficients of a polynomial p. 00072 //--------------------------------------------------------------------------- 00073 const int 00074 NoError = 0, // No error occurred 00075 ItFailed = 1; // Maximum number of iterations exceeded 00076 00077 //---------------------------------------------------------------------------- 00078 // Error messages depending on the error code. 00079 //---------------------------------------------------------------------------- 00080 char* RPolyEvalErrMsg ( int Err ) 00081 { 00082 static char Msg[80] = ""; 00083 00084 if (Err != NoError) { 00085 char Hlp[60]; 00086 00087 if (Err == ItFailed) 00088 sprintf(Hlp,"Maximum number of iterations (=%d) exceeded",kmax); 00089 else 00090 strcpy(Hlp,"Code not defined"); 00091 sprintf(Msg,"Error: %s!",Hlp); 00092 } 00093 return(Msg); 00094 } 00095 00096 //---------------------------------------------------------------------------- 00097 // Purpose: Determination of p(t) (a polynomial p with argument t) with 00098 // maximum accuracy. 00099 // Parameters: 00100 // In : 'p' : represents a real polynomial by its coefficients 00101 // 't' : specifies the point of evaluation 00102 // Out: 'z' : floating-point approximation computed by Horner's scheme 00103 // 'zz' : enclosure of p(t) 00104 // 'k' : number of iterations needed 00105 // 'Err': error flag 00106 // Description: The polynomial 'p' is defined by its coefficients and, 't' 00107 // denotes the evaluation point. Horner's scheme for evaluating p is equi- 00108 // valent to computing the solution of a linear system of equations 00109 // A(t)*x = p, where A(t) is a bidiagonal, Toeplitz matrix and 00110 // x =(x_n,...,x_0). The component x_0 of x is equal to p(t). The solution 00111 // x is enclosed in an interval vector [x] by iterative refinement. The 00112 // first element of [x] is an enclosure of p(t). It is returned in the 00113 // variable 'zz'. A floating-point approximation computed by Horner's 00114 // scheme is returned in 'z'. The number of iterations is returned in 'k' 00115 // and the state of success is returned in 'Err'. 00116 // Remark: The polynomial's coefficients and the argument are assumed to be 00117 // exact floating-point numbers! 00118 //---------------------------------------------------------------------------- 00119 void RPolyEval ( RPolynomial p, real t, real& z, 00120 interval& zz, int& k, int& Err ) 00121 { 00122 int i, j, n = Deg(p); // The j-th column of the matrix 00123 ivector rr(0,n), yy(0,n); // 'y' is used to store the j-th 00124 rvector x(0,n); // correction for the solution of 00125 rmatrix y(0,n,0,kmax); // A(t)*y = p. 00126 dotprecision Accu; //------------------------------- 00127 idotprecision IAccu; 00128 00129 k = 0; // Initialization 00130 Err = NoError; //--------------- 00131 00132 if (n == 0) // Constant polynomial 00133 { z = p[0]; zz = p[0]; } //-------------------- 00134 00135 else if (n == 1) { // Linear polynomial 00136 Accu = p[0]; //------------------ 00137 accumulate(Accu,t,p[1]); 00138 rnd(Accu,z); 00139 rnd(Accu,zz); 00140 } 00141 else { // n > 1 00142 x = 0.0; y = 0.0; // Initialization x = 0 and y = 0 00143 //-------------------------------- 00144 00145 x[n] = p[n]; // Computation of a first approxi- 00146 for (i = n-1; i >= 0; i--) // mation using Horner's scheme 00147 x[i] = x[i+1]*t + p[i]; //-------------------------------- 00148 z = x[0]; 00149 00150 y[Col(0)] = x; // Iterative refinement for the 00151 do { // solution of A*x = p 00152 //-------------------------------- 00153 if (k > 0) // If a residual was computed, its middle 00154 y[Col(k)] = mid(yy); // is stored as the next correction of 'y' 00155 00156 yy[n] = 0.0; // Computation of the residual [r] 00157 for (i = n-1; i >= 0; i--) { // and evaluation of the interval 00158 Accu = p[i]; // system A*[y] = [r] 00159 for(j = 0; j <= k; j++) Accu -= y[i][j]; 00160 for(j = 0; j <= k; j++) accumulate(Accu,t,y[i+1][j]); 00161 rnd(Accu,rr[i]); 00162 yy[i] = yy[i+1]*t + rr[i]; 00163 } 00164 00165 IAccu = yy[0]; // Determination of a new 00166 for (j = 0; j <= k; j++) IAccu += y[0][j]; // enclosure [z] of p(t) 00167 rnd(IAccu,zz); 00168 00169 k++; 00170 } while ( !(UlpAcc(zz,1) || (k > kmax)) ); 00171 00172 if ( !UlpAcc(zz,1) ) Err = ItFailed; 00173 } 00174 } 00175 00176 00177 00178
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]
00001 //============================================================================ 00002 // 00003 // Program/Module 00004 // from 00005 // C++ TOOLBOX FOR VERIFIED COMPUTING I 00006 // Basic Numerical Problems 00007 // 00008 // Copyright (c) 1995 Rolf Hammer, Matthias Hocks, Dietmar Ratz 00009 // 00010 // This program/module is free software for non-commercial use. For details 00011 // on theory, algorithms, and programs, see the book 00012 // 00013 // R. Hammer, M. Hocks, U. Kulisch, D. Ratz: C++ Toolbox for 00014 // Verified Computing I - Basic Numerical Problems. Springer-Verlag, 00015 // Heidelberg, New York, 1995. 00016 // 00017 // This program/module is distributed WITHOUT ANY WARRANTY. For details, 00018 // see the "Disclaimer / Legal Matters" of the book (page iv). 00019 // 00020 //============================================================================ 00021 //---------------------------------------------------------------------------- 00022 // Example: Evaluation of a real polynomial with maximum accuracy 00023 //---------------------------------------------------------------------------- 00024 #include <rpeval.hpp> // Polynomial evaluation 00025 00026 00027 using namespace cxsc; 00028 using namespace std; 00029 00030 00031 int main ( ) 00032 { 00033 int Err, No, n; 00034 real t, y; 00035 interval yy; 00036 00037 do { 00038 cout << "Enter the degree of the polynomial (>=0): "; 00039 cin >> n; cout << endl; 00040 } while (n < 0); 00041 00042 RPolynomial p(n); 00043 00044 cout << SetPrecision(23,15) << Scientific; // Output format 00045 00046 cout << "Enter the coefficients of p in increasing order:" << endl; 00047 cin >> p; cout << endl; 00048 cout << "Enter the argument t for evaluation:" << endl; 00049 cin >> t; cout << endl; 00050 00051 RPolyEval(p,t,y,yy,No,Err); 00052 00053 if (!Err) { 00054 cout << "Polynomial:" << endl << p <<endl; 00055 cout << "Floating-point evaluation of p(t) using Horner's scheme:" 00056 << endl << " " << y << endl << endl; 00057 cout << "Verified inclusion of p(t):" 00058 << endl << yy << endl << endl; 00059 cout << "Number of iterations needed: " << No << endl; 00060 } 00061 else 00062 cout << RPolyEvalErrMsg(Err) << endl; 00063 00064 return 0; 00065 }
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 .
1.4.6