This example shows how to design a function which takes a polynomial and a value and returns the sign of this polynomial at this point. This function is a filter: if the answer is not guaranteed, the functions says so. The reason of using a filter rather than a simple evaluation function is: computations with floating-point numbers will incur approximations and it can be enough to change the sign of the polynomial. So, in order to validate the result, the function will use interval arithmetic.
The first step is the inclusion of the appropriate headers. Because the function will handle floating-point bounds, the easiest solution is:
#include <boost/numeric/interval.hpp>
Now, let's begin the function. The polynomial is given by the array of its coefficients and its size (strictly greater to its degree). In order to simplify the code, two namespaces of the library are included.
int sign_polynomial(double x, double P[], int sz) {
using namespace boost::numeric;
using namespace interval_lib;
Then we can define the interval type. Since no special behavior is required, the default policies are enough:
typedef interval<double> I;
For the evaluation, let's just use the Horner scheme with interval
arithmetic. The library overloads all the arithmetic operators and provides
mixed operations, so the only difference between the code with and without
interval arithmetic lies in the type of the iterated value y:
I y = P[sz - 1];
for(int i = sz - 2; i >= 0; i--)
y = y * x + P[i];
The last step is the computation of the sign of y. It is done
by choosing an appropriate comparison scheme and then doing the comparison
with the usual operators:
using namespace compare::certain; if (y > 0.) return 1; if (y < 0.) return -1; return 0; }
The answer 0 does not mean the polynomial is zero at this
point. It only means the answer is not known since y contains
zero and thus does not have a precise sign.
Now we have the expected function. However, due to the poor implementations of floating-point rounding in most of the processors, it can be useful to say to optimize the code; or rather, to let the library optimize it. The main condition for this optimization is that the interval code should not be mixed with floating-point code. In this example, it is the case, since all the operations done in the functions involve the library. So the code can be rewritten:
int sign_polynomial(double x, double P[], int sz) {
using namespace boost::numeric;
using namespace interval_lib;
typedef interval<double> I_aux;
I_aux::traits_type::rounding rnd;
typedef unprotect<I_aux>::type I;
I y = P[sz - 1];
for(int i = sz - 2; i >= 0; i--)
y = y * x + P[i];
using namespace compare::certain;
if (y > 0.) return 1;
if (y < 0.) return -1;
return 0;
}
The difference between this code and the previous is the use of another
interval type. This new type I indicates to the library that all
the computations can be done without caring for the rounding mode. And
because of that, it is up to the function to care about it: a rounding object
need to be alive whenever the optimized type is used.
In libs/numeric/interval/test/ and
libs/numeric/interval/examples/ are some test and example
programs.. The examples illustrate a few uses of intervals. For a general
description and considerations on using this library, and some potential
domains of application, please read this mini-guide.
The test programs are as follows. Please note that they require the use of
the Boost.test library and can be automatically tested by using
bjam (except for interval_test.cpp).
add.cpp tests if the additive and subtractive operators and the respective _std and _opp rounding functions are correctly implemented. It is done by using symbolic expressions as a base type.
cmp.cpp, cmp_lex.cpp, cmp_set.cpp, and
cmp_tribool.cpp test if the operators <
> <= >= ==
!= behave correctly for the default, lexicographic, set, and
tristate comparisons. cmp_exp.cpp tests the explicit comparison
functions cer.. and pos.. behave correctly.
cmp_exn.cpp tests if the various policies correctly detect exceptional
cases. All these tests use some simple intervals ([1,2] and [3,4], [1,3] and
[2,4], [1,2] and [2,3], etc).
det.cpp tests if the _std and _opp
versions in protected and unprotected mode produce the same result when Gauss
scheme is used on an unstable matrix (in order to exercise rounding). The
tests are done for interval<float> and
interval<double>.
fmod.cpp defines a minimalistic version of
interval<int> and uses it in order to test
fmod on some specific interval values.
mul.cpp exercises the multiplication, the finite division, the square and the square root with some integer intervals leading to exact results.
pi.cpp tests if the interval value of π (for
int, float and double base types)
contains the number π (defined with 21 decimal digits) and if it is a
subset of [π±1ulp] (in order to ensure some precision).
pow.cpp tests if the pow function behaves correctly on
some simple test cases.
test_float.cpp exercises the arithmetic operations of the library for floating point base types.
interval_test.cpp tests if the interval library respects the
inclusion property of interval arithmetic by computing some functions and
operations for both double and
interval<double>.
filter.cpp contains filters for computational geometry able to find the sign of a determinant. This example is inspired by the article Interval arithmetic yields efficient dynamic filters for computational geometry by Brönnimann, Burnikel and Pion, 2001.
findroot_demo.cpp finds zeros of some functions by using dichotomy and even produces gnuplot data for one of them. The processor has to correctly handle elementary functions for this example to properly work.
horner.cpp is a really basic example of unprotecting the interval operations for a whole function (which computes the value of a polynomial by using Horner scheme).
io.cpp shows some stream input and output operators for intervals .The wide variety of possibilities explains why the library do not implement i/o operators and they are left to the user.
newton-raphson.cpp is an implementation of a specialized version of Newton-Raphson algorithm for finding the zeros of a function knowing its derivative. It exercises unprotecting, full division, some set operations and empty intervals.
transc.cpp implements the transcendental part of the rounding
policy for double by using an external library (the MPFR subset
of GMP in this case).
Revised: 2003-08-16
Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
Polytechnic University.
Copyright (c) Guillaume Melquiond, 2003.