commit 21e378f40b250205a086ef202d6d602e8e4f6bcf Author: Guillaume Melquiond Date: Tue Jan 21 17:20:32 2003 +0000 Import of Interval library [SVN r16980] diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..3e84d7c --- /dev/null +++ b/.gitattributes @@ -0,0 +1,96 @@ +* text=auto !eol svneol=native#text/plain +*.gitattributes text svneol=native#text/plain + +# Scriptish formats +*.bat text svneol=native#text/plain +*.bsh text svneol=native#text/x-beanshell +*.cgi text svneol=native#text/plain +*.cmd text svneol=native#text/plain +*.js text svneol=native#text/javascript +*.php text svneol=native#text/x-php +*.pl text svneol=native#text/x-perl +*.pm text svneol=native#text/x-perl +*.py text svneol=native#text/x-python +*.sh eol=lf svneol=LF#text/x-sh +configure eol=lf svneol=LF#text/x-sh + +# Image formats +*.bmp binary svneol=unset#image/bmp +*.gif binary svneol=unset#image/gif +*.ico binary svneol=unset#image/ico +*.jpeg binary svneol=unset#image/jpeg +*.jpg binary svneol=unset#image/jpeg +*.png binary svneol=unset#image/png +*.tif binary svneol=unset#image/tiff +*.tiff binary svneol=unset#image/tiff +*.svg text svneol=native#image/svg%2Bxml + +# Data formats +*.pdf binary svneol=unset#application/pdf +*.avi binary svneol=unset#video/avi +*.doc binary svneol=unset#application/msword +*.dsp text svneol=crlf#text/plain +*.dsw text svneol=crlf#text/plain +*.eps binary svneol=unset#application/postscript +*.gz binary svneol=unset#application/gzip +*.mov binary svneol=unset#video/quicktime +*.mp3 binary svneol=unset#audio/mpeg +*.ppt binary svneol=unset#application/vnd.ms-powerpoint +*.ps binary svneol=unset#application/postscript +*.psd binary svneol=unset#application/photoshop +*.rdf binary svneol=unset#text/rdf +*.rss text svneol=unset#text/xml +*.rtf binary svneol=unset#text/rtf +*.sln text svneol=native#text/plain +*.swf binary svneol=unset#application/x-shockwave-flash +*.tgz binary svneol=unset#application/gzip +*.vcproj text svneol=native#text/xml +*.vcxproj text svneol=native#text/xml +*.vsprops text svneol=native#text/xml +*.wav binary svneol=unset#audio/wav +*.xls binary svneol=unset#application/vnd.ms-excel +*.zip binary svneol=unset#application/zip + +# Text formats +.htaccess text svneol=native#text/plain +*.bbk text svneol=native#text/xml +*.cmake text svneol=native#text/plain +*.css text svneol=native#text/css +*.dtd text svneol=native#text/xml +*.htm text svneol=native#text/html +*.html text svneol=native#text/html +*.ini text svneol=native#text/plain +*.log text svneol=native#text/plain +*.mak text svneol=native#text/plain +*.qbk text svneol=native#text/plain +*.rst text svneol=native#text/plain +*.sql text svneol=native#text/x-sql +*.txt text svneol=native#text/plain +*.xhtml text svneol=native#text/xhtml%2Bxml +*.xml text svneol=native#text/xml +*.xsd text svneol=native#text/xml +*.xsl text svneol=native#text/xml +*.xslt text svneol=native#text/xml +*.xul text svneol=native#text/xul +*.yml text svneol=native#text/plain +boost-no-inspect text svneol=native#text/plain +CHANGES text svneol=native#text/plain +COPYING text svneol=native#text/plain +INSTALL text svneol=native#text/plain +Jamfile text svneol=native#text/plain +Jamroot text svneol=native#text/plain +Jamfile.v2 text svneol=native#text/plain +Jamrules text svneol=native#text/plain +Makefile* text svneol=native#text/plain +README text svneol=native#text/plain +TODO text svneol=native#text/plain + +# Code formats +*.c text svneol=native#text/plain +*.cpp text svneol=native#text/plain +*.h text svneol=native#text/plain +*.hpp text svneol=native#text/plain +*.ipp text svneol=native#text/plain +*.tpp text svneol=native#text/plain +*.jam text svneol=native#text/plain +*.java text svneol=native#text/plain diff --git a/doc/checking.htm b/doc/checking.htm new file mode 100644 index 0000000..6e82af5 --- /dev/null +++ b/doc/checking.htm @@ -0,0 +1,218 @@ + + + + + + Checking policies + + + +

Checking policies

+ +

A checking policy controls how the interval class will deal +with special cases like: empty intervals, infinite numbers, invalid +values.

+ +

For example, let's consider operator+(interval, T). The +second argument could be an invalid value (for a floating-point number, it is +a NaN). What to do in such a case? First, we could say that the second +argument can never be an invalid number. Second, we could also say such a +situation can arise but is forbidden. Third, we could allow such values and +generate an empty interval when encountered. And there is many other +possibilities.

+ +

It is the reason why such a policy is used: there is a lot of interesting +behaviors and it would be sad to arbitrarily select one of these.

+ +

Requirements

+ +

The checking class should satisfy the following requirement (in the form +of an interface):

+
/* requirements for checking policy */
+struct checking
+{
+  static T inf();
+  static T nan();
+  static bool is_nan(const T&);
+  static T empty_lower();
+  static T empty_upper();
+  static bool is_empty(const T&, const T&);
+};
+ +

The first function, inf is invoked each time the library need +to create the infinite bound of an interval. For example, +interval::whole computes interval(-checking::inf(), +checking::inf()). If infinite values are allowed and +std::numeric_limits<T>::infinity() returns a correct +value, such a value can be used.

+ +

Next comes nan. This function is used each time a function +need to return a value of type T but is unable to compute it. It +only happens when one of the arguments of the function is invalid. For +example, if you ask what the median value of an empty interval is, +nan will be used. But please remember: lower and +upper directly return the value stocked in the interval; so, if +the interval is empty, lower will not answer by a +call to checking::nan (but will return the same value than +checking::empty_lower could return).

+ +

empty_lower and empty_upper respectively return +the lower and upper bound of the empty interval. There is no requirements for +empty_lower and empty_upper to return the same +value than checking::nan. For example, if the type +T does not have any invalid value, the empty_ +functions can return the [1;0] interval.

+ +

is_nan is used to test if a value of type T is +invalid or not. is_empty tests if the interval formed by the two +arguments is empty or not. Such tests will generally be at the beginning of +each function which involves an argument of type T. If one of +the inputs is declared invalid, the the function will try to produce an +invalid value or an input interval.

+ +

Synopsis

+
namespace boost {
+  namespace interval_lib {
+
+template<class T>
+struct checking_base;
+template<class T, class Checking = checking_base<T>, class Exception = exception_create_empty<T> >
+struct checking_no_empty;
+template<class T, class Checking = checking_base<T> >
+struct checking_no_nan;
+template<class T, class Checking = checking_base<T>, class Exception = exception_invalid_number<T> >
+struct checking_catch_nan;
+
+template<class T> struct exception_create_empty { T operator()(); };
+template<class T> struct exception_invalid_number { void operator()(); };
+
+  } // namespace interval_lib
+} // namespace boost
+ +

Predefined classes

+ +

In order to simplify the customization of the policy, some templates are +already defined in the library.

+ +

First of all, there is checking_base. Thanks to the +information provided by std::numeric_limits<T>, this class +is able to generate a base for the policy. If T has quiet NaNs +(as said by numeric_limits::has_quiet_NaN), then the value is +used for nan, empty_lower, +empty_upper; and a basic test is used for is_nan +(it is x!=x). If T does not have quiet NaNs, then +nan is an assert(false), the empty interval is +[1,0], and is_nan always return false. As for +nan, inf returns +numeric_limits::infinity() if possible, or is an +assert(false) otherwise. Finally, is_empty(T l,T u) +is always defined by !(l<=u).

+ +

Next comes checking_no_empty. Using it means that each time +an empty interval should be produced (by empty_lower and +empty_upper), the function object given by the +Exception argument of the template is invoked and the value it +returns is propagated. So, if Exception is appropriately defined +(for example it could throw an exception, hence the name of the argument), +you can be sure no empty interval will ever be created. So +is_empty will always return false (since there is +no need to test for an empty interval). And as explained before, in that case +we can also replace nan by an assert(false); you +will be sure no invalid number will ever be produced. If this template is not +used, it implicitly means that all the functions can produce empty intervals +and they correctly deal with empty interval arguments.

+ +

Finally there are checking_no_nan and +checking_catch_nan. The first one expresses the functions of the +library will never get an invalid number as argument. So is_nan +will only return false. The other one means the arguments can be +an invalid number but in that case, is_nan will call the +function object Exception and return false. Indeed, +this template means invalid numbers should never make their way through to +the body of the function. If none of this two templates is used, it +implicitly means that all the functions can get invalid number arguments and +they will correctly deal with them.

+ +

exception_create_empty throws std::runtime_error +with the message "boost::interval: empty interval created" and +exception_invalid_number throws +std::invalid_argument with the message "boost::interval: +invalid number".

+ +

Customizing your own checking policy

+ +

In order to define a suitable policy, you need to correctly say what you +expect from your interval class. First of all, are you interested in getting +empty intervals at the end of a calculus? If you do not want to obtain empty +intervals, empty_lower and empty_upper have to fail +when invoked (they can throw an exception, set a flag, etc). However, if no +function is able to produce an empty interval, it is no more necessary to do +the test, so is_empty may always return false. In +this case, a good compiler will do a lot of optimizations.

+ +

You could also be interested in getting empty intervals at the end of the +calculus. For example, if you need to transform an array of unsure values (or +intervals) in a new array of intervals, you may not want to stop the +conversion at the first encountered problem. So empty_lower and +empty_upper need to return suitable values in order to define an +empty interval (you can use an upper bound which is not greater or equal than +the lower bound for example); and is_empty must be able to +distinguish empty intervals from the valid intervals.

+ +

Another important question is: is it possible that some base numbers +(objects of type T) are invalid? And if it is possible, are they +allowed or not ? If it is not possible, no test is necessary; +is_nan may always return false. In this case too, a +good compiler will do a lot of optimizations. If function arguments can hold +invalid numbers, two cases must be considered according to whether they are +allowed or not. If they are allowed, is_nan just has to test if +they are invalid or not. If they are forbidden, is_nan should +fail (exception, assert, etc.) when invoked on an invalid argument and return +false otherwise. The value returned by nan does not +have any interest since the interval functions are guaranteed not to produce +invalid interval bounds unless the user passes invalid numbers to the +constructors. So you can put an assert inside if you do not trust the +library. :-)

+ +

And finally, you need to decide what to do with nan if it has +not already been decided at the beginning, and with inf. These +functions should return a value or start an exceptional behavior (especially +if the base type does not have corresponding values).

+ +

Some examples

+ +
+ +

Revised: 2003-01-19
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/comparisons.htm b/doc/comparisons.htm new file mode 100644 index 0000000..a02dba7 --- /dev/null +++ b/doc/comparisons.htm @@ -0,0 +1,162 @@ + + + + + + Comparisons + + + +

Comparisons

+ +

As was said before, the definition of the comparison operators induces a +slight problem. There are many ways to define them, depending of the return +type or the expected order. It is the reason why the meaning of the operators +is not fixed once and for all.

+ +

The way the operators are defined could have been influenced by a policy, +as it is already the case for the rounding and the checking. However, +comparisons are more an external property of the the class rather than an +internal one. They are meant to be locally modified, independantly of the +type of the intervals.

+ +

The operators <, <=, >, +>=, ==, != are defined each time; +and like the arithmetic operators they can take an argument of the base type. +However, due to technical limitations, this base type can only be the second +argument; so the operators are unfortunately not fully symmetric. The return +type is not always bool, since some interesting results can be +achieved by using a tri-state return type. So here is the common signatures +of the operators:

+
template<class T, class Policies1, class Policies2>
+return_type operator== (const interval<T, Policies1>&, const interval<T, Policies2>&);
+
+template<class T, class Policies>
+return_type operator== (const interval<T, Policies>&, const T&);
+ +

Provided comparisons

+ +

Default comparison

+ +

If nothing is specified, the meaning of the comparison operators are an +extension of the operator on the base type. More precisely, if one of the +argument is invalid or empty, an exception is thrown. If the arguments are +valid, the following rules are applied to determine the result of +[a,b] op [c,d] (just consider +c == d if the second argument is of type +T):

+ + +

This comparison allows to replace base types by interval types without +changing the meaning of a program. Indeed, if no exception is thrown, the +result is the same as before; and if an exception is thrown, the previous +comparison was unsure and should have been rewritten.

+ +

Other comparisons

+ +

The other comparisons are selected by using a namespace. These namespaces +are located under boost::interval_lib::compare and are invoked +by:

+
using namespace boost::interval_lib::compare::the_comparison_to_select;
+ +

After this line, the default meaning of the operators will have been +replaced by the meaning located in the namespace. Please note that because of +C++ lookup rules, it is not possible to use two namespaces one after another +and they must be used in different block hierarchies. Otherwise the compiler +will complain about ambiguous operators. To summarize:

+
// example 1: BAD
+using namespace compare1;
+...
+using namespace compare2;
+...
+
+// example 2: GOOD
+{ using namespace compare1;
+  ...                       }
+{ using namespace compare2;
+  ...                       }
+
+// example 3: BAD
+using namespace compare1;
+...
+{ using namespace compare2;
+  ...                       }
+ +

Now comes the list of the provided comparisons. They all are located in +their respective header files under +<boost/numeric/interval/compare/>. And as for the default +comparison, the operators will generally complain by throwing an exception if +feed by invalid values.

+ + +

Exception

+ +

+
namespace boost {
+  namespace interval_lib {
+
+class comparison_error: std::runtime_error; // "boost::interval: uncertain comparison"
+
+  }
+}
+ +

Explicit comparison functions

+ +

In some situation, you may want to perform direct comparisons on the +bounds and avoid the indeterminate case that appears with default operators. +Some functions are provided for this purpose. They expect their arguments to +be valid and return a result after only one comparison. Their names are +composed by cer (for "certain", if the default comparison is +true, the result is true) or pos (for "possible", if the default +comparison is false, the result is false) followed by lt, +le, gt, ge, eq or +ne. They are located in +<boost/numeric/interval/compare/explicit.hpp>. Each of these +functions takes two parameters and returns a boolean; the parameters are +expected to be valid, undefined behavior may result otherwise. For example, +the definition of the "certainly less than" comparison is:

+
namespace boost {
+  namespace interval_lib {
+
+template<class T, class Policies1, class Policies2>
+bool cerlt(const interval<T, Policies1>& x, const interval<T, Policies2>& y);
+
+template<class T, class Policies>
+bool cerlt(const interval<T, Policies>& x, const T& y);
+
+template<class T, class Policies>
+bool cerlt(const T& x, const interval<T, Policies>& y);
+
+  } // namespace interval_lib
+} // namespace boost
+
+ +

Revised: 2002-10-13
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/examples.htm b/doc/examples.htm new file mode 100644 index 0000000..ad6c563 --- /dev/null +++ b/doc/examples.htm @@ -0,0 +1,93 @@ + + + + + + Tests and Examples + + + +

Tests and Examples

+ +

In libs/interval/test/ and +libs/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.

+ +

Tests

+ +

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 substractive 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 and cmp_lex.cpp test if the operators +< > <= >= +== != behave correctly for the default comparison +and the lexicographic comparison. cmp_exp.cpp tests the explicit +comparison functions cer.. and pos.. behave +correctly. 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 exercices 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>.

+ +

Examples

+ +

findroot_demo.cpp finds zeros of some functions by using dichotomy +and even produces gnuplot datas 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).

+ +

newton-raphson.cpp is a (bad) 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.

+ +

interval_speed.cpp is not really an example but more like a +benchmark for testing performance of the library. It needs to be compiled +with -DUSE_BOOST. The compilation option +-DUSE_BOOST_UNPROTECTED allows the program to use unprotected +local rounding, leading to better performance.

+
+ +

Revised: 2002-12-15
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/guide.htm b/doc/guide.htm new file mode 100644 index 0000000..cc71b53 --- /dev/null +++ b/doc/guide.htm @@ -0,0 +1,108 @@ + + + + + + Choosing Your Own Interval Type + + + +

Choosing Your Own Interval Type

+ +

First of all, you need to select your base type. In order to obtain an +useful interval type, the numbers should respect some requirements. Please +refer to this page in order to see them. When your +base type is robust enough, you can go to the next step: the choice of the +policies.

+ +

As you should already know if you did not come to this page by accident, +the interval class expect a policies argument describing the rounding and checking +policies. The first thing to do is to verify if the default policies are or +are not adapted to your case. If your base type is not float, +double, or long double, the default rounding policy +is probably not adapted. However, by specializing +interval_lib::rounded_math to your base type, the default +rounding policy will be suitable.

+ +

The default policies define an interval type that performs precise +computations (for float, double, long +double), detects invalid numbers and throws exception each times an +empty interval is created. This is a brief description and you should refer +to the corresponding sections for a more precise description of the default +policies. Unless you need some special behavior, this default type is usable +in a lot of situations.

+ +

After having completely defined the interval type (and its policies), the +only thing left to do is to verify that the constants are defined and +std::numeric_limits is correct (if needed). Now you can use your +brand new interval type.

+ +

Some Examples

+ +

Solving systems

+ +

If you use the interval library in order to solve equation and inequation +systems by bisection, something like +boost::interval<double> is probably what you need. The +computations are precise, and they may be fast if enclosed in a protected +rounding mode block (see the performance +section). The comparison are "certain"; it is probably the most used type of +comparison, and the other comparisons are still accessible by the explicit +comparison functions. The checking forbid empty interval; they are not needed +since there would be an empty interval at end of the computation if an empty +interval is created during the computation, and no root would be inside. The +checking also forbid invalid numbers (NaN for floating-point numbers). It can +be a minor performance hit if you only use exact floating-point constants +(which are clearly not NaNs); however, if performance really does matter, you +will probably use a good compiler which knows how to inline functions and all +these annoying little tests will magically disappear (if not, it is time to +upgrade your compiler).

+ +

Manipulating wide intervals

+ +

You may want to use the library on intervals with imprecise bounds or on +inexact numbers. In particular, it may be an existing algorithm that you want +to rewrite and simplify by using the library. In that case, you are not +really interested by the inclusion property; you are only interested by the +computation algorithms the library provides. So you do not need to use any +rounding; the checking also may not be useful. Use an "exact computation" +rounding (you are allowed to think the name stangely applies to the +situation) and a checking that never tests for any invalid numbers or empty +intervals. By doing that, you will obtain library functions reduced to their +minimum (an addition of two intervals will only be two additions of +numbers).

+ +

Computing ranges

+ +

The inputs of your program may be empty intervals or invalid values (for +example, a database can allow undefined values in some field) and the core of +your program could also do some non-arithmetic computations that do not +always propagate empty intervals. For example, in the library, the +hull function can happily receive an empty interval but not +generate an empty interval if the other input is valid. The +intersect function is also able to produce empty intervals if +the intervals do not overlap. In that case, it is not really interesting if +an exception is thrown each time an empty interval is produced or an invalid +value is used; it would be better to generate and propagate empty intervals. +So you need to change the checking policy to something like +interval_lib::checking_base<T>.

+ +

Switching interval types

+ +

This example does not deal with a full case, but with a situation that can +occur often. Sometimes, it can be useful to change the policies of an +interval by converting it to another type. For example, this happens when you +use an unprotected version of the interval type in order to speed up the +computations; it is a change of the rounding policy. It also happens when you +want to temporarily allow empty intervals to be created; it is a change of +the checking policy. These changes should not be prohibited: they can greatly +enhance a program (lisibility, interest, performance).

+
+ +

Revised: 2003-01-15
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/includes.htm b/doc/includes.htm new file mode 100644 index 0000000..376c502 --- /dev/null +++ b/doc/includes.htm @@ -0,0 +1,169 @@ + + + + + + Headers Inclusion + + + +

Headers Inclusion

+ +

The easiest way to access to the library is by including the main +header:

+
#include <boost/numeric/interval.hpp>
+ +

This header will include almost all the other headers (except the ones +listed as extensions). However, you may not want to access all the +functionalities of the library. So this page stands as a reminder for the +whole structure of the library. <boost/numeric/interval.hpp> +is the only header to be located directly under boost/numeric/; +all the other headers are located in the subdirectory +boost/numeric/interval/. And each time this documentation will +refer to interval/something.hpp, it is +<boost/numeric/interval/something.hpp>.

+ +

Please also note that all the following headers are independent and can +easily be pre-compiled if necessary (for compilers which support pre-compiled +headers of course).

+ +

Function definitions

+ +

The following headers contain the definition of the interval +class and all the friendly functions and operators.

+ +

interval/interval.hpp

+ +

This header contains the definition and the declaration of the +interval class. However, this class is templated and the default +template parameters are not available by this header. In particular, this +header does not provide the default specialization of the +interval class for the floating-point types +(interval<float>, interval<double> and +interval<long double>). So, unless you use your own +policies, this header is not really useful on its own.

+ +

interval/utility.hpp

+ +

In this header are all the functions that do not expect any arithmetic +property from the base number type. It only expects the bounds to be ordered; +but it should not surprise you since it is a requirement of the whole +library. You will find in this header the definitions of access and related +functions: lower, upper, +checked_lower, checked_upper, median, +width, widen. There are also the set-like +functions: in, in_zero, empty, +subset, proper_subset, overlap, +singleton, equal, intersect, +hull, bisect. Finally, abs, +min and max are defined.

+ +

interval/arith.hpp

+ +

Here are the binary operators +, -, +*, / and the unary operator -.

+ +

interval/arith2.hpp

+ +

This header defines fmod, square, +sqrt and pow.

+ +

interval/transc.hpp

+ +

It is the last of the three headers with mathematical functions; it +provides the following functions: cos, sin, +tan, acos, asin, atan, +cosh, sinh, tanh, acosh, +asinh, atanh, exp and +log.

+ +

Policies

+ +

The following headers define some policies. They may be needed if you use +the default policies.

+ +

interval/rounded_arith.hpp

+ +

This header defines the three provided rounding policies for the +arithmetic functions: rounded_arith_std, +rounded_arith_opp, rounded_arith_exact.

+ +

interval/rounded_transc.hpp

+ +

This header defines the three provided rounding policies for the +transcendental functions: rounded_transc_std, +rounded_transc_opp, rounded_transc_exact. It is +separated from rounded_arith.hpp since the transcendental part +of the rounding policy is probably less useful than the arithmetic part.

+ +

interval/hw_rounding.hpp

+ +

Here are full rounding policies for the basic floating-point types. The +policies are processor-dependent; and to allow the user code to be portable, +they only define the common subset of the hardware available functions, which +are the arithmetic functions of the rounding policy.

+ +

interval/checking.hpp

+ +

This header provides the predefined checking policies: +checking_base, checking_no_empty, +checking_no_nan, checking_catch_nan, +checking_strict.

+ +

interval/policies.hpp

+ +

Here are defined the helpers for manipulating policies. It contains +policies (and so is needed for using default policies), +change_rounding, change_checking, +unprotect, etc.

+ +

Comparisons

+ +

interval/compare.hpp

+ +

This header includes all the following headers. They provide some +predefined comparison namespaces.

+ +

interval/compare/explicit.hpp

+ +

The explicit comparison functions cerlt, posge, +etc are defined in this header.

+ +

interval/compare/lexicographic.hpp

+ +

This header provides compare::lexicographic.

+ +

interval/compare/set.hpp

+ +

This header provides compare::set.

+ +

Extensions

+ +

The following headers are not included by interval.hpp and +will usually provide not always desirable capabilities.

+ +

interval/io.hpp

+ +

Here are defined basic stream operators << and +>>. They should only be used as a first approach and later +be replaced by a customized version.

+ +

interval/compare/tribool.hpp

+ +

This header provides a comparison namespace compare::tribool +especially adapted to a tristate boolean.

+ +

interval/ext/x86_fast_rounding_control.hpp

+ +

This header defines a new rounding policy allowing to workaround the +precision problem of the x86 processors (and so speeding up the +computations). However, it only is a partial solution and it shouldn't be +used when there is a possibility of underflow.

+
+ +

Revised: 2003-01-21
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/index.html b/doc/index.html new file mode 100644 index 0000000..305e74e --- /dev/null +++ b/doc/index.html @@ -0,0 +1,9 @@ + + + + + +Automatic redirection failed, please go to +interval.htm. + + diff --git a/doc/interval.htm b/doc/interval.htm new file mode 100644 index 0000000..c692386 --- /dev/null +++ b/doc/interval.htm @@ -0,0 +1,864 @@ + + + + + + Boost Interval Arithmetic Library + + + +

c++boost.gif (8819 bytes) Interval Arithmetic Library

+ +
+ + + + + + + +
Contents of this page:
+ Introduction
+ Synopsis
+ Template class interval
+ Operations and functions
+ Interval support library
+ Compilation notes
+ Common pitfalls and dangers
+ Rationale
+ History and Acknowledgments		Other pages associated with this page:
+ Rounding policies
+ Checking policies
+ Policies manipulation
+ Comparisons
+ Base number type requirements
+ Choosing your own interval type
+ Test and example programs
+ Headers inclusion
+
+ +

Introduction and Overview

+ +

As implied by its name, this library is intended to help manipulating +mathematical intervals. It consists of a single header <boost/numeric/interval.hpp> +and principally +a type which can be used as interval<T>. In fact, this +interval template is declared as interval<T,Policies> where +Policies is a policy class that controls the various behaviours of +the interval class; interval<T> just happens to pick the +default policies for the type T.

+ +

Interval Arithmetic

+ +

An interval is a pair of numbers which represents all the numbers between +these two. (Intervals are considered close so the bounds are included.) The +purpose of this library is to extend the usual arithmetic functions to +intervals. These intervals will be written [a,b] to represent +all the numbers between a and b (included). a and +b can be infinite (but they can not be the same infinite) and a +≤ b.

+ +

The fundamental property of interval arithmetic is the +inclusion property:

+
+
``if f is a function on a set of numbers, f can be + extended to a new function defined on intervals. This new function + f takes one interval argument and returns an interval result + such as: ∀ x ∈ [a,b], + f(x) ∈ f([a,b]).''
+
+ +

Such a property is not limited to functions with only one argument. +Whenever possible, the interval result should be the smallest one able to +satisfy the property (it is not really useful if the new functions always +answer [-∞,+∞]).

+ +

There are at least two reasons a user would like to use this library. The +obvious one is when the user has to compute with intervals. One example is +when input data have some builtin imprecision: instead of a number, an input +variable can be passed as an interval. Another example application is to +solve equations, by bisecting an interval until the interval width is small +enough. A third example application is in computer graphics, where +computations with boxes, segments or rays can be reduced to computations with +points via intervals.

+ +

Another common reason to use interval arithmetic is when the computer +doesn't produce exact results: by using intervals, it is possible to quantify +the propagation of rounding errors. This approach is used often in numerical +computation. For example, let's assume the computer stores numbers with ten +decimal significant digits. To the question 1 + 1E-100 - 1, the computer will +answer 0 although the correct answer would be 1E-100. With the help of +interval arithmetic, the computer will answer [0,1E-9]. This is quite a huge +interval for such a little result, but the precision is now known, without +having to compute error propagation.

+ +

Numbers, rounding, and exceptional behavior

+ +

The base number type is the type that holds the +bounds of the interval. In order to successfully use interval arithmetic, the +base number type must present some characteristics. Firstly, due to the definition of an +interval, the base numbers have to be totally ordered so, for instance, +complex<T> is not usable as base number type for +intervals. The mathematical functions for the base number type should also +be compatible with the total order (for instance if x>y and z>t, then +it should also hold that x+z > y+t), so modulo types are not usable +either.

+ +

Secondly, the computations must be exact or provide some rounding methods +(for instance, toward minus or plus infinity) if we want to guarantee the +inclusion property. Note that we also may explicitely specify no rounding, +for instance if the base number type is exact, i.e. the result of a +mathematic operation is always computed and represented without loss of +precision. If the number type is not exact, we may still explicitely specify +no rounding, with the obvious consequence that the inclusion property is no +longuer guaranteed.

+ +

Finally, because heavy loss of precision is always possible, some numbers +have to represent infinities or an exceptional behavior must be defined. The +same situation also occurs for NaN (Not a Number).

+ +

Given all this, one may want to limit the template argument T of the class +template interval to the floating point types +float, double, and long double, as +defined by the IEEE-754 Standard. Indeed, if the interval arithmetic is +intended to replace the arithmetic provided by the floating point unit of a +processor, these types are the best choice. Unlike std::complex, +however, we don't want to limit T to these types. This is why we +allow the rounding and exceptional behaviors to be given by the two policies +(rounding and checking). We do nevertheless provide highly optimized rounding +and checking class specializations for the above-mentioned floating point +types.

+ +

Operations and functions

+ +

It is straightforward to define the elementary arithmetic operations on +intervals, being guided by the inclusion property. For instance, if [a,b] and +[c,d] are intervals, [a,b]+[c,d] can be computed by taking the smallest +interval that contains all the numbers x+y for x in [a,b] and y in [c,d]; in +this case, rounding a+b down and c+d up will suffice. Other operators and +functions are similarly defined (see their definitions below).

+ +

Comparisons

+ +

It is also possible to define some comparison operators. Given two +intervals, the result is a tri-state boolean type +{false,true,indeterminate}. The answers false and +true are easy to manipulate since they can directly be mapped on the +boolean true and false. But it is not the case for the answer +indeterminate since comparison operators are supposed to be boolean +functions. So, what to do in order to obtain boolean answers?

+ +

One solution consists of deciding to adopt an exceptional behavior, such +as a failed assertion or raising an exception. In this case, the exceptional +behavior will be triggered when the result is indeterminate.

+ +

Another solution is to map indeterminate always to false, +or always to true. If false is chosen, the comparison will be +called "certain;" indeed, the result of [a,b] < +[c,d] will be true if and only if: ∀ x +∈ [a,b] ∀ y ∈ [c,d], +x < y. If true is chosen, the comparison will be +called "possible;" indeed, the result of [a,b] < +[c,d] will be true if and only if: ∃ x +∈ [a,b] ∃ y ∈ [c,d], +x < y.

+ +

Since any of these solution has a clearly defined semantics, it is not +clear that we should enforce either of them. For this reason, the default +behavior consists to mimic the real comparisons by throwing an exception in +the indeterminate case. Other behaviors can be selected bu using specific +comparison namespace. There is also a bunch of explicitely named comparison +functions. See comparisons pages for further +details.

+ +

Overview of the library, and usage

+ +

This library provides two quite distinct levels of usage. One is to use +the basic class template interval<T> without specifying +the policy. This only requires to know and understand the concepts developed +above and the content of the namespace boost. In addition to the class +interval<T>, this level of usage provides arithmetic +operators (+, -, *, /), +algebraic and piecewise-algebraic functions (abs, +square, sqrt, pow), transcendental and +trigonometric functions (exp, log, +sin, cos, tan, asin, +acos, atan, sinh, cosh, +tanh, asinh, acosh, +atanh), and the standard comparison operators +(<, <=, >, +>=, ==, !=), as well as several +interval-specific functions (min, max, which have a +different meaning than std::min and std::max; +lower, upper, width, +median, empty, singleton, equal, +in, in_zero, subset, +proper_subset, overlap, intersection, +hull, bisect).

+ +

For some functions which take several parameters of type +interval<T>, all combinations of argument types +T and interval<T> which contain at least one +interval<T>, are considered in order to avoid a conversion +from the arguments of type T to a singleton of type +interval<T>. This is done for efficiency reasons (the fact +that an argument is a singleton sometimes renders some tests unnecessary).

+ +

A somewhat more advanced usage of this library is to hand-pick the +policies Rounding and Checking and pass them to +interval<T, Policies> through the use of Policies := +boost::interval_lib::policies<Rounding,Checking>. Appropriate +policies can be fabricated by using the various classes provided in the +namespace boost::interval_lib as detailed in section Interval Support Library. It is also possible to +choose the comparison scheme by overloading operators through namespaces.

+ +

Synopsis

+
namespace boost {
+namespace numeric {
+
+namespace interval_lib {
+
+/* this declaration is necessary for the declaration of interval */
+template <class T> struct default_policies;
+
+/* ... ; the full synopsis of namespace interval_lib can be found here */
+
+} // namespace interval_lib
+
+
+/* template interval_policies; class definition can be found here */
+template<class Rounding, class Checking>
+struct interval_policies;
+
+/* template class interval; class definition can be found here */
+template<class T, class Policies = typename interval_lib::default_policies<T> >
+class interval;
+
+/* arithmetic operators involving intervals */
+template <class T, class Traits>  interval<T, Traits> operator+(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> operator-(const interval<T, Traits>& x);
+
+template <class T, class Traits>  interval<T, Traits> operator+(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> operator+(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> operator+(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  interval<T, Traits> operator-(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> operator-(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> operator-(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  interval<T, Traits> operator*(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> operator*(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> operator*(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  interval<T, Traits> operator/(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> operator/(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> operator/(const T& r, const interval<T, Traits>& x);
+
+/* algebraic functions: sqrt, abs, square, pow */
+template <class T, class Traits>  interval<T, Traits> abs(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> sqrt(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> square(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> pow(const interval<T, Traits>& x, int y);
+
+/* transcendental functions: exp, log */
+template <class T, class Traits>  interval<T, Traits> exp(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> log(const interval<T, Traits>& x);
+
+/* fmod, for trigonometric function argument reduction (see below) */
+template <class T, class Traits>  interval<T, Traits> fmod(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> fmod(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> fmod(const T& x, const interval<T, Traits>& y);
+
+/* trigonometric functions */
+template <class T, class Traits>  interval<T, Traits> sin(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> cos(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> tan(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> asin(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> acos(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> atan(const interval<T, Traits>& x);
+
+/* hyperbolic trigonometric functions */
+template <class T, class Traits>  interval<T, Traits> sinh(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> cosh(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> tanh(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> asinh(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> acosh(const interval<T, Traits>& x);
+template <class T, class Traits>  interval<T, Traits> atanh(const interval<T, Traits>& x);
+
+/* min, max external functions (NOT std::min/max, see below) */
+template <class T, class Traits>  interval<T, Traits> max(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> max(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> max(const T& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> min(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> min(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> min(const T& x, const interval<T, Traits>& y);
+
+/* bounds-related interval functions */
+template <class T, class Traits>  T lower(const interval<T, Traits>& x);
+template <class T, class Traits>  T upper(const interval<T, Traits>& x);
+template <class T, class Traits>  T width(const interval<T, Traits>& x);
+template <class T, class Traits>  T median(const interval<T, Traits>& x);
+
+/* bounds-related interval functions */
+template <class T, class Traits>  bool empty(const interval<T, Traits>& b);
+template <class T, class Traits>  bool singleton(const interval<T, Traits>& x);
+template <class T, class Traits>  bool equal(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool in(const T& r, const interval<T, Traits>& b);
+template <class T, class Traits>  bool in_zero(const interval<T, Traits>& b);
+template <class T, class Traits>  bool subset(const interval<T, Traits>& a, const interval<T, Traits>& b);
+template <class T, class Traits>  bool proper_subset(const interval<T, Traits>& a, const interval<T, Traits>& b);
+template <class T, class Traits>  bool overlap(const interval<T, Traits>& x, const interval<T, Traits>& y);
+
+/* set manipulation interval functions */
+template <class T, class Traits>  interval<T, Traits> intersection(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> hull(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> hull(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  interval<T, Traits> hull(const T& x, const interval<T, Traits>& y);
+template <class T, class Traits>  interval<T, Traits> hull(const T& x, const T& y);
+template <class T, class Traits>  std::pair<interval<T, Traits>, interval<T, Traits> > bisect(const interval<T, Traits>& x);
+
+/* interval comparison operators */
+template<class T, class Traits>  bool operator<(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template<class T, class Traits>  bool operator<(const interval<T, Traits>& x, const T& y);
+template<class T, class Traits>  bool operator<(const T& x, const interval<T, Traits>& y);
+
+template<class T, class Traits>  bool operator<=(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template<class T, class Traits>  bool operator<=(const interval<T, Traits>& x, const T& y);
+template<class T, class Traits>  bool operator<=(const T& x, const interval<T, Traits>& y);
+
+template<class T, class Traits>  bool operator>(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template<class T, class Traits>  bool operator>(const interval<T, Traits>& x, const T& y);
+template<class T, class Traits>  bool operator>(const T& x, const interval<T, Traits>& y);
+
+template<class T, class Traits>  bool operator>=(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template<class T, class Traits>  bool operator>=(const interval<T, Traits>& x, const T& y);
+template<class T, class Traits>  bool operator>=(const T& x, const interval<T, Traits>& y);
+
template<class T, class Traits>  bool operator==(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template<class T, class Traits>  bool operator==(const interval<T, Traits>& x, const T& y);
+template<class T, class Traits>  bool operator==(const T& x, const interval<T, Traits>& y);
+
+template<class T, class Traits>  bool operator!=(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template<class T, class Traits>  bool operator!=(const interval<T, Traits>& x, const T& y);
+template<class T, class Traits>  bool operator!=(const T& x, const interval<T, Traits>& y);
+
+} // namespace numeric
+} // namespace boost
+ +

Template class interval

+The public interface of the template class interval itself is kept at a +simplest minimum: +
template <class T, class Policies = typename interval_lib::default_policies<T>::type>
+class interval
+{
+  public:
+  typedef T      base_type;
+  typedef Traits traits_type;
+
+  interval(const T& v = T());
+  interval(const T& l, const T& u);
+
+  template<class Traits2>
+  interval(const interval<T,Traits2>& r);
+
+  // compiler-generated copy constructor and assignment operator are fine
+
+  interval& operator=(const T& x);
+  void assign(const T& l, const T& u);
+
+  const T& lower() const;
+  const T& upper() const;
+
+  static interval empty();
+  static interval whole();
+  static interval hull(const T& x, const T& y);
+
+  interval& operator+= (const T& r);
+  interval& operator+= (const interval& r);
+  interval& operator-= (const T& r);
+  interval& operator-= (const interval& r);
+  interval& operator*= (const T& r);
+  interval& operator*= (const interval& r);
+  interval& operator/= (const T& r);
+  interval& operator/= (const interval& r);
+};
+ +

The two main constructors accept one and two arguments of type T. The +first one produces a singleton interval [v,v]; and the second one an interval +[l,u]. The second one will use the checking policy if u<l. If you are not +sure whether the bounds are inverted or not, it is better to use the function +call hull(l,u). This will not create a problem if u<l.

+ +

There is a template constructor in order to change the traits parameter of +an interval. But there is no constructor that changes the base type of an +interval (it could lead to violating the inclusion property if the old base +type is not contained in the new one).

+ +

There is an assign function in order to directly change the bounds of an +interval. It behaves like the second constructor if the bounds are inverted. +There is no assign function that directly takes an interval or only one +number as a parameter; just use the assignment operator in that case.

+ +

Assignment operator, copy constructor and destructor are the ones the +compiler automatically generates. There is also an assignment operator for +the type T.

+ +

The static functions empty and whole produces +the corresponding intervals. They are static member functions rather than +global functions because they cannot guess their return types. Likewise for +hull.

+ +

Operations and Functions

+ +

Some of the following functions expect min and +max to be defined for the base type. Those are the only +requirements for the interval class (but the policies can have +other requirements).

+ +

Operators + - * / ++= -= *= /=

+ +

The basic operations are the unary minus and the binary + +- * /. The unary minus takes an +interval and returns an interval. The binary operations take two intervals, +or one interval and a number, and return an interval. If an argument is a +number instead of an interval, you can expect the result to be the same as if +the number was first converted to an interval. This property will be true for +all the following functions and operators.

+ +

There are also some assignment operators += -= +*= /=. There is not much to say: x op= +y is equivalent to x = x op y. If an exception is thrown +during the computations, the l-value is not modified (but it may be corrupt +if an exception is thrown by the base type during an assignment).

+ +

The operators / and /= will try to produce an +empty interval if the denominator is exactly zero. If the denominator +contains zero (but not only zero), the result will be the smallest interval +containing the set of division results; so one of its bound will be infinite, +but it may not be the whole interval.

+ +

lower upper median

+ +

lower, upper, median respectively +compute the lower bound, the upper bound, and the median number of an +interval ((lower+upper)/2 rounded to nearest). +width computes the width of an interval +(upper-lower rounded toward plus infinity).

+ +

min max abs square +pow division_part?

+ +

The functions min, max and abs are +also defined. Please do not mistake them for the functions defined in the +standard library (aka a<b?a:b, a>b?a:b, +a<0?-a:a). These functions are compatible with the elementary +property of interval arithmetic. For example, max([a,b], +[c,d]) = {max(x,y) such that x in +[a,b] and y in [c,d]}. They are not +defined in the std namespace but in the boost namespace in order +to avoid conflict with the other definitions.

+ +

The square function is quite particular. As you can expect +from its name, it computes the square of its argument. The reason this +function is provided is: square(x) is not x*x but +only a subset when x contains zero. For example, [-2,2]*[-2,2] = +[-4,4] but [-2,2]˛ = [0,4]; the result is a smaller interval. Consequently, +square(x) should be used instead of x*x because of +its better accuracy and a small performance improvement.

+ +

As for square, pow provides an efficient and +more accurate way to compute the integer power of an interval. Please note: +when the power is 0 and the interval is not empty, the result is 1, even if +the input interval contains 0.

+ +

The functions division_part1 and division_part2 +are useful when the user expects the division to return disjoint intervals if +necessary. For example, the narrowest closed set containg [2,3] / [-2,1] is +not ]-∞,∞[ but the union of ]-∞,-1] and [2,∞[. +When the result of the division is representable by only one interval, +division_part1 returns this interval and sets the boolean +reference to false. However, if the result needs two intervals, +division_part1 returns the negative part and sets the boolean +reference to true; a call to division_part2 is now +needed to get the positive part.

+ +

intersect hull overlap +in in_zero subset +proper_subset empty singleton +equal

+ +

intersect computes the set intersection of two closed sets, +hull computes the smallest interval which contains the two +parameters; those parameters can be numbers or intervals. If one of the +arguments is an invalid number or an empty interval, the function will only +use the other argument to compute the resulting interval (if allowed by the +checking policy).

+ +

There is no union function since the union of two intervals is not an +interval if they do not overlap. If they overlap, the hull +function computes the union.

+ +

The function overlap tests if two intervals have some common +subset. in tests if a number is in an interval; +in_zero is a variant which tests if zero is in the interval. +subset tests if the first interval is a subset of the second; +and proper_subset tests if it is a proper subset. +empty and singleton test if an interval is empty or +is a singleton. Finally, equal tests if two intervals are +equal.

+ +

sqrt log exp sin +cos tan asin acos +atan sinh cosh tanh +asinh acosh atanh +fmod

+ +

The functions sqrt, log, exp, +sin, cos, tan, asin, +acos, atan, sinh, cosh, +tanh, asinh, acosh, atanh +are also defined. There is not much to say; these functions extend the +traditional functions to the intervals and respect the basic property of +interval arithmetic. They use the checking policy to +produce empty intervals when the input interval is strictly outside of the +domain of the function.

+ +

The function fmod(interval x, interval y) expects the lower +bound of y to be strictly positive and returns an interval +z such as 0 <= z.lower() < y.upper() and such +as z is a superset of x-n*y (with n +being an integer). So, if the two arguments are positive singletons, this +function fmod(interval, interval) will behave like the +traditional function fmod(double, double).

+ +

Please note that fmod does not respect the inclusion property +of arithmetic interval. For example, the result of +fmod([13,17],[7,8]) should be [0,8] (since it must contain [0,3] +and [5,8]). But this answer is not really useful when the purpose is to +restrict an interval in order to compute a periodic function. It is the +reason why fmod will answer [5,10].

+ +

Constants

+ +

Some constants are hidden in the boost::interval_lib +namespace. They need to be explicitely templated by the interval type. The +functions are pi<I>(), pi_half<I>() and +pi_twice<I>(), and they return an object of interval type +I. Their respective values are π, π/2 and +2π.

+ +

Exception throwing

+ +

The interval class and all the functions defined around this class never +throw any exceptions by themselves. However, it does not mean that an +operation will never throw an exception. For example, let's consider the copy +constructor. As explained before, it is the default copy constructor +generated by the compiler. So it will not throw an exception if the copy +constructor of the base type does not throw an exception.

+ +

The same situation applies to all the functions: exceptions will only be +thrown if the base type or one of the three policies throws an exception.

+ +

Interval Support Library

+ +

The interval support library consists of a collection of classes that can +be used and combined to fabricate almost various commonly-needed interval +policies. In contrast to the basic classes and functions which are used in +conjunction with interval<T> (and the default policies as +the implicit second template parameter in this type), which belong simply to +the namespace boost, these components belong to the namespace +boost::numeric::interval_lib.

+ +

We merely give the synopsis here and defer each section to a separate web +page since it is only intended for the advanced user. This allows to expand +on each topic with examples, without unduly stretching the limits of this +document.

+ +

Synopsis

+
namespace boost {
+namespace numeric {
+namespace interval_lib {
+
+/* built-in rounding policy and its specializations */
+template <class T>  struct rounded_math;
+template <>         struct rounded_math<float>;
+template <>         struct rounded_math<double>;
+template <>         struct rounded_math<long double>;
+
+/* built-in rounding construction blocks */
+template <class T>  struct rounding_control;
+
+template <class T, class Rounding = rounding_control<T> >  struct rounded_arith_exact;
+template <class T, class Rounding = rounding_control<T> >  struct rounded_arith_std;
+template <class T, class Rounding = rounding_control<T> >  struct rounded_arith_opp;
+
+template <class T, class Rounding>  struct rounded_transc_dummy;
+template <class T, class Rounding = rounded_arith_exact<T> >  struct rounded_transc_exact;
+template <class T, class Rounding = rounded_arith_std<T> >  struct rounded_transc_std;
+template <class T, class Rounding = rounded_arith_opp<T> >  struct rounded_transc_opp;
+
+template <class Rounding> struct save_state;
+template <class Rounding> struct save_state_nothing;
+
+/* built-in checking policies */
+template <class T> struct checking_nothing;
+template <class T> struct checking_lax;
+template <class T> struct checking_strict;
+
+/* some metaprogramming to manipulate interval policies */
+template <class Rounding, class Checking> struct policies;
+template <class OldInterval, class NewRounding> struct change_rounding;
+template <class OldInterval, class NewChecking> struct change_checking;
+template <class OldInterval> struct unprotect;
+
+/* constants, need to be explicitly templated */
+template<class I> I pi();
+template<class I> I pi_half();
+template<class I> I pi_twice();
+
+/* interval explicit comparison functions:
+ * the mode can be cer=certainly or pos=possibly,
+ * the function lt=less_than, gt=greater_than, le=less_than_or_equal_to, ge=greater_than_or_equal_to
+ *   eq=equal_to, ne= not_equal_to */
+template <class T, class Traits>  bool cerlt(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool cerlt(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool cerlt(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool cerle(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool cerle(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool cerle(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool cergt(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool cergt(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool cergt(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool cerge(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool cerge(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool cerge(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool cereq(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool cereq(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool cereq(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool cerne(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool cerne(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool cerne(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool poslt(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool poslt(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool poslt(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool posle(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool posle(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool posle(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool posgt(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool posgt(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool posgt(const T& x, const interval<T, Traits> & y);
+
+template <class T, class Traits>  bool posge(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool posge(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool posge(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool poseq(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool poseq(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool poseq(const T& x, const interval<T, Traits>& y);
+
+template <class T, class Traits>  bool posne(const interval<T, Traits>& x, const interval<T, Traits>& y);
+template <class T, class Traits>  bool posne(const interval<T, Traits>& x, const T& y);
+template <class T, class Traits>  bool posne(const T& x, const interval<T, Traits>& y);
+
+/* comparison namespaces */
+namespace compare {
+  namespace lexicographic;
+  namespace set;
+  namespace tribool;
+} // namespace compare
+
+} // namespace interval_lib
+} // namespace numeric
+} // namespace boost
+ +

Each component of the interval support library is detailed in its own +page.

+ + +

Compilation Notes

+ +

Compilation options

+We have made many efforts to reduce the number of compilation options to a +minimum, but depending on your platform, you may find that you need one or +two. Some versions of gcc miss proper specialization of +numeric_limits and using the flag -DBOOST_NO_LIMITS +may help. Also removing the warnings for older gcc compilers about the lack +of long double support can be done with the option +-Wno-long-double. The tests and examples have successfully been +compiled with gcc 2.95 on x86 (Debian, FreeBSD), on sparc (Debian); with gcc +2.95 and stlport on x86 (Solaris), sparc (Solaris), and powerpc (Debian); +with gcc 2.96, 3.04, 3.1.1, and 3.2 on x86 (Debian); with Intel cc 6.0 and +7.0 on x86 (Debian); and with SunPRO 5.3 (Solaris). There are some issues +with Borland on x86 (Windows), and Visual C++ which we're working on. + +

Common Pitfalls and Dangers

+ +

Comparisons

+ +

One of the biggest problems is problably the correct use of the comparison +functions and operators. First, functions and operators do not try to know if +two intervals are the same mathematical object. So, if the comparison used is +"certain", then x != x is always true unless x is a +singleton interval; and the same problem arises with cereq and +cerne.

+ +

Another misleading interpretation of the comparison is: you cannot always +expect [a,b] < [c,d] to be !([a,b] >= [c,d]) since the comparison is +not necessarily total. Equality and less comparison should be seen as two +distincts relational operators.

+ +

Interval values and references

+ +

This problem is a corollary of the previous problem with x != +x. All the functions of the library only consider the value of an +interval and not the reference of an interval. In particular, you should not +expect (unless x is a singleton) the following values to be +equal: x/x and 1, x*x and square(x), +x-x and 0, etc. So the main cause of wide intervals is that +interval arithmetic does not identify different occurences of the same +variable. So, whenever possible, the user has to rewrite the formulas to +eliminate multiple occurences of the same variable. For example, +square(x)-2*x is far less precise than +square(x-1)-1.

+ +

Unprotected rounding

+ +

As explained in this section, a good way +to speed up computations when the base type is a basic floating-point type is +to unprotect the intervals at the hot spots of the algorithm. This method is +safe and really an improvement for interval computations. But please remember +that any basic floating-point operation executed inside the unprotection +blocks will probably have an undefined behavior (but only for the current +thread). And do not forget to create a rounding object as explained in the example.

+ +

Rationale

+ +

The purpose of this library is to provide an efficient and generalized way +to deal with interval arithmetic through the use of a templatized class +boost::interval. The big contention for which we provide a +rationale is the format of this class template.

+ +

It would have been easier to provide a class interval whose base type is +double. Or to follow std::complex and allow only specializations +for float, double, and long double. We +decided not to do this to allow intervals on custom types, e.g. +fixed-precision bigfloat library types (MPFR, etc.)

+ +

Policy design. Although it was tempting to make it a +class template with only one template argument, the diversity of uses for an +interval arithmetic practically forced us to use policies. The behavior of +this class can be fixed by three policies. These policies are packaged into a +single policy class, rather than making interval with three +template parameters. This is both for ease of use (the policy class can be +picked by default) and for readability.

+ +

The first policy provides all the mathematical functions on the base type +needed to define the functions on the interval type. The second one sets the +way exceptional cases encountered during computations are handled.

+ +

We could foresee situations where any combination of these policies would +be appropriate. Moreover, we wanted to enable the user of the library to +reuse the interval class template while at the same time +choosing his own behavior. See this page for some +examples.

+ +

Rounding policy. The library provides specialized +implementations of the second policy for the primitive types float and +double. In order for these implementations to be correct and fast, the +library needs to work a lot with rounding modes. Some processors are directly +dealt with and some mecanisms are provided in order to speed up the +computations. It seems to be heavy and hazardous optimizations for a gain of +only a few computer cycles; but in reality, the speed-up factor can easily go +past 2 or 3 depending on the computer. Moreover, these optimizations do not +impact the interface in any major way (with the design we have chosen, +everything can be added by specialization or by passing different template +parameters).

+ +

Pred/succ. In a previous version, two functions +pred and succ, with various corollaries like +widen, were supplied. The intent was to enlarge the interval by +one ulp (as little as possible), e.g. to ensure the inclusion property. Since +making interval a template of T, we could not define ulp for a random +parameter. In turn, rounding policies let us eliminate entirely the use of +ulp while making the intervals tighter (if a result is a representable +singleton, there is no use to widen the interval). We decided to drop those +functions.

+ +

Specialization of std::less. Since the +operator < depends on the comparison namespace chosen locally +chosen by the user, it is not possible to correctly specialize +std::less. So you have to explicitely provide such a class to +all the algorithms and templates that could require it (for example, +std::map).

+ +

Input/output. The basic header +boost/numeric/interval.hpp does not include the I/O operators for the +class interval. Somewhat default operators can be had by including +boost/numeric/interval/io.hpp, and they make very few effort to +guarantee the inclusion property (they do so for float, +double, and long double). But they cannot guarantee +it for any base number type, especially since it is not know how the type +will be displayed (for example, there could be trouble if the type is a +symbolic type that uses the character `[` or `]'). For this, we would need an +extra i/o policy. It was not felt worth the trouble. For one thing, it's +always possible to display the bounds. For another, intervals are used +somewhat internally and there is seldom need to display them. But if this +should be needed, the header is there.

+ +

Conversion from interval<T1> to +interval<T>. There is no such conversion in the +library. Specifically, there is no way to change the base number type without +explicitly casting both the lower and upper bounds from T1 to +T2 and using a constructor. We did not want to provide this +conversion because there is no easy way to make sure the inclusion property +is guaranteed: a simple casting with loss of precision can lead to invalid or +empty intervals. Moreover it is not clear who should be responsible for the +conversion: T2 is not required to know about all the possible +types T1. We also do not want to throw an exception is there is +a loss of precision (as is done for instance in numeric_cast); +instead, the result should be correctly rounded to ensure the inclusion +property. So we decided not to provide the conversion and leave it up to the +user to do it himself or herself.

+ +

History and Acknowledgments

+ +

This library was mostly inspired by previous work from Jens Maurer. Some +discussions about his work are reproduced here and the +work itself can be found here. +Jeremy Siek and Maarten Keijzer provided some rounding control for MSVC and +Sparc platforms.

+ +

Guillaume Melquiond, Hervé Brönnimann and Sylvain Pion started from the +library left by Jens and added the policy design. Guillaume and Sylvain +worked hard on the code, especially the porting and mostly tuning of the +rounding modes to the different architectures. Guillaume mostly did most of +the coding, while Sylvain and Hervé have provided some useful comments in +order for this library to be written. Hervé reorganized and wrote chapters of +the documentation based on Guillaume's great starting point.

+
+ +

Revised: 2003-01-21
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University, 2002.

+ + diff --git a/doc/numbers.htm b/doc/numbers.htm new file mode 100644 index 0000000..8024a24 --- /dev/null +++ b/doc/numbers.htm @@ -0,0 +1,129 @@ + + + + + + Numbers Requirements + + + +

Numbers Requirements

+ +

What we call "number" is the base type of the interval class. +The interval library expect a lot of properties from this base type in order +to respect the inclusion property. All these properties are already detailed +in the other sections of this documentation; but we will try to summarize +them here.

+ +

Ordering

+ +

The numbers need to be supplied with an ordering. This ordering expresses +itself by the operators < <= => > == !=. It must be +a total order (reflexivity, antisymmetry, transitivity, and each pair of +numbers is ordered). So complex<T> will not be a good +candidate for the base type; if you need the inclusion property of interval +property, you should use complex<interval<T> > in +place of interval<complex<T> > (but unfortunately, +complex only provides specialization).

+ +

Please note that invalid numbers are not concerned by the order; it can +even be conceptually better if a comparison with these invalid numbers is +always false (except for !=). If your checking +policy uses interval_lib::checking_base and your base type +contains invalid numbers, then this property is needed: nan!=nan +(here nan is an invalid number). If this property is not +present, then you should not use checking_base directly.

+ +

Numeric limits

+ +

Another remark about the checking policy. It normally is powerful enough +to handle the exceptional behavior that the basic type could induce; in +particular infinite and invalid numbers (thanks to the three functions +inf, nan and is_nan). However, if you +use interval_lib::checking_base (and the default checking policy +uses it), your base type should have a correctly specialized +std::numeric_limits<T>. In particular, the values +has_infinity and has_quiet_NaN, and the functions +infinity and quiet_NaN should be accordingly +defined.

+ +

So, to summarize, if you do not rely on the default policy and do not use +interval_lib::checking_base, it is not necessary to have a +specialization of the numeric limits for your base type.

+ +

Mathematical properties

+ +

Ensuring the numbers are correctly ordered is not enough. The basic +operators should also respect some properties depending on the order. Here +they are:

+ + +

The previous properties are also used (and enough) for abs, +square and pow. For all the transcendental +functions (including sqrt), other properties are needed. These +functions should have the same properties than the corresponding real +functions. For example, the expected properties for cos are:

+ + +

Rounding

+ +

If you work with a base type and no inexact result is ever computed, you +can skip the rest of this paragraph. You can also skip it if you are not +interested in the inclusion property (if approximate results are enough). If +you are still reading, it is probably because you want to know the basic +properties the rounding policy should validate.

+ +

Whichever operation or function you consider, the following property +should be present: f_down(x,y) <= f(x,y) <= f_up(x,y). +Here, f denotes the infinitely precise function computed and +f_down and f_up are functions which return possibly +inexact values but of the correct type (the base type). If possible, they +should try to return the nearest representable value, but it is not always +easy.

+ +

Constants

+ +

In order for the trigonometric functions to correctly work, the library +need to know the value of the π constant (and also π/2 and +2π). Since these constants may not be representable in the base type, +the library does not have to know the exact value: a lower bound and an upper +bound are enough. If these values are not provided by the user, the default +values will be used: they are integer values (so π is bounded by 3 and +4).

+ +

Operators and conversions

+ +

As explained at the beginning, the comparison operators should be defined +for the base type. The rounding policy defines a lot of functions used by the +interval library. So the arithmetic operators do not need to be defined for +the base type (unless required by one of the predefined classes). However, +there is an exception: the unary minus need to be defined. Moreover, this +operator should only provide exact results; it is the reason why the rounding +policy does not provide some negation functions.

+ +

The conversion from int to the base type need to be defined. +The conversion the other way around is provided by the rounding policy; and +no other conversion is needed.

+
+ +

Revised: 2003-01-21
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/policies.htm b/doc/policies.htm new file mode 100644 index 0000000..684a727 --- /dev/null +++ b/doc/policies.htm @@ -0,0 +1,64 @@ + + + + + + Policies + + + +

Policies

+ +

The interval template requires two arguments. The first +corresponds to the base type chosen for the bounds. And the second defines +the rounding and checking behaviors of the newly constructed class. +This second argument is not mandatory but may need some customizations. In +order to ease the manipulations, some helper templates are provided in +interval/policies.hpp.

+
namespace boost {
+namespace numeric {
+namespace interval_lib {
+
+template<class Rounding, class Checking>
+struct policies {
+  typedef Rounding rounding;
+  typedef Checking checking;
+};
+
+template<class OldInterval, class NewRounding>
+struct change_rounding {
+  typedef ... type;
+};
+
+template<class OldInterval, class NewChecking>
+struct change_checking {
+  typedef ... type;
+};
+
+template<class OldInterval>
+struct unprotect {
+  typedef ... type;
+};
+
+} // namespace interval_lib
+} // namespace numeric
+} // namespace boost
+ +

The policies template should be used whenever the user needs +to define a policy structure for an interval class. +change_rounding and change_checking can be used to +get the type of a new interval by changing one of the policies of an old +interval; the new type is available thanks to the type definition +type. Finally, unprotect looks like +change_rounding and directly changes the rounding of an interval +to its unprotected version (a better explanation is available here).

+
+ +

Revised: 2003-01-21
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/rounding.htm b/doc/rounding.htm new file mode 100644 index 0000000..e2ff70e --- /dev/null +++ b/doc/rounding.htm @@ -0,0 +1,559 @@ + + + + + + Rounding Policies + + + +

Rounding Policies

+ +

In order to be as general as possible, the library uses a class to compute +all the necessary functions rounded upward or downward. This class is the +first parameter of policies, it is also the type named +rounding in the policy definition of interval.

+ +

By default, it is interval_lib::rounded_math<T>. +The class interval_lib::rounded_math is already +specialized for the standard floating types (float , +double and long double). So if the base type of +your intervals is not one of these, a good solution would probably be to +provide a specialization of this class. But if the default specialization of +rounded_math<T> for float, +double, or long double is not what you seek, or you +do not want to specialize +interval_lib::rounded_math<T> (say because you +prefer to work in your own namespace) you can also define your own rounding +policy and pass it directly to interval_lib::policies.

+ +

Requirements

+ +

Here comes what the class is supposed to provide. The domains are written +next to their respective functions (as you can see, the functions do not have +to worry about invalid values, but they have to handle infinite +arguments).

+
/* Rounding requirements */
+struct rounding {
+  // defaut constructor, destructor
+  rounding();
+  ~rounding();
+  // mathematical operations
+  T add_down(T, T); // [-∞;+∞][-∞;+∞]
+  T add_up  (T, T); // [-∞;+∞][-∞;+∞]
+  T sub_down(T, T); // [-∞;+∞][-∞;+∞]
+  T sub_up  (T, T); // [-∞;+∞][-∞;+∞]
+  T mul_down(T, T); // [-∞;+∞][-∞;+∞]
+  T mul_up  (T, T); // [-∞;+∞][-∞;+∞]
+  T div_down(T, T); // [-∞;+∞]([-∞;+∞]-{0})
+  T div_up  (T, T); // [-∞;+∞]([-∞;+∞]-{0})
+  T sqrt_down(T);   // ]0;+∞]
+  T sqrt_up  (T);   // ]0;+∞]
+  T exp_down(T);    // [-∞;+∞]
+  T exp_up  (T);    // [-∞;+∞]
+  T log_down(T);    // ]0;+∞]
+  T log_up  (T);    // ]0;+∞]
+  T cos_down(T);    // [0;2π]
+  T cos_up  (T);    // [0;2π]
+  T tan_down(T);    // ]-π/2;π/2[
+  T tan_up  (T);    // ]-π/2;π/2[
+  T asin_down(T);   // [-1;1]
+  T asin_up  (T);   // [-1;1]
+  T acos_down(T);   // [-1;1]
+  T acos_up  (T);   // [-1;1]
+  T atan_down(T);   // [-∞;+∞]
+  T atan_up  (T);   // [-∞;+∞]
+  T sinh_down(T);   // [-∞;+∞]
+  T sinh_up  (T);   // [-∞;+∞]
+  T cosh_down(T);   // [-∞;+∞]
+  T cosh_up  (T);   // [-∞;+∞]
+  T tanh_down(T);   // [-∞;+∞]
+  T tanh_up  (T);   // [-∞;+∞]
+  T asinh_down(T);  // [-∞;+∞]
+  T asinh_up  (T);  // [-∞;+∞]
+  T acosh_down(T);  // [1;+∞]
+  T acosh_up  (T);  // [1;+∞]
+  T atanh_down(T);  // [-1;1]
+  T atanh_up  (T);  // [-1;1] 
+  T median(T, T);   // [-∞;+∞][-∞;+∞]
+  T int_down(T);    // [-∞;+∞]
+  T int_up  (T);    // [-∞;+∞]
+  // unprotected rounding class
+  typedef ... unprotected_rounding;
+};
+ +

The constructor and destructor of the rounding class have a very important +semantic requirement: they are responsible for setting and resetting the +rounding modes of the computation on T. For instance, if T is a standard +floating point type and floating point computation is performed according to +the Standard IEEE 754, the constructor can save the current rounding state, +each _up (resp. _down) function will round up +(resp. down), and the destructor will restore the saved rounding state. +Indeed this is the behavior of the default rounding policy.

+ +

The meaning of all the mathematical functions up until +atanh_up is clear: each function returns number representable in +the type T which is a lower bound (for _down) or +upper bound (for _up) on the true mathematical result of the +corresponding function. The function median computes the average +of its two arguments rounded to its nearest representable number. The +functions int_down and int_up compute the nearest +integer smaller or bigger than their argument.

+ +

The type unprotected_rounding allows to remove all controls. +For reasons why one might to do this, see the protection paragraph below.

+ +

Overview of the provided classes

+ +

A lot of classes are provided. The classes are organized by level. At the +bottom is the class rounding_control. At the next level come +rounded_arith_exact, rounded_arith_std and +rounded_arith_opp. Then there are +rounded_transc_dummy, rounded_transc_exact, +rounded_transc_std and rounded_transc_opp. And +finally are save_state and save_state_nothing. Each +of these classes provide a set of members that are required by the classes of +the next level. For example, a rounded_transc_... class needs +the members of a rounded_arith_... class.

+ +

When they exist in two versions _std and _opp, +the first one does switch the rounding mode each time, and the second one +tries to keep it oriented toward plus infinity. The main purpose of the +_opp version is to speed up the computations through the use of +the "opposite trick" (see the performance notes). This +version requires the rounding mode to be upward before entering any +computation functions of the class. It guarantees that the rounding mode will +still be upward at the exit of the functions.

+ +

Please note that it is really a very bad idea to mix the _opp +version with the _std since they do not have compatible +properties.

+ +

There is a third version named _exact which computes the +functions without changing the rounding mode. It is an "exact" version +because it is intended for a base type that produces exact results.

+ +

The last version is the _dummy version. It does not do any +computations but still produces compatible results.

+ +

Please note that it is possible to use the "exact" version for an inexact +base type, e.g. float or double. In that case, the +inclusion property is no longer guaranteed, but this can be useful to speed +up the computation when the inclusion property is not desired strictly. For +instance, in computer graphics, a small error due to floating-point roundoff +is acceptable as long as an approximate version of the inclusion property +holds.

+ +

Here comes what each class defines. Later, when they will be described +more thoroughly, these members will not be repeated. Please come back here in +order to see them. Inheritance is also used to avoid repetitions.

+
template <class T>
+struct rounding_control
+{
+  typedef ... rounding_mode;
+  void set_rounding_mode(rounding_mode);
+  void get_rounding_mode(rounding_mode&);
+  void downward ();
+  void upward   ();
+  void to_nearest();
+  T to_int(T);
+  T force_rounding(T);
+};
+
+template <class T, class Rounding>
+struct rounded_arith_... : Rounding
+{
+  void init();
+  T add_down(T, T);
+  T add_up  (T, T);
+  T sub_down(T, T);
+  T sub_up  (T, T);
+  T mul_down(T, T);
+  T mul_up  (T, T);
+  T div_down(T, T);
+  T div_up  (T, T);
+  T sqrt_down(T);
+  T sqrt_up  (T);
+  T median(T, T);
+  T int_down(T);
+  T int_up  (T);
+};
+
+template <class T, class Rounding>
+struct rounded_transc_... : Rounding
+{
+  T exp_down(T);
+  T exp_up  (T);
+  T log_down(T);
+  T log_up  (T);
+  T cos_down(T);
+  T cos_up  (T);
+  T tan_down(T);
+  T tan_up  (T);
+  T asin_down(T);
+  T asin_up  (T);
+  T acos_down(T);
+  T acos_up  (T);
+  T atan_down(T);
+  T atan_up  (T);
+  T sinh_down(T);
+  T sinh_up  (T);
+  T cosh_down(T);
+  T cosh_up  (T);
+  T tanh_down(T);
+  T tanh_up  (T);
+  T asinh_down(T);
+  T asinh_up  (T);
+  T acosh_down(T);
+  T acosh_up  (T);
+  T atanh_down(T);
+  T atanh_up  (T);
+};
+
+template <class Rounding>
+struct save_state_... : Rounding
+{
+  save_state_...();
+  ~save_state_...();
+  typedef ... unprotected_rounding;
+};
+ +

Synopsis.

+
namespace boost {
+namespace numeric {
+namespace interval_lib {
+
+/* basic rounding control */
+template <class T>  struct rounding_control;
+
+/* arithmetic functions rounding */
+template <class T, class Rounding = rounding_control<T> > struct rounded_arith_exact;
+template <class T, class Rounding = rounding_control<T> > struct rounded_arith_std;
+template <class T, class Rounding = rounding_control<T> > struct rounded_arith_opp;
+
+/* transcendental functions rounding */
+template <class T, class Rounding> struct rounded_transc_dummy;
+template <class T, class Rounding = rounded_arith_exact<T> > struct rounded_transc_exact;
+template <class T, class Rounding = rounded_arith_std<T> > struct rounded_transc_std;
+template <class T, class Rounding = rounded_arith_opp<T> > struct rounded_transc_opp;
+
+/* rounding-state-saving classes */
+template <class Rounding> struct save_state;
+template <class Rounding> struct save_state_nothing;
+
+/* default policy for type T */
+template <class T>  struct rounded_math;
+template <>  struct rounded_math<float>;
+template <>  struct rounded_math<double>;
+
+/* some metaprogramming to convert a protected to unprotected rounding */
+template <class I> struct unprotect;
+
+} // namespace interval_lib
+} // namespace numeric
+} // namespace boost
+ +

Description of the provided classes

+ +

We now describe each class in the order they appear in the definition of a +rounding policy (this outermost-to-innermost order is the reverse order from +the synopsis).

+ +

Protection control

+ +

Protection refers to the fact that the interval operations will be +surrounded by rounding mode controls. Unprotecting a class means to remove +all the rounding controls. Each rounding policy provides a type +unprotected_rounding. The required type +unprotected_rounding gives another rounding class that enables +to work when nested inside rounding. For example, the first three lines below +should all produce the same result (because the first operation is the +rounding constructor, and the last is its destructor, which take care of +setting the rounding modes); and the last line is allowed to have an +undefined behavior (since no rounding constructor or destructor is ever +called).

+
T c; { rounding rnd; c = rnd.add_down(a, b); }
+T c; { rounding rnd1; { rounding rnd2; c = rnd2.add_down(a, b); } }
+T c; { rounding rnd1; { rounding::unprotected_rounding rnd2; c = rnd2.add_down(a, b); } }
+T d; { rounding::unprotected_rounding rnd; d = rnd.add_down(a, b); }
+ +

Naturally rounding::unprotected_rounding may simply be +rounding itself. But it can improve performance if it is a +simplified version with empty constructor and destructor. In order to avoid +undefined behaviors, in the library, an object of type +rounding::unprotected_rounding is guaranteed to be created only +when an object of type rounding is already alive. See the performance notes for some additional details.

+ +

The support library defines a metaprogramming class template +unprotect which takes an interval type I and +returns an interval type unprotect<I>::type where the +rounding policy has been unprotected. Some information about the types: +interval<T, interval_lib::policies<Rounding, _> +>::traits_type::rounding is the same type as +Rounding, and unprotect<interval<T, +interval_lib::policies<Rounding, _> > >::type is +the same type as interval<T, +interval_lib::policies<Rounding::unprotected, _> >.

+ +

State saving

+ +

First comes save_state. This class is responsible for saving +the current rounding mode and calling init in its constructor, and for +restoring the saved rounding mode in its destructor. This class also defines +the unprotected_rounding type.

+ +

If the rounding mode does not require any state-saving or initialization, +save_state_nothing can be used instead of +save_state.

+ +

Transcendental functions

+ +

The classes rounded_transc_exact, +rounded_transc_std and rounded_transc_opp expect +the std namespace to provide the functions exp log cos tan acos asin atan +cosh sinh tanh acosh asinh atanh. For the _std and +_opp versions, all these functions should respect the current +rounding mode fixed by a call to downward or upward.

+ +

Please note: Unfortunately, the latter is rarely the +case. It is the reason why a class rounded_transc_dummy is +provided which does not depend on the functions from the std namespace. There +is no magic, however. The functions of rounded_transc_dummy do +not compute anything. They only return valid values. For example, +cos_down always returns -1. In this way, we do verify the +inclusion property for the default implementation, even if this has strictly +no value for the user. In order to have useful values, another policy should +be used explicitely, which will most likely lead to a violation of the +inclusion property. In this way, we ensure that the violation is clearly +pointed out to the user who then knows what he stands against. This class +could have been used as the default transcendental rounding class, but it was +decided it would be better for the compilation to fail due to missing +declarations rather than succeed thanks to valid but unusable functions.

+ +

Basic arithmetic functions

+ +

The classes rounded_arith_std and +rounded_arith_opp expect the operators + - * / and the function +std::sqrt to respect the current rounding mode.

+ +

The class rounded_arith_exact requires +std::floor and std::ceil to be defined since it can +not rely on to_int.

+ +

Rounding control

+ +

The functions defined by each of the previous classes did not need any +explanation. For example, the behavior of add_down is to compute +the sum of two numbers rounded downward. For rounding_control, +the situation is a bit more complex.

+ +

The basic function is force_rounding which returns its +argument correctly rounded accordingly to the current rounding mode if it was +not already the case. This function is necessary to handle delayed rounding. +Indeed, depending on the way the computations are done, the intermediate +results may be internaly stored in a more precise format and it can lead to a +wrong rounding. So the function enforces the rounding. Here is an example of what happens when the rounding is +not enforced.

+ +

The function get_rounding_mode returns the current rounding +mode, set_rounding_mode sets the rounding mode back to a +previous value returned by get_rounding_mode. +downward, upward and to_nearest sets +the rounding mode in one of the three directions. This rounding mode should +be global to all the functions that use the type T. For example, +after a call to downward, force_rounding(x+y) is +expected to return the sum rounded toward -∞.

+ +

The function to_int computes the nearest integer accordingly +to the current rounding mode.

+ +

The non-specialized version of rounding_control does not do +anything. The functions for the rounding mode are empty, and +to_int and force_rounding are identity functions. +The pi_ constant functions return suitable integers (for +example, pi_up returns T(4)).

+ +

The class template rounding_control is specialized for +float, double and long double in order +to best use the floating point unit of the computer.

+ +

Template class rounded_math

+ +

The default policy (aka rounded_math<T>) is simply +defined as:

+
template <class T> struct rounded_math<T> : save_state_nothing<rounded_arith_exact<T> > {};
+ +

and the specializations for float, double and +long double use rounded_arith_opp, as in:

+
template <> struct rounded_math<float>       : save_state<rounded_arith_opp<float> >       {};
+template <> struct rounded_math<double>      : save_state<rounded_arith_opp<double> >      {};
+template <> struct rounded_math<long double> : save_state<rounded_arith_opp<long double> > {};
+ +

Performance Issues

+ +

This paragraph deals mostly with the performance of the library with +intervals using the floating-point unit (FPU) of the computer. Let's consider +the sum of [a,b] and [c,d] as an example. The +result is [down(a+c), +up(b+d)], where down and +up indicate the rounding mode needed.

+ +

Rounding Mode Switch

+ +

If the FPU is able to use a different rounding mode for each operation, +there is no problem. For example, it's the case for the Alpha processor: each +floating-point instruction can specify a different rounding mode. However, +the IEEE-754 Standard does not require such a behavior. So most of the FPUs +only provide some instructions to set the rounding mode for all subsequent +operations. And generally, these instructions need to flush the pipeline of +the FPU.

+ +

In this situation, the time needed to sum [a,b] and +[c,d] is far worse than the time needed to calculate +a+b and c+d since the two additions cannot be +parallelized. Consequently, the objective is to diminish the number of +rounding mode switches.

+ +

If this library is not used to provide exact computations, but only for +pair arithmetic, the solution is quite simple: do not use rounding. In that +case, doing the sum [a,b] and [c,d] will be as +fast as computing a+b and c+d. Everything is +perfect.

+ +

However, if exact computations are required, such a solution is totally +unthinkable. So, are we penniless? No, there is still a trick available. +Indeed, down(a+c) = -up(-a-c) if the unary minus +is an exact operation. It is now possible to calculate the whole sum with the +same rounding mode. Generally, the cost of the mode switching is worse than +the cost of the sign changes.

+ +

Let's recapitulate. Before, when doing an addition, there were three +rounding mode switches (down, up and restore). Now, with this little trick, +there are only two switches (up and restore). It is better, but still a +bottleneck when many operations are nested. Indeed, the generated code for +[a,b] + [c,d] + [e,f] will probably +look like:

+
up();
+t1 = -(-a - c);
+t2 = b + d;
+restore();
+up();
+x = -(-t1 - e);
+y = t2 + f;
+restore();
+ +

If you think it is possible to do much better, you are right. For example, +this is better (and probably optimal):

+
up();
+x = -(-a - c - e);
+y = b + d + f;
+restore();
+ +

Such a code will be generated by a compiler if the computations are made +without initialization and restoration of the rounding mode. However, it +would be far too easy if there were no drawback: because the rounding mode is +not restored in the meantime, operations on floating-point numbers must be +prohibited. This method can only be used if all the operations are operations +on intervals (or operations between an interval and a floating point +number).

+ +

Example

+ +

Here is an example of the Horner scheme to compute the value of a polynom. +The rounding mode switches are disabled for the whole computation.

+
// I is an interval class, the polynom is a simple array
+template<class I>
+I horner(const I& x, const I p[], int n) {
+
+  // save and initialize the rounding mode
+  typename I::traits_type::rounding rnd;
+
+  // define the unprotected version of the interval type
+  typedef typename boost::numeric::interval_lib::unprotect<I>::type R;
+  
+  const R& a = x;
+  R y = p[n - 1];
+  for(int i = n - 2; i >= 0; i--) {
+    y = y * a + (const R&)(p[i]);
+  }
+  return y;
+
+  // restore the rounding mode with the destruction of rnd
+}
+ +

Please note that a rounding object is especially created in order to +compensate for the protection loss. Each interval of type I is converted in +an interval of type R before any operations. If this conversion is not done, +the result is still correct, but the interest of this whole optimization has +disappeared. Whenever possible, it is good to convert to const +R& instead of R: indeed, the function could already +be called inside an unprotection block so the types R and +I would be the same interval, no need for a conversion.

+ +

Uninteresting remark

+ +

It was said at the beginning that the Alpha processors can use a specific +rounding mode for each operation. However, due to the instruction format, the +rounding toward plus infinity is not available. Only the rounding toward +minus infinity can be used. So the trick using the change of sign becomes +essential, but there is no need to save and restore the rounding mode on both +sides of an operation.

+ +

Extended Registers

+ +

There is another problem besides the cost of the rounding mode switch. +Some FPUs use extended registers (for example, float computations will be +done with double registers, or double computations with long double +registers). Consequently, many problems can arise.

+ +

The first one is due to to the extended precision of the mantissa. The +rounding is also done on this extended precision. And consequently, we still +have down(a+b) = -up(-a-b) in the extended +registers. But back to the standard precision, we now have +down(a+b) < -up(-a-b) instead of an equality. +A solution could be not to use this method. But there still are other +problems, with the comparisons between numbers for example.

+ +

Naturally, there is also a problem with the extended precision of the +exponent. To illustrate this problem, let m be the biggest number +before +inf. If we calculate 2*[m,m], the answer should +be [m,inf]. But due to the extended registers, the FPU will +first store [2m,2m] and then convert it to +[inf,inf] at the end of the calculus (when the rounding mode is +toward +inf). So the answer is no more accurate.

+ +

There is only one solution: to force the FPU to convert the extended +values back to standard precision after each operation. Some FPUs provide an +instruction able to do this conversion (for example the PowerPC processors). +But for the FPUs that do not provide it (the x86 processors), the only +solution is to write the values to memory and read them back. Such an +operation is obviously very expensive.

+ +

Some Examples

+ +

Here come several cases:

+ +
+ +

Revised: 2003-01-21
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/doc/todo.htm b/doc/todo.htm new file mode 100644 index 0000000..32194ea --- /dev/null +++ b/doc/todo.htm @@ -0,0 +1,45 @@ + + + + + + Interval-TODO.htm + + + +

TODO list for the Interval Arithmetic library

+ +

Comments from the review process

+ + +

Various items

+ +
+ +

Revised: 2003-01-19
+Copyright (c) Guillaume Melquiond, Sylvain Pion, Hervé Brönnimann, 2002.
+Polytechnic University.

+ + diff --git a/examples/findroot_demo.cpp b/examples/findroot_demo.cpp new file mode 100644 index 0000000..0446df9 --- /dev/null +++ b/examples/findroot_demo.cpp @@ -0,0 +1,170 @@ +/* boost findroot_demo.cc find zero points of some function + * + * Copyright Jens Maurer 2000 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without free provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * Jens Maurer makes no representations about the suitability of this + * software for any purpose. It is provided "as is" without express or + * implied warranty. + * + * $Id$ + * + * + * The idea and the 2D function are based on RVInterval, + * which contains the following following copyright notice: + + This file is copyrighted 1996 by Ronald Van Iwaarden. + + Permission is hereby granted, without written agreement and + without license or royalty fees, to use, copy, modify, and + distribute this software and its documentation for any + purpose, subject to the following conditions: + + The above license notice and this permission notice shall + appear in all copies or substantial portions of this software. + + The name "RVInterval" cannot be used for any modified form of + this software that does not originate from the authors. + Nevertheless, the name "RVInterval" may and should be used to + designate the optimization software implemented and described + in this package, even if embedded in any other system, as long + as any portion of this code remains. + + The authors specifically disclaim any warranties, including, + but not limited to, the implied warranties of merchantability + and fitness for a particular purpose. The software provided + hereunder is on an "as is" basis, and the authors have no + obligation to provide maintenance, support, updates, + enhancements, or modifications. In no event shall the authors + be liable to any party for direct, indirect, special, + incidental, or consequential damages arising out of the use of + this software and its documentation. +*/ + +#include // must be first for workaround +#include +#include +#include +#include +#include +#include + + +template +struct test_func2d +{ + T operator()(T x, T y) const + { + return sin(x)*cos(y) - exp(x*y)/45.0 * (pow(x+y, 2)+100.0) - + cos(sin(y))*y/4.0; + } +}; + +template +struct test_func1d +{ + T operator()(T x) const + { + return sin(x)/(x*x+1.0); + } +}; + +template +struct test_func1d_2 +{ + T operator()(T x) const + { + using std::sqrt; + return sqrt(x*x-1.0); + } +}; + +template +void find_zeros(std::ostream & os, Function f, I searchrange) +{ + std::list l, done; + l.push_back(searchrange); + while(!l.empty()) { + I range = l.front(); + l.pop_front(); + I val = f(range); + if (in_zero(val)) { + if(width(range) < 1e-6) { + os << range << '\n'; + continue; + } + // there's still a solution hidden somewhere + std::pair p = bisect(range); + l.push_back(p.first); + l.push_back(p.second); + } + } +} + +template +inline std::ostream & operator<<(std::ostream & os, const std::pair & x) +{ + os << "(" << x.first << ", " << x.second << ")"; + return os; +} + +static const double epsilon = 5e-3; + +template +void find_zeros(std::ostream & os, Function f, I rx, I ry) +{ + typedef std::pair rectangle; + typedef std::deque container; + container l, done; + // l.reserve(50); + l.push_back(std::make_pair(rx, ry)); + for(int i = 1; !l.empty(); ++i) { + rectangle rect = l.front(); + l.pop_front(); + I val = f(rect.first, rect.second); + if (in_zero(val)) { + if(width(rect.first) < epsilon && width(rect.second) < epsilon) { + os << median(rect.first) << " " << median(rect.second) << " " + << lower(rect.first) << " " << upper(rect.first) << " " + << lower(rect.second) << " " << upper(rect.second) + << '\n'; + } else { + if(width(rect.first) > width(rect.second)) { + std::pair p = bisect(rect.first); + l.push_back(std::make_pair(p.first, rect.second)); + l.push_back(std::make_pair(p.second, rect.second)); + } else { + std::pair p = bisect(rect.second); + l.push_back(std::make_pair(rect.first, p.first)); + l.push_back(std::make_pair(rect.first, p.second)); + } + } + } + if(i % 10000 == 0) + std::cerr << "\rIteration " << i << ", l.size() = " << l.size(); + } + std::cerr << '\n'; +} + +int main() +{ + using namespace boost; + using namespace numeric; + using namespace interval_lib; + + typedef interval >, + checking_base > > I; + + std::cout << "Zero points of sin(x)/(x*x+1)\n"; + find_zeros(std::cout, test_func1d(), I(-11, 10)); + std::cout << "Zero points of sqrt(x*x-1)\n"; + find_zeros(std::cout, test_func1d_2(), I(-5, 6)); + std::cout << "Zero points of Van Iwaarden's 2D function\n"; + std::ofstream f("func2d.data"); + find_zeros(f, test_func2d(), I(-20, 20), I(-20, 20)); + std::cout << "Use gnuplot, command 'plot \"func2d.data\" with dots' to plot\n"; +} diff --git a/examples/horner.cpp b/examples/horner.cpp new file mode 100644 index 0000000..9b5c734 --- /dev/null +++ b/examples/horner.cpp @@ -0,0 +1,31 @@ +#include +#include +#include + +// I is an interval class, the polynom is a simple array +template +I horner(const I& x, const I p[], int n) { + + // initialize and restore the rounding mode + typename I::traits_type::rounding rnd; + + // define the unprotected version of the interval type + typedef typename boost::numeric::interval_lib::unprotect::type R; + + const R& a = x; + R y = p[n - 1]; + for(int i = n - 2; i >= 0; i--) { + y = y * a + (const R&)(p[i]); + } + return y; + + // restore the rounding mode with the destruction of rnd +} + +int main() { + typedef boost::numeric::interval I; + I p[3] = { -1.0, 0, 1.0 }; + I x = 1.0; + std::cout << horner(x, p, 3) << std::endl; + return 0; +} diff --git a/examples/newton-raphson.cpp b/examples/newton-raphson.cpp new file mode 100644 index 0000000..e47eb7c --- /dev/null +++ b/examples/newton-raphson.cpp @@ -0,0 +1,125 @@ +#include +#include +#include +#include +#include +#include +#include + +template I f(const I& x) +{ return x * (x - 1.) * (x - 2.) * (x - 3.) * (x - 4.); } +template I f_diff(const I& x) +{ return (((5. * x - 40.) * x + 105.) * x - 100.) * x + 24.; } + +static const double max_width = 1e-10; +static const double alpha = 0.75; + +using namespace boost; +using namespace numeric; +using namespace interval_lib; + +// First method: no empty intervals + +typedef interval I1_aux; +typedef unprotect::type I1; + +std::vector newton_raphson(const I1& xs) { + std::vector l, res; + I1 vf, vd, x, x1, x2; + l.push_back(xs); + while (!l.empty()) { + x = l.back(); + l.pop_back(); + bool x2_used; + double xx = median(x); + vf = f(xx); + vd = f_diff(x); + if (in_zero(vf) && in_zero(vd)) { + x1 = I1::whole(); + x2_used = false; + } else { + x1 = xx - division_part1(vf, vd, x2_used); + if (x2_used) x2 = xx - division_part2(vf, vd); + } + if (overlap(x1, x)) x1 = intersect(x, x1); + else if (x2_used) { x1 = x2; x2_used = false; } + else continue; + if (x2_used) + if (overlap(x2, x)) x2 = intersect(x, x2); + else x2_used = false; + if (x2_used && width(x2) > width(x1)) std::swap(x1, x2); + if (!in_zero(f(x1))) + if (x2_used) { x1 = x2; x2_used = false; } + else continue; + if (width(x1) < max_width) res.push_back(x1); + else if (width(x1) > alpha * width(x)) { + std::pair p = bisect(x); + if (in_zero(f(p.first))) l.push_back(p.first); + x2 = p.second; + x2_used = true; + } else l.push_back(x1); + if (x2_used && in_zero(f(x2))) + if (width(x2) < max_width) res.push_back(x1); + else l.push_back(x2); + } + return res; +} + +// Second method: with empty intervals + +typedef change_checking >::type I2_aux; +typedef unprotect::type I2; + +std::vector newton_raphson(const I2& xs) { + std::vector l, res; + I2 vf, vd, x, x1, x2; + l.push_back(xs); + while (!l.empty()) { + x = l.back(); + l.pop_back(); + double xx = median(x); + vf = f(xx); + vd = f_diff(x); + if (in_zero(vf) && in_zero(vd)) { + x1 = x; + x2 = I2::empty(); + } else { + bool x2_used; + x1 = intersect(x, xx - division_part1(vf, vd, x2_used)); + x2 = x2_used ? intersect(x, xx - division_part2(vf, vd)) : I2::empty(); + } + if (width(x2) > width(x1)) std::swap(x1, x2); + if (empty(x1) || !in_zero(f(x1))) + if (!empty(x2)) { x1 = x2; x2 = I2::empty(); } + else continue; + if (width(x1) < max_width) res.push_back(x1); + else if (width(x1) > alpha * width(x)) { + std::pair p = bisect(x); + if (in_zero(f(p.first))) l.push_back(p.first); + x2 = p.second; + } else l.push_back(x1); + if (!empty(x2) && in_zero(f(x2))) + if (width(x2) < max_width) res.push_back(x1); + else l.push_back(x2); + } + return res; +} + +int main() { + { + I1_aux::traits_type::rounding rnd; + std::vector res = newton_raphson(I1(-1, 5.1)); + std::cout << "Results: " << std::endl << std::setprecision(12); + for(std::vector::const_iterator i = res.begin(); i != res.end(); ++i) + std::cout << " " << *i << std::endl; + std::cout << std::endl; + } + { + I2_aux::traits_type::rounding rnd; + std::vector res = newton_raphson(I2(-1, 5.1)); + std::cout << "Results: " << std::endl << std::setprecision(12); + for(std::vector::const_iterator i = res.begin(); i != res.end(); ++i) + std::cout << " " << *i << std::endl; + std::cout << std::endl; + } +} diff --git a/examples/rational.cpp b/examples/rational.cpp new file mode 100644 index 0000000..edd0a2e --- /dev/null +++ b/examples/rational.cpp @@ -0,0 +1,37 @@ +// Example program of how to use interval<> of rational<> numbers. +// +// Sylvain Pion, Guillaume Melquiond, 2002. + +// it would have been enough to only include: +// +// but it's a bit overkill to include processor intrinsics +// and transcendental functions, so we do it by ourselves + +#include // base class +#include // default arithmetic rounding policy +#include // default checking policy +#include // += *= -= etc +#include // default policy + +#include +#include + +typedef boost::rational Rat; +typedef boost::numeric::interval Interval; + +std::ostream& operator<<(std::ostream& os, const Interval& r) { + os << "[" << r.lower() << "," << r.upper() << "]"; + return os; +} + +int main() { + Rat p(2, 3), q(3, 4); + Interval z(4, 5); + Interval a(p, q); + a += z; + z *= q; + a -= p; + a /= q; + std::cout << z << std::endl; + std::cout << a << std::endl; +} diff --git a/include/boost/numeric/interval/arith.hpp b/include/boost/numeric/interval/arith.hpp new file mode 100644 index 0000000..aef3179 --- /dev/null +++ b/include/boost/numeric/interval/arith.hpp @@ -0,0 +1,314 @@ +/* Boost interval/arith.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP +#define BOOST_NUMERIC_INTERVAL_ARITH_HPP + +#include +#include +#include +#include +#include + +namespace boost { +namespace numeric { + +/* + * Basic arithmetic operators + */ + +template inline +const interval& operator+(const interval& x) +{ + return x; +} + +template inline +interval operator-(const interval& x) +{ + if (interval_lib::detail::test_input(x)) + return interval::empty(); + return interval(-x.upper(), -x.lower(), true); +} + +template inline +interval& interval::operator+=(const interval& r) +{ + if (interval_lib::detail::test_input(*this, r)) + set_empty(); + else { + typename Policies::rounding rnd; + set(rnd.add_down(low, r.low), rnd.add_up(up, r.up)); + } + return *this; +} + +template inline +interval& interval::operator+=(const T& r) +{ + if (interval_lib::detail::test_input(*this, r)) + set_empty(); + else { + typename Policies::rounding rnd; + set(rnd.add_down(low, r), rnd.add_up(up, r)); + } + return *this; +} + +template inline +interval& interval::operator-=(const interval& r) +{ + if (interval_lib::detail::test_input(*this, r)) + set_empty(); + else { + typename Policies::rounding rnd; + set(rnd.sub_down(low, r.up), rnd.sub_up(up, r.low)); + } + return *this; +} + +template inline +interval& interval::operator-=(const T& r) +{ + if (interval_lib::detail::test_input(*this, r)) + set_empty(); + else { + typename Policies::rounding rnd; + set(rnd.sub_down(low, r), rnd.sub_up(up, r)); + } + return *this; +} + +template inline +interval& interval::operator*=(const interval& r) +{ + return *this = *this * r; +} + +template inline +interval& interval::operator*=(const T& r) +{ + return *this = r * *this; +} + +template inline +interval& interval::operator/=(const interval& r) +{ + return *this = *this / r; +} + +template inline +interval& interval::operator/=(const T& r) +{ + return *this = *this / r; +} + +template inline +interval operator+(const interval& x, + const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + return interval(rnd.add_down(x.lower(), y.lower()), + rnd.add_up (x.upper(), y.upper()), true); +} + +template inline +interval operator+(const T& x, const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + return interval(rnd.add_down(x, y.lower()), + rnd.add_up (x, y.upper()), true); +} + +template inline +interval operator+(const interval& x, const T& y) +{ return y + x; } + +template inline +interval operator-(const interval& x, + const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + return interval(rnd.sub_down(x.lower(), y.upper()), + rnd.sub_up (x.upper(), y.lower()), true); +} + +template inline +interval operator-(const T& x, const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + return interval(rnd.sub_down(x, y.upper()), + rnd.sub_up (x, y.lower()), true); +} + +template inline +interval operator-(const interval& x, const T& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + return interval(rnd.sub_down(x.lower(), y), + rnd.sub_up (x.upper(), y), true); +} + +template inline +interval operator*(const interval& x, + const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + BOOST_NUMERIC_INTERVAL_using_max(max); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + typename Policies::rounding rnd; + const T& xl = x.lower(); + const T& xu = x.upper(); + const T& yl = y.lower(); + const T& yu = y.upper(); + + if (interval_lib::detail::is_neg(xl)) + if (interval_lib::detail::is_pos(xu)) + if (interval_lib::detail::is_neg(yl)) + if (interval_lib::detail::is_pos(yu)) // M * M + return I(min(rnd.mul_down(xl, yu), rnd.mul_down(xu, yl)), + max(rnd.mul_up (xl, yl), rnd.mul_up (xu, yu)), true); + else // M * N + return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true); + else + if (interval_lib::detail::is_pos(yu)) // M * P + return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true); + else // M * Z + return I(0, 0, true); + else + if (interval_lib::detail::is_neg(yl)) + if (interval_lib::detail::is_pos(yu)) // N * M + return I(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true); + else // N * N + return I(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true); + else + if (interval_lib::detail::is_pos(yu)) // N * P + return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true); + else // N * Z + return I(0, 0, true); + else + if (interval_lib::detail::is_pos(xu)) + if (interval_lib::detail::is_neg(yl)) + if (interval_lib::detail::is_pos(yu)) // P * M + return I(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true); + else // P * N + return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true); + else + if (interval_lib::detail::is_pos(yu)) // P * P + return I(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true); + else // P * Z + return I(0, 0, true); + else // Z * ? + return I(0, 0, true); +} + +template inline +interval operator*(const T& x, const interval& y) +{ + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + typename Policies::rounding rnd; + const T& yl = y.lower(); + const T& yu = y.upper(); + // x is supposed not to be infinite + if (interval_lib::detail::is_neg(x)) + return I(rnd.mul_down(x, yu), rnd.mul_up(x, yl), true); + else if (interval_lib::detail::is_zero(x)) + return I(0, 0, true); + else + return I(rnd.mul_down(x, yl), rnd.mul_up(x, yu), true); +} + +template inline +interval operator*(const interval& x, const T& y) +{ return y * x; } + +template inline +interval operator/(const interval& x, + const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + if (in_zero(y)) + if (!interval_lib::detail::is_zero(y.lower())) + if (!interval_lib::detail::is_zero(y.upper())) + return interval_lib::detail::div_zero(x); + else + return interval_lib::detail::div_negative(x, y.lower()); + else + if (!interval_lib::detail::is_zero(y.upper())) + return interval_lib::detail::div_positive(x, y.upper()); + else + return interval::empty(); + else + return interval_lib::detail::div_non_zero(x, y); +} + +template inline +interval operator/(const T& x, const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + if (in_zero(y)) + if (!interval_lib::detail::is_zero(y.lower())) + if (!interval_lib::detail::is_zero(y.upper())) + return interval_lib::detail::div_zero(x); + else + return interval_lib::detail::div_negative(x, y.lower()); + else + if (!interval_lib::detail::is_zero(y.upper())) + return interval_lib::detail::div_positive(x, y.upper()); + else + return interval::empty(); + else + return interval_lib::detail::div_non_zero(x, y); +} + +template inline +interval operator/(const interval& x, const T& y) +{ + if (interval_lib::detail::test_input(x, y) || interval_lib::detail::is_zero(y)) + return interval::empty(); + typename Policies::rounding rnd; + const T& xl = x.lower(); + const T& xu = x.upper(); + if (interval_lib::detail::is_neg(y)) + return interval(rnd.div_down(xu, y), rnd.div_up(xl, y), true); + else + return interval(rnd.div_down(xl, y), rnd.div_up(xu, y), true); +} + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP diff --git a/include/boost/numeric/interval/arith2.hpp b/include/boost/numeric/interval/arith2.hpp new file mode 100644 index 0000000..34aae32 --- /dev/null +++ b/include/boost/numeric/interval/arith2.hpp @@ -0,0 +1,198 @@ +#ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP +#define BOOST_NUMERIC_INTERVAL_ARITH2_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace numeric { + +template inline +interval fmod(const interval& x, + const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect >::type I; + const T& yb = interval_lib::detail::is_neg(x.lower()) ? y.lower() : y.upper(); + T n = rnd.int_down(rnd.div_down(x.lower(), yb)); + return (const I&)x - n * (const I&)y; +} + +template inline +interval fmod(const interval& x, const T& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect >::type I; + T n = rnd.int_down(rnd.div_down(x.lower(), y)); + return (const I&)x - n * I(y); +} + +template inline +interval fmod(const T& x, const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect >::type I; + const T& yb = interval_lib::detail::is_neg(x) ? y.lower() : y.upper(); + T n = rnd.int_down(rnd.div_down(x, yb)); + return x - n * (const I&)y; +} + +namespace interval_lib { + +template inline +interval division_part1(const interval& x, + const interval& y, bool& b) +{ + typedef interval I; + b = false; + if (detail::test_input(x, y)) + return I::empty(); + if (in_zero(y)) + if (!detail::is_zero(y.lower())) + if (!detail::is_zero(y.upper())) + return detail::div_zero_part1(x, y, b); + else + return detail::div_negative(x, y.lower()); + else + if (!detail::is_zero(y.upper())) + return detail::div_positive(x, y.upper()); + else + return I::empty(); + else + return detail::div_non_zero(x, y); +} + +template inline +interval division_part2(const interval& x, + const interval& y) +{ + return detail::div_zero_part2(x, y); +} + +template inline +interval multiplicative_inverse(const interval& x) +{ + typedef interval I; + if (detail::test_input(x)) + return I::empty(); + T one = static_cast(1); + typename Policies::rounding rnd; + if (in_zero(x)) { + typedef typename I::checking checking; + if (!detail::is_zero(x.lower())) + if (!detail::is_zero(x.upper())) + return I::whole(); + else + return I(-checking::inf(), rnd.div_up(one, x.lower()), true); + else + if (!detail::is_zero(x.upper())) + return I(rnd.div_down(one, x.upper()), checking::inf(), true); + else + return I::empty(); + } else + return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true); +} + +namespace detail { + +template inline +T pow_aux(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive +{ + T x = x_; + T y = (pwr & 1) ? x_ : 1; + pwr >>= 1; + while (pwr > 0) { + x = rnd.mul_up(x, x); + if (pwr & 1) y = rnd.mul_up(x, y); + pwr >>= 1; + } + return y; +} + +} // namespace detail +} // namespace interval_lib + +template inline +interval pow(const interval& x, int pwr) +{ + using interval_lib::detail::pow_aux; + typedef interval I; + + if (interval_lib::detail::test_input(x)) + return I::empty(); + + if (pwr == 0) + if (interval_lib::detail::is_zero(x.lower()) + && interval_lib::detail::is_zero(x.upper())) + return I::empty(); + else + return I(1); + else if (pwr < 0) + return multiplicative_inverse(pow(x, -pwr)); + + typename Policies::rounding rnd; + + if (interval_lib::detail::is_neg(x.upper())) { // [-2,-1] + T yl = pow_aux(-x.upper(), pwr, rnd); + T yu = pow_aux(-x.lower(), pwr, rnd); + if (pwr & 1) // [-2,-1]^1 + return I(-yu, -yl, true); + else // [-2,-1]^2 + return I(yl, yu, true); + } else if (interval_lib::detail::is_neg(x.lower())) { // [-1,1] + if (pwr & 1) { // [-1,1]^1 + return I(-pow_aux(-x.lower(), pwr, rnd), pow_aux(x.upper(), pwr, rnd), true); + } else { // [-1,1]^2 + BOOST_NUMERIC_INTERVAL_using_max(max); + return I(0, pow_aux(max(-x.lower(), x.upper()), pwr, rnd), true); + } + } else { // [1,2] + return I(pow_aux(x.lower(), pwr, rnd), pow_aux(x.upper(), pwr, rnd), true); + } +} + +template inline +interval sqrt(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x) || interval_lib::detail::is_neg(x.upper())) + return I::empty(); + typename Policies::rounding rnd; + T l = (x.lower() <= static_cast(0)) ? 0 : rnd.sqrt_down(x.lower()); + return I(l, rnd.sqrt_up(x.upper()), true); +} + +template inline +interval square(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + const T& xl = x.lower(); + const T& xu = x.upper(); + if (interval_lib::detail::is_neg(xu)) + return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true); + else if (interval_lib::detail::is_pos(x.lower())) + return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true); + else + return I(0, (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true); +} + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP diff --git a/include/boost/numeric/interval/checking.hpp b/include/boost/numeric/interval/checking.hpp new file mode 100644 index 0000000..aa39968 --- /dev/null +++ b/include/boost/numeric/interval/checking.hpp @@ -0,0 +1,141 @@ +/* Boost interval/checking.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_CHECKING_HPP +#define BOOST_NUMERIC_INTERVAL_CHECKING_HPP + +#include +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { + +struct exception_create_empty +{ + void operator()() + { + throw std::runtime_error("boost::interval: empty interval created"); + } +}; + +struct exception_invalid_number +{ + void operator()() + { + throw std::invalid_argument("boost::interval: invalid number"); + } +}; + +template +struct checking_base +{ + static T inf() + { + assert(std::numeric_limits::has_infinity); + return std::numeric_limits::infinity(); + } + static T nan() + { + assert(std::numeric_limits::has_quiet_NaN); + return std::numeric_limits::quiet_NaN(); + } + static bool is_nan(const T& x) + { + return std::numeric_limits::has_quiet_NaN && (x != x); + } + static T empty_lower() + { + return (std::numeric_limits::has_quiet_NaN ? + std::numeric_limits::quiet_NaN() : T(1)); + } + static T empty_upper() + { + return (std::numeric_limits::has_quiet_NaN ? + std::numeric_limits::quiet_NaN() : T(0)); + } + static bool is_empty(const T& l, const T& u) + { + return !(l <= u); // safety for partial orders + } +}; + +template, + class Exception = exception_create_empty> +struct checking_no_empty: Checking +{ + static T nan() + { + assert(false); + return Checking::nan(); + } + static T empty_lower() + { + Exception()(); + return Checking::empty_lower(); + } + static T empty_upper() + { + Exception()(); + return Checking::empty_upper(); + } + static bool is_empty(const T&, const T&) + { + return false; + } +}; + +template > +struct checking_no_nan: Checking +{ + static bool is_nan(const T&) + { + return false; + } +}; + +template, + class Exception = exception_invalid_number> +struct checking_catch_nan: Checking +{ + static bool is_nan(const T& x) + { + if (Checking::is_nan(x)) Exception()(); + return false; + } +}; + +template +struct checking_strict: + checking_catch_nan > +{}; + +namespace detail { + +template inline bool is_nan(const T& x) { return x != x; } + +} // namespace detail + +} // namespace interval +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_CHECKING_HPP diff --git a/include/boost/numeric/interval/compare.hpp b/include/boost/numeric/interval/compare.hpp new file mode 100644 index 0000000..3832078 --- /dev/null +++ b/include/boost/numeric/interval/compare.hpp @@ -0,0 +1,28 @@ +/* Boost interval/compare.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_COMPARE_HPP +#define BOOST_NUMERIC_INTERVAL_COMPARE_HPP + +#include +#include +#include + +#endif // BOOST_NUMERIC_INTERVAL_COMPARE_HPP diff --git a/include/boost/numeric/interval/compare/explicit.hpp b/include/boost/numeric/interval/compare/explicit.hpp new file mode 100644 index 0000000..87368ce --- /dev/null +++ b/include/boost/numeric/interval/compare/explicit.hpp @@ -0,0 +1,238 @@ +#ifndef BOOST_NUMERIC_INTERVAL_COMPARE_EXPLICIT_HPP +#define BOOST_NUMERIC_INTERVAL_COMPARE_EXPLICIT_HPP + +#include + +namespace boost { +namespace numeric { +namespace interval_lib { + +/* + * Certainly... operations + */ + +template inline +bool cerlt(const interval& x, const interval& y) +{ + return x.upper() < y.lower(); +} + +template inline +bool cerlt(const interval& x, const T& y) +{ + return x.upper() < y; +} + +template inline +bool cerlt(const T& x, const interval& y) +{ + return x < y.lower(); +} + +template inline +bool cerle(const interval& x, const interval& y) +{ + return x.upper() <= y.lower(); +} + +template inline +bool cerle(const interval& x, const T& y) +{ + return x.upper() <= y; +} + +template inline +bool cerle(const T& x, const interval& y) +{ + return x <= y.lower(); +} + +template inline +bool cergt(const interval& x, const interval& y) +{ + return x.lower() > y.upper(); +} + +template inline +bool cergt(const interval& x, const T& y) +{ + return x.lower() > y; +} + +template inline +bool cergt(const T& x, const interval& y) +{ + return x > y.upper(); +} + +template inline +bool cerge(const interval& x, const interval& y) +{ + return x.lower() >= y.upper(); +} + +template inline +bool cerge(const interval& x, const T& y) +{ + return x.lower() >= y; +} + +template inline +bool cerge(const T& x, const interval& y) +{ + return x >= y.upper(); +} + +template inline +bool cereq(const interval& x, const interval& y) +{ + return x.lower() == y.upper() && y.lower() == x.upper(); +} + +template inline +bool cereq(const interval& x, const T& y) +{ + return x.lower() == y && x.upper() == y; +} + +template inline +bool cereq(const T& x, const interval& y) +{ + return x == y.lower() && x == y.upper(); +} + +template inline +bool cerne(const interval& x, const interval& y) +{ + return x.upper() < y.lower() || y.upper() < x.lower(); +} + +template inline +bool cerne(const interval& x, const T& y) +{ + return x.upper() < y || y < x.lower(); +} + +template inline +bool cerne(const T& x, const interval& y) +{ + return x < y.lower() || y.upper() < x; +} + +/* + * Possibly... comparisons + */ + +template inline +bool poslt(const interval& x, const interval& y) +{ + return x.lower() < y.upper(); +} + +template inline +bool poslt(const interval& x, const T& y) +{ + return x.lower() < y; +} + +template inline +bool poslt(const T& x, const interval& y) +{ + return x < y.upper(); +} + +template inline +bool posle(const interval& x, const interval& y) +{ + return x.lower() <= y.upper(); +} + +template inline +bool posle(const interval& x, const T& y) +{ + return x.lower() <= y; +} + +template inline +bool posle(const T& x, const interval& y) +{ + return x <= y.upper(); +} + +template inline +bool posgt(const interval& x, const interval& y) +{ + return x.upper() > y.lower(); +} + +template inline +bool posgt(const interval& x, const T& y) +{ + return x.upper() > y; +} + +template inline +bool posgt(const T& x, const interval & y) +{ + return x > y.lower(); +} + +template inline +bool posge(const interval& x, const interval& y) +{ + return x.upper() >= y.lower(); +} + +template inline +bool posge(const interval& x, const T& y) +{ + return x.upper() >= y; +} + +template inline +bool posge(const T& x, const interval& y) +{ + return x >= y.lower(); +} + +template inline +bool poseq(const interval& x, const interval& y) +{ + return x.upper() >= y.lower() && y.upper() >= x.lower(); +} + +template inline +bool poseq(const interval& x, const T& y) +{ + return x.upper() >= y && y >= x.lower(); +} + +template inline +bool poseq(const T& x, const interval& y) +{ + return x >= y.lower() && y.upper() >= x; +} + +template inline +bool posne(const interval& x, const interval& y) +{ + return x.upper() != y.lower() || y.upper() != x.lower(); +} + +template inline +bool posne(const interval& x, const T& y) +{ + return x.upper() != y || y != x.lower(); +} + +template inline +bool posne(const T& x, const interval& y) +{ + return x != y.lower() || y.upper() != x; +} + +} // namespace interval_lib +} // namespace numeric +} //namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_COMPARE_EXPLICIT_HPP diff --git a/include/boost/numeric/interval/compare/lexicographic.hpp b/include/boost/numeric/interval/compare/lexicographic.hpp new file mode 100644 index 0000000..8b3663e --- /dev/null +++ b/include/boost/numeric/interval/compare/lexicographic.hpp @@ -0,0 +1,113 @@ +#ifndef BOOST_NUMERIC_INTERVAL_COMPARE_LEXICOGRAPHIC_HPP +#define BOOST_NUMERIC_INTERVAL_COMPARE_LEXICOGRAPHIC_HPP + +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace compare { +namespace lexicographic { + +template inline +bool operator<(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + const T& xl = x.lower(); + const T& yl = y.lower(); + return xl < yl || (xl == yl && x.upper() < y.upper()); +} + +template inline +bool operator<(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + return x.lower() < y; +} + +template inline +bool operator<=(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + const T& xl = x.lower(); + const T& yl = y.lower(); + return xl < yl || (xl == yl && x.upper() <= y.upper()); +} + +template inline +bool operator<=(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + const T& xl = x.lower(); + return xl < y || (xl == y && x.upper() <= y); +} + +template inline +bool operator>(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + const T& xl = x.lower(); + const T& yl = y.lower(); + return xl > yl || (xl == yl && x.upper() > y.upper()); +} + +template inline +bool operator>(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + const T& xl = x.lower(); + return xl > y || (xl == y && x.upper() > y); +} + +template inline +bool operator>=(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + const T& xl = x.lower(); + const T& yl = y.lower(); + return xl > yl || (xl == yl && x.upper() >= y.upper()); +} + +template inline +bool operator>=(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + return x.lower() >= y; +} + +template inline +bool operator==(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + return x.lower() == y.lower() && x.upper() == y.upper(); +} + +template inline +bool operator==(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + return x.lower() == y && x.upper() == y; +} + +template inline +bool operator!=(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + return x.lower() != y.lower() || x.upper() != y.upper(); +} + +template inline +bool operator!=(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + return x.lower() != y || x.upper() != y; +} + +} // namespace lexicographic +} // namespace compare +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_COMPARE_LEXICOGRAPHIC_HPP diff --git a/include/boost/numeric/interval/compare/set.hpp b/include/boost/numeric/interval/compare/set.hpp new file mode 100644 index 0000000..b3d848d --- /dev/null +++ b/include/boost/numeric/interval/compare/set.hpp @@ -0,0 +1,55 @@ +#ifndef BOOST_NUMERIC_INTERVAL_COMPARE_SET_HPP +#define BOOST_NUMERIC_INTERVAL_COMPARE_SET_HPP + +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace compare { +namespace set { + +template inline +bool operator<(const interval& x, const interval& y) +{ + return proper_subset(x, y); +} + +template inline +bool operator<=(const interval& x, const interval& y) +{ + return subset(x, y); +} + +template inline +bool operator>(const interval& x, const interval& y) +{ + return proper_subset(y, x); +} + +template inline +bool operator>=(const interval& x, const interval& y) +{ + return subset(y, x); +} + +template inline +bool operator==(const interval& x, const interval& y) +{ + return equal(y, x); +} + +template inline +bool operator!=(const interval& x, const interval& y) +{ + return !equal(y, x); +} + +} // namespace set +} // namespace compare +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_COMPARE_SET_HPP diff --git a/include/boost/numeric/interval/compare/tribool.hpp b/include/boost/numeric/interval/compare/tribool.hpp new file mode 100644 index 0000000..c283caf --- /dev/null +++ b/include/boost/numeric/interval/compare/tribool.hpp @@ -0,0 +1,129 @@ +#ifndef BOOST_NUMERIC_INTERVAL_COMPARE_TRIBOOL_HPP +#define BOOST_NUMERIC_INTERVAL_COMPARE_TRIBOOL_HPP + +#include +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace compare { +namespace tribool { + +template inline +tribool operator<(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() < y.lower()) return true; + if (x.lower() >= y.upper()) return false; + return indeterminate; +} + +template inline +tribool operator<(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() < y) return true; + if (x.lower() >= y) return false; + return indeterminate; +} + +template inline +tribool operator<=(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() <= y.lower()) return true; + if (x.lower() > y.upper()) return false; + return indeterminate; +} + +template inline +tribool operator<=(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() <= y) return true; + if (x.lower() > y) return false; + return indeterminate; +} + +template inline +tribool operator>(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.lower() > y.upper()) return true; + if (x.upper() <= y.lower()) return false; + return indeterminate; +} + +template inline +tribool operator>(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.lower() > y) return true; + if (x.upper() <= y) return false; + return indeterminate; +} + +template inline +tribool operator>=(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.lower() >= y.upper()) return true; + if (x.upper() < y.lower()) return false; + return indeterminate; +} + +template inline +tribool operator>=(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.lower() >= y) return true; + if (x.upper() < y) return false; + return indeterminate; +} + +template inline +tribool operator==(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() == y.lower() && x.lower() == y.upper()) return true; + if (x.upper() < y.lower() || x.lower() > y.upper()) return false; + return indeterminate; +} + +template inline +tribool operator==(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() == y && x.lower() == y) return true; + if (x.upper() < y || x.lower() > y) return false; + return indeterminate; +} + +template inline +tribool operator!=(const interval& x, const interval& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() < y.lower() || x.lower() > y.upper()) return true; + if (x.upper() == y.lower() && x.lower() == y.upper()) return false; + return indeterminate; +} + +template inline +tribool operator!=(const interval& x, const T& y) +{ + if (detail::test_input(x, y)) throw comparison_error(); + if (x.upper() < y || x.lower() > y) return true; + if (x.upper() == y && x.lower() == y) return false; + return indeterminate; +} + +} // namespace tribool +} // namespace compare +} // namespace interval_lib +} // namespace numeric +} // namespace boost + + +#endif // BOOST_NUMERIC_INTERVAL_COMPARE_TRIBOOL_HPP diff --git a/include/boost/numeric/interval/constants.hpp b/include/boost/numeric/interval/constants.hpp new file mode 100644 index 0000000..0ef7bb9 --- /dev/null +++ b/include/boost/numeric/interval/constants.hpp @@ -0,0 +1,96 @@ +/* Boost interval/constants.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_CONSTANTS_HPP +#define BOOST_NUMERIC_INTERVAL_CONSTANTS_HPP + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace constants { + +// These constants should be exactly computed. +// Decimal representations wouldn't do it since the standard doesn't +// specify the rounding (even nearest) that should be used. + +static const float pi_f_l = 13176794.0f/(1<<22); +static const float pi_f_u = 13176795.0f/(1<<22); +static const double pi_d_l = (3373259426.0 + 273688.0 / (1<<21)) / (1<<30); +static const double pi_d_u = (3373259426.0 + 273689.0 / (1<<21)) / (1<<30); + +template inline T pi_lower() { return 3; } +template inline T pi_upper() { return 4; } +template inline T pi_half_lower() { return 1; } +template inline T pi_half_upper() { return 2; } +template inline T pi_twice_lower() { return 6; } +template inline T pi_twice_upper() { return 7; } + +template<> inline float pi_lower() { return pi_f_l; } +template<> inline float pi_upper() { return pi_f_u; } +template<> inline float pi_half_lower() { return pi_f_l / 2; } +template<> inline float pi_half_upper() { return pi_f_u / 2; } +template<> inline float pi_twice_lower() { return pi_f_l * 2; } +template<> inline float pi_twice_upper() { return pi_f_u * 2; } + +template<> inline double pi_lower() { return pi_d_l; } +template<> inline double pi_upper() { return pi_d_u; } +template<> inline double pi_half_lower() { return pi_d_l / 2; } +template<> inline double pi_half_upper() { return pi_d_u / 2; } +template<> inline double pi_twice_lower() { return pi_d_l * 2; } +template<> inline double pi_twice_upper() { return pi_d_u * 2; } + +template<> inline long double pi_lower() { return pi_d_l; } +template<> inline long double pi_upper() { return pi_d_u; } +template<> inline long double pi_half_lower() { return pi_d_l / 2; } +template<> inline long double pi_half_upper() { return pi_d_u / 2; } +template<> inline long double pi_twice_lower() { return pi_d_l * 2; } +template<> inline long double pi_twice_upper() { return pi_d_u * 2; } + +} // namespace constants + +template inline +I pi() +{ + typedef typename I::base_type T; + return I(constants::pi_lower(), + constants::pi_upper(), true); +} + +template inline +I pi_half() +{ + typedef typename I::base_type T; + return I(constants::pi_half_lower(), + constants::pi_half_upper(), true); +} + +template inline +I pi_twice() +{ + typedef typename I::base_type T; + return I(constants::pi_twice_lower(), + constants::pi_twice_upper(), true); +} + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_CONSTANTS_HPP diff --git a/include/boost/numeric/interval/detail/bcc_rounding_control.hpp b/include/boost/numeric/interval/detail/bcc_rounding_control.hpp new file mode 100644 index 0000000..e59e4c1 --- /dev/null +++ b/include/boost/numeric/interval/detail/bcc_rounding_control.hpp @@ -0,0 +1,67 @@ +/* Boost interval/detail/bcc_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_BCC_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_BCC_ROUNDING_CONTROL_HPP + +#ifndef __BORLANDC__ +# error This header is only intended for Borland C++. +#endif + +#ifndef _M_IX86 +# error This header only works on x86 CPUs. +#endif + +#include // Borland C++ rounding control + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +#ifndef BOOST_NUMERIC_INTERVAL_KEEP_EXCEPTIONS_FOR_BCC +extern "C" { unsigned int _RTLENTRY _fm_init(void); } + +struct borland_workaround { + borland_workaround() { _fm_init(); } +}; + +static borland_workaround borland_workaround_exec; +#endif // BOOST_NUMERIC_INTERVAL_KEEP_EXCEPTIONS_FOR_BCC + +__inline double rint(double) +{ __emit__(0xD9); __emit__(0xFC); /* asm FRNDINT */ } + +struct x86_rounding +{ + typedef unsigned int rounding_mode; + static void get_rounding_mode(rounding_mode& mode) + { mode = _control87(0, 0); } + static void set_rounding_mode(const rounding_mode mode) + { _control87(mode, 0xffff); } + static double to_int(const double& x) { return rint(x); } +}; + +} // namespace detail +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif /* BOOST_NUMERIC_INTERVAL_DETAIL_BCC_ROUNDING_CONTROL_HPP */ diff --git a/include/boost/numeric/interval/detail/bugs.hpp b/include/boost/numeric/interval/detail/bugs.hpp new file mode 100644 index 0000000..4132b0e --- /dev/null +++ b/include/boost/numeric/interval/detail/bugs.hpp @@ -0,0 +1,83 @@ +/* Boost interval/detail/bugs.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_BUGS +#define BOOST_NUMERIC_INTERVAL_DETAIL_BUGS + +#include + +#if defined(__GLIBC__) && !defined(__GLIBCPP__) && (defined(__USE_MISC) || defined(__USE_XOPEN_EXTENDED) || defined(__USE_ISOC99)) && !defined(__ICC) +# define BOOST_HAVE_INV_HYPERBOLIC +#endif + +#ifndef BOOST_HAVE_INV_HYPERBOLIC +# define BOOST_NUMERIC_INTERVAL_using_ahyp(a) +#endif + +#if defined(BOOST_NO_STDC_NAMESPACE) +# define BOOST_NUMERIC_INTERVAL_using_max(a) ::a +# define BOOST_NUMERIC_INTERVAL_using_math(a) ::a +# ifndef BOOST_NUMERIC_INTERVAL_using_ahyp +# define BOOST_NUMERIC_INTERVAL_using_ahyp(a) ::a +# endif +#else +# define BOOST_NUMERIC_INTERVAL_using_max(a) using std::a +# define BOOST_NUMERIC_INTERVAL_using_math(a) using std::a +# ifndef BOOST_NUMERIC_INTERVAL_using_ahyp +# define BOOST_NUMERIC_INTERVAL_using_ahyp(a) using std::a +# endif +#endif + +#if __GNUC__ <= 2 +// cf PR c++/1981 for a description of the bug +#include +#include +namespace boost { +namespace numeric { + using std::min; + using std::max; + using std::sqrt; + using std::exp; + using std::log; + using std::cos; + using std::tan; + using std::asin; + using std::acos; + using std::atan; + using std::ceil; + using std::floor; + using std::sinh; + using std::cosh; + using std::tanh; +# undef BOOST_NUMERIC_INTERVAL_using_max +# undef BOOST_NUMERIC_INTERVAL_using_math +# define BOOST_NUMERIC_INTERVAL_using_max(a) +# define BOOST_NUMERIC_INTERVAL_using_math(a) +# if defined(BOOST_HAVE_INV_HYPERBOLIC) + using std::asinh; + using std::acosh; + using std::atanh; +# undef BOOST_NUMERIC_INTERVAL_using_ahyp +# define BOOST_NUMERIC_INTERVAL_using_ahyp(a) +# endif +} // namespace numeric +} // namespace boost +#endif + +#endif // BOOST_NUMERIC_INTERVAL_DETAIL_BUGS diff --git a/include/boost/numeric/interval/detail/c99_rounding_control.hpp b/include/boost/numeric/interval/detail/c99_rounding_control.hpp new file mode 100644 index 0000000..05e5722 --- /dev/null +++ b/include/boost/numeric/interval/detail/c99_rounding_control.hpp @@ -0,0 +1,54 @@ +/* Boost interval/detail/c99_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_C99_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_C99_ROUNDING_CONTROL_HPP + +#include + +namespace boost { +namespace nuemric { +namespace interval_lib { +namespace detail { + +struct c99_rounding_control +{ + template + static T force_rounding(const T& r) { volatile T r_ = r; return r_; } +}; + +} // namespace detail + +template<> +struct rounding_control: + detail::c99_rounding_control { }; + +template<> +struct rounding_control: + detail::c99_rounding_control { }; + +template<> +struct rounding_control: + detail::c99_rounding_control { }; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_DETAIL_C99_ROUNDING_CONTROL_HPP diff --git a/include/boost/numeric/interval/detail/c99sub_rounding_control.hpp b/include/boost/numeric/interval/detail/c99sub_rounding_control.hpp new file mode 100644 index 0000000..d4c7dcf --- /dev/null +++ b/include/boost/numeric/interval/detail/c99sub_rounding_control.hpp @@ -0,0 +1,50 @@ +/* Boost interval/detail/c99sub_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_C99SUB_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_C99SUB_ROUNDING_CONTROL_HPP + +#include // ISO C 99 rounding mode control + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +extern "C" { double rint(double); } + +struct c99_rounding +{ + typedef int rounding_mode; + + static void set_rounding_mode(const rounding_mode mode) { fesetround(mode); } + static void get_rounding_mode(rounding_mode &mode) { mode = fegetround(); } + static void downward() { set_rounding_mode(FE_DOWNWARD); } + static void upward() { set_rounding_mode(FE_UPWARD); } + static void to_nearest() { set_rounding_mode(FE_TONEAREST); } + static void toward_zero() { set_rounding_mode(FE_TOWARDZERO); } + + template + static T to_int(const T& r) { return rint(r); } +}; + +} // namespace detail +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_DETAIL_C99SUB_ROUBDING_CONTROL_HPP diff --git a/include/boost/numeric/interval/detail/division.hpp b/include/boost/numeric/interval/detail/division.hpp new file mode 100644 index 0000000..15f9b0a --- /dev/null +++ b/include/boost/numeric/interval/detail/division.hpp @@ -0,0 +1,183 @@ +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_DIVISION_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_DIVISION_HPP + +#include +#include +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +template inline +interval div_non_zero(const interval& x, + const interval& y) +{ + // assert(!in_zero(y)); + typename Policies::rounding rnd; + typedef interval I; + const T& xl = x.lower(); + const T& xu = x.upper(); + const T& yl = y.lower(); + const T& yu = y.upper(); + if (is_neg(xu)) + if (is_neg(yu)) + return I(rnd.div_down(xu, yl), rnd.div_up(xl, yu), true); + else + return I(rnd.div_down(xl, yl), rnd.div_up(xu, yu), true); + else if (is_neg(xl)) + if (is_neg(yu)) + return I(rnd.div_down(xu, yu), rnd.div_up(xl, yu), true); + else + return I(rnd.div_down(xl, yl), rnd.div_up(xu, yl), true); + else + if (is_neg(yu)) + return I(rnd.div_down(xu, yu), rnd.div_up(xl, yl), true); + else + return I(rnd.div_down(xl, yu), rnd.div_up(xu, yl), true); +} + +template inline +interval div_non_zero(const T& x, const interval& y) +{ + // assert(!in_zero(y)); + typename Policies::rounding rnd; + typedef interval I; + const T& yl = y.lower(); + const T& yu = y.upper(); + if (is_neg(x)) + return I(rnd.div_down(x, yl), rnd.div_up(x, yu), true); + else + return I(rnd.div_down(x, yu), rnd.div_up(x, yl), true); +} + +template inline +interval div_positive(const interval& x, const T& yu) +{ + // assert(yu > T(0)); + if (is_zero(x)) return x; + typename Policies::rounding rnd; + typedef interval I; + const T& xl = x.lower(); + const T& xu = x.upper(); + typedef typename I::checking checking; + const T& inf = checking::inf(); + if (is_neg(xu)) + return I(-inf, rnd.div_up(xu, yu), true); + else if (is_neg(xl)) + return I(-inf, inf, true); + else + return I(rnd.div_down(xl, yu), inf, true); +} + +template inline +interval div_positive(const T& x, const T& yu) +{ + // assert(yu > T(0)); + typedef interval I; + if (is_zero(x)) return I(0, 0, true); + typename Policies::rounding rnd; + typedef typename I::checking checking; + const T& inf = checking::inf(); + if (is_neg(x)) + return I(-inf, rnd.div_up(x, yu), true); + else + return I(rnd.div_down(x, yu), inf, true); +} + +template inline +interval div_negative(const interval& x, const T& yl) +{ + // assert(yl < T(0)); + if (is_zero(x.lower()) && is_zero(x.upper())) + return x; + typename Policies::rounding rnd; + typedef interval I; + const T& xl = x.lower(); + const T& xu = x.upper(); + typedef typename I::checking checking; + const T& inf = checking::inf(); + if (is_neg(xu)) + return I(rnd.div_down(xu, yl), inf, true); + else if (is_neg(xl)) + return I(-inf, inf, true); + else + return I(-inf, rnd.div_up(xl, yl), true); +} + +template inline +interval div_negative(const T& x, const T& yl) +{ + // assert(yl < T(0)); + typedef interval I; + if (is_zero(x)) return I(0, 0, true); + typename Policies::rounding rnd; + typedef typename I::checking checking; + const T& inf = checking::inf(); + if (is_neg(x)) + return I(rnd.div_down(x, yl), inf, true); + else + return I(-inf, rnd.div_up(x, yl), true); +} + +template inline +interval div_zero(const interval& x) +{ + if (is_zero(x.lower()) && is_zero(x.upper())) + return x; + else return interval::whole(); +} + +template inline +interval div_zero(const T& x) +{ + if (is_zero(x)) return interval(0, 0, true); + else return interval::whole(); +} + +template inline +interval div_zero_part1(const interval& x, + const interval& y, bool& b) +{ + // assert(y.lower() < 0 && y.upper() > 0); + if (is_zero(x.lower()) && is_zero(x.upper())) + { b = false; return x; } + typename Policies::rounding rnd; + typedef interval I; + const T& xl = x.lower(); + const T& xu = x.upper(); + const T& yl = y.lower(); + const T& yu = y.upper(); + typedef typename I::checking checking; + const T& inf = checking::inf(); + if (is_neg(xu)) + { b = true; return I(-inf, rnd.div_up(xu, yu), true); } + else if (is_neg(xl)) + { b = false; return I(-inf, inf, true); } + else + { b = true; return I(-inf, rnd.div_up(xl, yl), true); } +} + +template inline +interval div_zero_part2(const interval& x, + const interval& y) +{ + // assert(y.lower() < 0 && y.upper() > 0 && (div_zero_part1(x, y, b), b)); + typename Policies::rounding rnd; + typedef interval I; + typedef typename I::checking checking; + const T& inf = checking::inf(); + if (is_neg(x.upper())) + return I(rnd.div_down(x.upper(), y.lower()), inf, true); + else + return I(rnd.div_down(x.lower(), y.upper()), inf, true); +} + +} // namespace detail +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_DETAIL_DIVISION_HPP diff --git a/include/boost/numeric/interval/detail/interval_prototype.hpp b/include/boost/numeric/interval/detail/interval_prototype.hpp new file mode 100644 index 0000000..e84241f --- /dev/null +++ b/include/boost/numeric/interval/detail/interval_prototype.hpp @@ -0,0 +1,32 @@ +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_INTERVAL_PROTOTYPE_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_INTERVAL_PROTOTYPE_HPP + +namespace boost { +namespace numeric { + +namespace interval_lib { + +template struct rounded_math; +template struct checking_strict; +class comparison_error; +template struct policies; + +/* + * default policies class + */ + +template +struct default_policies +{ + typedef policies, checking_strict > type; +}; + +} // namespace interval_lib + +template::type > +class interval; + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_DETAIL_INTERVAL_PROTOTYPE_HPP diff --git a/include/boost/numeric/interval/detail/msvc_rounding_control.hpp b/include/boost/numeric/interval/detail/msvc_rounding_control.hpp new file mode 100644 index 0000000..eedd08e --- /dev/null +++ b/include/boost/numeric/interval/detail/msvc_rounding_control.hpp @@ -0,0 +1,52 @@ +/* Boost interval/detail/msvc_rounding_control.hpp file + * + * Copyright Maarten Keijzer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_MSVC_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_MSVC_ROUNDING_CONTROL_HPP + +#ifndef _MSC_VER +# error This header is only intended for MSVC, but might work for Borland as well +#endif + +#include // MSVC rounding control + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +extern "C" { double rint(double); } + +struct x86_rounding +{ + typedef unsigned int rounding_mode; + static void get_rounding_mode(rounding_mode& mode) + { mode = _control87(0, 0); } + static void set_rounding_mode(const rounding_mode mode) + { _control87(mode, 0xffff); } + static double to_int(const double& x) { return rint(x); } +}; + +} // namespace detail +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif /* BOOST_NUMERIC_INTERVAL_DETAIL_MSVC_ROUNDING_CONTROL_HPP */ diff --git a/include/boost/numeric/interval/detail/ppc_rounding_control.hpp b/include/boost/numeric/interval/detail/ppc_rounding_control.hpp new file mode 100644 index 0000000..bb4d1b6 --- /dev/null +++ b/include/boost/numeric/interval/detail/ppc_rounding_control.hpp @@ -0,0 +1,103 @@ +/* Boost interval/detail/ppc_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_PPC_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_PPC_ROUNDING_CONTROL_HPP + +#ifndef __GNUC__ +#error This header only works with GNU CC. +#endif + +#if !defined(powerpc) && !defined(__powerpc__) && !defined(__ppc__) +#error This header only works on PPC CPUs. +#endif + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +typedef union { + long long int imode; + double dmode; +} rounding_mode_struct; + +static const rounding_mode_struct mode_upward = { 0xFFF8000000000002LL }; +static const rounding_mode_struct mode_downward = { 0xFFF8000000000003LL }; +static const rounding_mode_struct mode_to_nearest = { 0xFFF8000000000001LL }; +static const rounding_mode_struct mode_toward_zero = { 0xFFF8000000000000LL }; + +struct ppc_rounding_control +{ + typedef double rounding_mode; + + static void set_rounding_mode(const rounding_mode mode) + { __asm__ __volatile__ ("mtfsf 255,%0" : : "f"(mode)); } + + static void get_rounding_mode(rounding_mode& mode) + { __asm__ __volatile__ ("mffs %0" : "=f"(mode)); } + + static void downward() { set_rounding_mode(mode_downward.dmode); } + static void upward() { set_rounding_mode(mode_upward.dmode); } + static void to_nearest() { set_rounding_mode(mode_to_nearest.dmode); } + static void toward_zero() { set_rounding_mode(mode_toward_zero.dmode); } +}; + +} // namespace detail + +extern "C" { + float rintf(float); + double rint(double); +} + +template<> +struct rounding_control: + detail::ppc_rounding_control +{ + static float force_rounding(const float r) + { + float tmp; + __asm__ __volatile__ ("frsp %0, %1" : "=f" (tmp) : "f" (r)); + return tmp; + } + static float to_int(const float& x) { return rintf(x); } +}; + +template<> +struct rounding_control: + detail::ppc_rounding_control +{ + static const double & force_rounding(const double& r) { return r; } + static double to_int(const double& r) { return rint(r); } +}; + +template<> +struct rounding_control: + detail::ppc_rounding_control +{ + static const long double & force_rounding(const long double& r) { return r; } + static long double to_int(const long double& r) { return rint(r); } +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif /* BOOST_NUMERIC_INTERVAL_DETAIL_PPC_ROUNDING_CONTROL_HPP */ diff --git a/include/boost/numeric/interval/detail/sparc_rounding_control.hpp b/include/boost/numeric/interval/detail/sparc_rounding_control.hpp new file mode 100644 index 0000000..f8619bb --- /dev/null +++ b/include/boost/numeric/interval/detail/sparc_rounding_control.hpp @@ -0,0 +1,120 @@ +/* Boost interval/detail/sparc_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + + /* The basic code in this file was kindly provided by Jeremy Siek. */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_SPARC_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_SPARC_ROUNDING_CONTROL_HPP + +#if !defined(sparc) && !defined(__sparc__) +# error This header is only intended for SPARC CPUs. +#endif + +#ifdef __SUNPRO_CC +# include +#endif + + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +struct sparc_rounding_control +{ + typedef unsigned int rounding_mode; + + static void set_rounding_mode(const rounding_mode& mode) + { +# if defined(__GNUC__) + __asm__ __volatile__("ld %0, %%fsr" : : "m"(mode)); +# elif defined (__SUNPRO_CC) + fpsetround(fp_rnd(mode)); +# elif defined(__KCC) + asm("sethi %hi(mode), %o1"); + asm("ld [%o1+%lo(mode)], %fsr"); +# else +# error Unsupported compiler for Sparc rounding control. +# endif + } + + static void get_rounding_mode(rounding_mode& mode) + { +# if defined(__GNUC__) + __asm__ __volatile__("st %%fsr, %0" : "=m"(mode)); +# elif defined (__SUNPRO_CC) + mode = fpgetround(); +# elif defined(__KCC) +# error KCC on Sun SPARC get_round_mode: please fix me + asm("st %fsr, [mode]"); +# else +# error Unsupported compiler for Sparc rounding control. +# endif + } + +#if defined(__SUNPRO_CC) + static void downward() { set_rounding_mode(FP_RM); } + static void upward() { set_rounding_mode(FP_RP); } + static void to_nearest() { set_rounding_mode(FP_RN); } + static void toward_zero() { set_rounding_mode(FP_RZ); } +#else + static void downward() { set_rounding_mode(0xc0000000); } + static void upward() { set_rounding_mode(0x80000000); } + static void to_nearest() { set_rounding_mode(0x00000000); } + static void toward_zero() { set_rounding_mode(0x40000000); } +#endif +}; + +} // namespace detail + +extern "C" { + float rintf(float); + double rint(double); +} + +template<> +struct rounding_control: + detail::sparc_rounding_control +{ + static const float& force_rounding(const float& x) { return x; } + static float to_int(const float& x) { return rintf(x); } +}; + +template<> +struct rounding_control: + detail::sparc_rounding_control +{ + static const double& force_rounding(const double& x) { return x; } + static double to_int(const double& x) { return rint(x); } +}; + +template<> +struct rounding_control: + detail::sparc_rounding_control +{ + static const long double& force_rounding(const long double& x) { return x; } + static long double to_int(const long double& x) { return rint(x); } +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif /* BOOST_NUMERIC_INTERVAL_DETAIL_SPARC_ROUNDING_CONTROL_HPP */ diff --git a/include/boost/numeric/interval/detail/test_input.hpp b/include/boost/numeric/interval/detail/test_input.hpp new file mode 100644 index 0000000..fc1c4e5 --- /dev/null +++ b/include/boost/numeric/interval/detail/test_input.hpp @@ -0,0 +1,63 @@ +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_TEST_INPUT_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_TEST_INPUT_HPP + +#include + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +template inline +bool is_neg(const T& x) { return x < static_cast(0); } + +template inline +bool is_zero(const T& x) { return x == static_cast(0); } + +template inline +bool is_pos(const T& x) { return x > static_cast(0); } + +template inline +bool test_input(const interval& x) { + typedef typename Policies::checking checking; + return checking::is_empty(x.lower(), x.upper()); +} + +template inline +bool test_input(const interval& x, const interval& y) { + typedef typename Policies1::checking checking1; + typedef typename Policies2::checking checking2; + return checking1::is_empty(x.lower(), x.upper()) || + checking2::is_empty(y.lower(), y.upper()); +} + +template inline +bool test_input(const T& x, const interval& y) { + typedef typename Policies::checking checking; + return checking::is_nan(x) || checking::is_empty(y.lower(), y.upper()); +} + +template inline +bool test_input(const interval& x, const T& y) { + typedef typename Policies::checking checking; + return checking::is_empty(x.lower(), x.upper()) || checking::is_nan(y); +} + +template inline +bool test_input(const T& x) { + typedef typename Policies::checking checking; + return checking::is_nan(x); +} + +template inline +bool test_input(const T& x, const T& y) { + typedef typename Policies::checking checking; + return checking::is_nan(x) || checking::is_nan(y); +} + +} // namespace detail +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_DETAIL_TEST_INPUT_HPP diff --git a/include/boost/numeric/interval/detail/x86_rounding_control.hpp b/include/boost/numeric/interval/detail/x86_rounding_control.hpp new file mode 100644 index 0000000..5179fe0 --- /dev/null +++ b/include/boost/numeric/interval/detail/x86_rounding_control.hpp @@ -0,0 +1,93 @@ +/* Boost interval/detail/x86_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_X86_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_X86_ROUNDING_CONTROL_HPP + +#ifdef __GNUC__ +# include +#elif defined(__BORLANDC__) +# include +#elif defined(_MSC_VER) +# include +#elif defined(__MWERKS__) || defined(__ICC) +# define BOOST_NUMERIC_INTERVAL_USE_C99_SUBSYSTEM +# include +#else +# error Unsupported C++ compiler. +#endif + +namespace boost { +namespace numeric { +namespace interval_lib { + +namespace detail { + +#ifdef BOOST_NUMERIC_INTERVAL_USE_C99_SUBSYSTEM +typedef c99_rounding x86_rounding_control; +#undef BOOST_NUMERIC_INTERVAL_USE_C99_SUBSYSTEM +#else +struct fpu_rounding_modes +{ + unsigned short to_nearest; + unsigned short downward; + unsigned short upward; + unsigned short toward_zero; +}; + +// exceptions masked, extended precision +// hardware default is 0x037f (0x1000 only has a meaning on 287) +static const fpu_rounding_modes rnd_mode = { 0x137f, 0x177f, 0x1b7f, 0x1f7f }; + +struct x86_rounding_control: x86_rounding +{ + static void to_nearest() { set_rounding_mode(rnd_mode.to_nearest); } + static void downward() { set_rounding_mode(rnd_mode.downward); } + static void upward() { set_rounding_mode(rnd_mode.upward); } + static void toward_zero() { set_rounding_mode(rnd_mode.toward_zero); } +}; +#endif // BOOST_NUMERIC_INTERVAL_USE_C99_SUBSYSTEM + +} // namespace detail + +template<> +struct rounding_control: detail::x86_rounding_control +{ + static float force_rounding(const float& r) + { volatile float r_ = r; return r_; } +}; + +template<> +struct rounding_control: detail::x86_rounding_control +{ + static double force_rounding(const double& r) + { volatile double r_ = r; return r_; } +}; + +template<> +struct rounding_control: detail::x86_rounding_control +{ + static const long double& force_rounding(const long double& r) { return r; } +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif /* BOOST_NUMERIC_INTERVAL_DETAIL_X86_ROUNDING_CONTROL_HPP */ diff --git a/include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp b/include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp new file mode 100644 index 0000000..102d8d8 --- /dev/null +++ b/include/boost/numeric/interval/detail/x86gcc_rounding_control.hpp @@ -0,0 +1,61 @@ +/* Boost interval/detail/x86gcc_rounding_control.hpp file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_DETAIL_X86GCC_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_DETAIL_X86GCC_ROUNDING_CONTROL_HPP + +#ifndef __GNUC__ +# error This header only works with GNU CC. +#endif + +#ifndef __i386__ +# error This header only works on x86 CPUs. +#endif + +namespace boost { +namespace numeric { +namespace interval_lib { +namespace detail { + +struct x86_rounding +{ + typedef unsigned short rounding_mode; + + static void set_rounding_mode(const rounding_mode& mode) + { __asm__ __volatile__ ("fldcw %0" : : "m"(mode)); } + + static void get_rounding_mode(rounding_mode& mode) + { __asm__ __volatile__ ("fnstcw %0" : "=m"(mode)); } + + template + static T to_int(T r) + { + T r_; + __asm__ ("frndint" : "=&t"(r_) : "0"(r)); + return r_; + } +}; + +} // namespace detail +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif /* BOOST_NUMERIC_INTERVAL_DETAIL_X86GCC_ROUNDING_CONTROL_HPP */ diff --git a/include/boost/numeric/interval/ext/x86_fast_rounding_control.hpp b/include/boost/numeric/interval/ext/x86_fast_rounding_control.hpp new file mode 100644 index 0000000..18b97dd --- /dev/null +++ b/include/boost/numeric/interval/ext/x86_fast_rounding_control.hpp @@ -0,0 +1,54 @@ +#ifndef BOOST_NUMERIC_INTERVAL_EXT_X86_FAST_ROUNDING_CONTROL_HPP +#define BOOST_NUMERIC_INTERVAL_EXT_X86_FAST_ROUNDING_CONTROL_HPP + +namespace boost { +namespace numeric { +namespace interval_lib { + +namespace detail { + +// exceptions masked, expected precision (the mask is 0x0300) +static const fpu_rounding_modes rnd_mode_f = { 0x107f, 0x147f, 0x187f, 0x1c7f }; +static const fpu_rounding_modes rnd_mode_d = { 0x127f, 0x167f, 0x1a7f, 0x1e7f }; +static const fpu_rounding_modes rnd_mode_l = { 0x137f, 0x177f, 0x1b7f, 0x1f7f }; + +} // namespace detail + +template +struct x86_fast_rounding_control; + +template<> +struct x86_fast_rounding_control: detail::x86_rounding +{ + static void to_nearest() { set_rounding_mode(detail::rnd_mode_f.to_nearest); } + static void downward() { set_rounding_mode(detail::rnd_mode_f.downward); } + static void upward() { set_rounding_mode(detail::rnd_mode_f.upward); } + static void toward_zero() { set_rounding_mode(detail::rnd_mode_f.toward_zero); } + static const float& force_rounding(const float& r) { return r; } +}; + +template<> +struct x86_fast_rounding_control: detail::x86_rounding +{ + static void to_nearest() { set_rounding_mode(detail::rnd_mode_d.to_nearest); } + static void downward() { set_rounding_mode(detail::rnd_mode_d.downward); } + static void upward() { set_rounding_mode(detail::rnd_mode_d.upward); } + static void toward_zero() { set_rounding_mode(detail::rnd_mode_d.toward_zero); } + static const double& force_rounding(const double& r) { return r; } +}; + +template<> +struct x86_fast_rounding_control: detail::x86_rounding +{ + static void to_nearest() { set_rounding_mode(detail::rnd_mode_l.to_nearest); } + static void downward() { set_rounding_mode(detail::rnd_mode_l.downward); } + static void upward() { set_rounding_mode(detail::rnd_mode_l.upward); } + static void toward_zero() { set_rounding_mode(detail::rnd_mode_l.toward_zero); } + static const long double& force_rounding(const long double& r) { return r; } +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_EXT_X86_FAST_ROUNDING_CONTROL_HPP diff --git a/include/boost/numeric/interval/hw_rounding.hpp b/include/boost/numeric/interval/hw_rounding.hpp new file mode 100644 index 0000000..0cb23f1 --- /dev/null +++ b/include/boost/numeric/interval/hw_rounding.hpp @@ -0,0 +1,47 @@ +#ifndef BOOST_NUMERIC_INTERVAL_HW_ROUNDING_HPP +#define BOOST_NUMERIC_INTERVAL_HW_ROUNDING_HPP + +#include +#include + +// define appropriate specialization of rounding_control for built-in types +#if defined(__i386__) || defined(__BORLANDC__) || defined(BOOST_MSVC) +# include +#elif defined(powerpc) || defined(__powerpc__) || defined(__ppc__) +# include +#elif defined(sparc) || defined(__sparc__) +# include +#elif defined(__USE_ISOC99) +# include +#else +# error Boost::interval: Please specify rounding control mechanism. +#endif + +namespace boost { +namespace numeric { +namespace interval_lib { + + /* + * Three specializations of rounded_math + */ + +template<> +struct rounded_math + : save_state > +{}; + +template<> +struct rounded_math + : save_state > +{}; + +template<> +struct rounded_math + : save_state > +{}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_HW_ROUNDING_HPP diff --git a/include/boost/numeric/interval/interval.hpp b/include/boost/numeric/interval/interval.hpp new file mode 100644 index 0000000..fff6d97 --- /dev/null +++ b/include/boost/numeric/interval/interval.hpp @@ -0,0 +1,345 @@ +#ifndef BOOST_NUMERIC_INTERVAL_INTERVAL_HPP +#define BOOST_NUMERIC_INTERVAL_INTERVAL_HPP + +#include +#include + +namespace boost { +namespace numeric { + +namespace interval_lib { + +class comparison_error + : public std::runtime_error +{ +public: + comparison_error() + : std::runtime_error("boost::interval: uncertain comparison") + { } +}; + +} // namespace interval_lib + +/* + * interval class + */ + +template +class interval +{ +private: + struct interval_holder; + struct number_holder; +public: + typedef T base_type; + typedef Policies traits_type; + + interval(const T& v = static_cast(0)); + interval(const T& l, const T& u); + + // The following needs to be defined in the class body for VC++. + template + interval(const interval& r) + : low(r.lower()), up(r.upper()) + { + typedef typename Policies2::checking checking2; + if (checking2::is_empty(low, up)) set_empty(); + } + + // compiler-generated copy constructor and assignment operator are fine + + interval& operator=(const T& x); + void assign(const T& l, const T& u); + + static interval empty(); + static interval whole(); + static interval hull(const T& x, const T& y); + + const T& lower() const; + const T& upper() const; + + interval& operator+= (const T& r); + interval& operator+= (const interval& r); + interval& operator-= (const T& r); + interval& operator-= (const interval& r); + interval& operator*= (const T& r); + interval& operator*= (const interval& r); + interval& operator/= (const T& r); + interval& operator/= (const interval& r); + + bool operator< (const interval_holder& r) const; + bool operator> (const interval_holder& r) const; + bool operator<= (const interval_holder& r) const; + bool operator>= (const interval_holder& r) const; + bool operator== (const interval_holder& r) const; + bool operator!= (const interval_holder& r) const; + + bool operator< (const number_holder& r) const; + bool operator> (const number_holder& r) const; + bool operator<= (const number_holder& r) const; + bool operator>= (const number_holder& r) const; + bool operator== (const number_holder& r) const; + bool operator!= (const number_holder& r) const; + + // the following is for internal use only, it is not a published interface + // nevertheless, it's public because friends don't always work correctly. + interval(const T& l, const T& u, bool): low(l), up(u) {} + void set_empty(); + void set_whole(); + void set(const T& l, const T& u); + +private: + struct interval_holder { + template + interval_holder(const interval& r) + : low(r.lower()), up(r.upper()) + { + typedef typename Policies2::checking checking2; + if (checking2::is_empty(low, up)) + throw interval_lib::comparison_error(); + } + + const T& low; + const T& up; + }; + + struct number_holder { + number_holder(const T& r) : val(r) + { + if (checking::is_nan(r)) + throw interval_lib::comparison_error(); + } + + const T& val; + }; + + typedef typename Policies::checking checking; + typedef typename Policies::rounding rounding; + + T low; + T up; +}; + +template inline +interval::interval(const T& v): low(v), up(v) +{ + if (checking::is_nan(v)) set_empty(); +} + +template inline +interval::interval(const T& l, const T& u): low(l), up(u) +{ + if (checking::is_nan(l) || checking::is_nan(u) || !(l <= u)) + set_empty(); +} + +template inline +interval& interval::operator=(const T& x) +{ + if (checking::is_nan(x)) set_empty(); + else low = up = x; + return *this; +} + +template inline +void interval::assign(const T& l, const T& u) +{ + if (checking::is_nan(l) || checking::is_nan(u) || !(l <= u)) + set_empty(); + else set(l, u); +} + +template inline +void interval::set(const T& l, const T& u) +{ + low = l; + up = u; +} + +template inline +void interval::set_empty() +{ + low = checking::empty_lower(); + up = checking::empty_upper(); +} + +template inline +void interval::set_whole() +{ + const T& inf = checking::inf(); + low = -inf; + up = inf; +} + +template inline +interval interval::hull(const T& x, const T& y) +{ + bool bad_x = checking::is_nan(x); + bool bad_y = checking::is_nan(y); + if (bad_x) + if (bad_y) return interval::empty(); + else return interval(y, y, true); + else + if (bad_y) return interval(x, x, true); + if (x < y) return interval(x, y, true); + else return interval(y, x, true); +} + +template inline +interval interval::empty() +{ + return interval(checking::empty_lower(), + checking::empty_upper(), true); +} + +template inline +interval interval::whole() +{ + const T& inf = checking::inf(); + return interval(-inf, inf, true); +} + +template inline +const T& interval::lower() const +{ + return low; +} + +template inline +const T& interval::upper() const +{ + return up; +} + +/* + * interval/interval comparisons + */ + +template inline +bool interval::operator< (const interval_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up < r.low) return true; + else if (low >= r.up) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator> (const interval_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (low > r.up) return true; + else if (up <= r.low) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator<= (const interval_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up <= r.low) return true; + else if (low > r.up) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator>= (const interval_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (low >= r.up) return true; + else if (up < r.low) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator== (const interval_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up == r.low && low == r.up) return true; + else if (up < r.low || low > r.up) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator!= (const interval_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up < r.low || low > r.up) return true; + else if (up == r.low && low == r.up) return false; + } + throw interval_lib::comparison_error(); +} + +/* + * interval/number comparisons + */ + +template inline +bool interval::operator< (const number_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up < r.val) return true; + else if (low >= r.val) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator> (const number_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (low > r.val) return true; + else if (up <= r.val) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator<= (const number_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up <= r.val) return true; + else if (low > r.val) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator>= (const number_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (low >= r.val) return true; + else if (up < r.val) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator== (const number_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up == r.val && low == r.val) return true; + else if (up < r.val || low > r.val) return false; + } + throw interval_lib::comparison_error(); +} + +template inline +bool interval::operator!= (const number_holder& r) const +{ + if (!checking::is_empty(low, up)) { + if (up < r.val || low > r.val) return true; + else if (up == r.val && low == r.val) return false; + } + throw interval_lib::comparison_error(); +} + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_INTERVAL_HPP diff --git a/include/boost/numeric/interval/io.hpp b/include/boost/numeric/interval/io.hpp new file mode 100644 index 0000000..3826bca --- /dev/null +++ b/include/boost/numeric/interval/io.hpp @@ -0,0 +1,121 @@ +/* Boost interval/io.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_IO_HPP +#define BOOST_NUMERIC_INTERVAL_IO_HPP + +#include +#include // for non-explicit string constructor +// #include +#include + + +namespace boost { +namespace numeric { + +/* + * Input/Output + */ + +#if 0 +template +std::basic_ostream& +operator<<(std::basic_ostream& os, const interval& r) +{ + typename std::basic_ostream::sentry sentry(os); + if (sentry) { + T l = r.lower(), u = r.upper(); + os << '[' << l << ',' << u << ']'; + } + return os; +} +#else +// Ensure generation of outer bounds in all cases +template inline +std::ostream& operator<<(std::ostream& os, const interval& r) +{ + std::streamsize p = os.precision(); // decimal notation + // FIXME poor man's power of 10, only up to 1E-15 + p = (p > 15) ? 15 : p - 1; + double eps = 1.0; while (p>0) { eps /= 10.0; --p; } + // widen the interval so output is correct + interval r_wide = widen(r, static_cast(eps / 2.0)); + os << "[" << r_wide.lower() << "," << r_wide.upper() << "]"; + return os; +} +#endif + + +#if 0 // FIXME: make the code work for g++ 3.* +template +std::basic_istream& +operator>>(std::basic_istream& is, const interval& r) +{ + T l = 0, u = 0; + std::locale loc = is.getloc(); + const std::ctype c_type = std::use_facet >(loc); + const char punct[] = "[,]"; + Ch wpunct[3]; + c_type.widen(punct, punct+3, wpunct); + Ch c; + is >> c; + if (ChTr::eq(c, wpunct[0])) { + is >> l >> c; + if (ChTr::eq(c, wpunct[1])) + is >> u >> c; + if (!ChTr::eq(c, wpunct[2])) + is.setstate(is.failbit); + } else { + is.putback(c); + is >> l; + u = l; + } + if (is) + r = succ(interval(l, r)); // round outward by 1 ulp + return is; +} +#else +template inline +std::istream& operator>>(std::istream& is, interval& r) +{ + T l, u; + char c = 0; + // we assume that l and u are representable numbers + if(is.peek() == '[') { + is >> c; + is >> l; + if(is.peek() == ',') { + is >> c; + is >> u; + is >> c; + } + } + if (c != ']') + is.setstate(is.failbit); + else if (is.good()) + r.assign(l, u); + return is; +} +#endif + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_IO_HPP diff --git a/include/boost/numeric/interval/limits.hpp b/include/boost/numeric/interval/limits.hpp new file mode 100644 index 0000000..15fea81 --- /dev/null +++ b/include/boost/numeric/interval/limits.hpp @@ -0,0 +1,41 @@ +#ifndef BOOST_NUMERIC_INTERVAL_LIMITS_HPP +#define BOOST_NUMERIC_INTERVAL_LIMITS_HPP + +#ifndef BOOST_NO_TEMPLATE_PARTIAL_SPECIALIZATION + +#include +#include +#include + +namespace std { + +template +class numeric_limits > + : public numeric_limits +{ +private: + typedef boost::numeric::interval I; + typedef numeric_limits bl; +public: + static I min() throw() { return I(bl::min(), bl::min()); } + static I max() throw() { return I(bl::max(), bl::max()); } + static I epsilon() throw() { return I(bl::epsilon(), bl::epsilon()); } + + BOOST_STATIC_CONSTANT(float_round_style, round_style = round_indeterminate); + BOOST_STATIC_CONSTANT(bool, is_iec559 = false); + + static I infinity () throw() { return I::whole(); } + static I quiet_NaN() throw() { return I::empty(); } + static I signaling_NaN() throw() + { return I(bl::signaling_NaN(), bl::signaling_Nan()); } + static I denorm_min() throw() + { return I(bl::denorm_min(), bl::denorm_min()); } +private: + static I round_error(); // hide this on purpose, not yet implemented +}; + +} // namespace std + +#endif + +#endif // BOOST_NUMERIC_INTERVAL_LIMITS_HPP diff --git a/include/boost/numeric/interval/policies.hpp b/include/boost/numeric/interval/policies.hpp new file mode 100644 index 0000000..dc46c6d --- /dev/null +++ b/include/boost/numeric/interval/policies.hpp @@ -0,0 +1,66 @@ +#ifndef BOOST_NUMERIC_INTERVAL_POLICIES_HPP +#define BOOST_NUMERIC_INTERVAL_POLICIES_HPP + +#include + +namespace boost { +namespace numeric { +namespace interval_lib { + +/* + * policies class + */ + +template +struct policies +{ + typedef Rounding rounding; + typedef Checking checking; +}; + +/* + * policies switching classes + */ + +template +class change_rounding +{ + typedef typename OldInterval::base_type T; + typedef typename OldInterval::traits_type p; + typedef typename p::checking checking; +public: + typedef interval > type; +}; + +template +class change_checking +{ + typedef typename OldInterval::base_type T; + typedef typename OldInterval::traits_type p; + typedef typename p::rounding rounding; +public: + typedef interval > type; +}; + +/* + * Protect / unprotect: control whether the rounding mode is set/reset + * at each operation, rather than once and for all. + */ + +template +class unprotect +{ + typedef typename OldInterval::base_type T; + typedef typename OldInterval::traits_type p; + typedef typename p::rounding r; + typedef typename r::unprotected_rounding newRounding; +public: + typedef typename change_rounding::type type; +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + + +#endif // BOOST_NUMERIC_INTERVAL_POLICIES_HPP diff --git a/include/boost/numeric/interval/rounded_arith.hpp b/include/boost/numeric/interval/rounded_arith.hpp new file mode 100644 index 0000000..2cb4862 --- /dev/null +++ b/include/boost/numeric/interval/rounded_arith.hpp @@ -0,0 +1,121 @@ +/* Boost interval/rounded_arith.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_ROUNDED_ARITH_HPP +#define BOOST_NUMERIC_INTERVAL_ROUNDED_ARITH_HPP + +#include +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { + +/* + * Three classes of rounding: exact, std, opp + * See documentation for details. + */ + +template +struct rounded_arith_exact: Rounding { + void init() { } + T add_down (const T& x, const T& y) { return x + y; } + T add_up (const T& x, const T& y) { return x + y; } + T sub_down (const T& x, const T& y) { return x - y; } + T sub_up (const T& x, const T& y) { return x - y; } + T mul_down (const T& x, const T& y) { return x * y; } + T mul_up (const T& x, const T& y) { return x * y; } + T div_down (const T& x, const T& y) { return x / y; } + T div_up (const T& x, const T& y) { return x / y; } + T median (const T& x, const T& y) { return (x + y) / 2; } + T sqrt_down(const T& x) + { BOOST_NUMERIC_INTERVAL_using_math(sqrt); return sqrt(x); } + T sqrt_up (const T& x) + { BOOST_NUMERIC_INTERVAL_using_math(sqrt); return sqrt(x); } + T int_down (const T& x) + { BOOST_NUMERIC_INTERVAL_using_math(floor); return floor(x); } + T int_up (const T& x) + { BOOST_NUMERIC_INTERVAL_using_math(ceil); return ceil(x); } +}; + +template +struct rounded_arith_std: Rounding { + #define BOOST_DN(EXPR) (downward(), force_rounding( EXPR )) + #define BOOST_NR(EXPR) (to_nearest(), force_rounding( EXPR )) + #define BOOST_UP(EXPR) (upward(), force_rounding( EXPR )) + void init() { } + T add_down (const T& x, const T& y) { return BOOST_DN( x + y ); } + T add_up (const T& x, const T& y) { return BOOST_UP( x + y ); } + T sub_down (const T& x, const T& y) { return BOOST_DN( x - y ); } + T sub_up (const T& x, const T& y) { return BOOST_UP( x - y ); } + T mul_down (const T& x, const T& y) { return BOOST_DN( x * y ); } + T mul_up (const T& x, const T& y) { return BOOST_UP( x * y ); } + T div_down (const T& x, const T& y) { return BOOST_DN( x / y ); } + T div_up (const T& x, const T& y) { return BOOST_UP( x / y ); } + T median (const T& x, const T& y) { return BOOST_NR( (x+y)/2 ); } + T sqrt_down(const T& x) { BOOST_NUMERIC_INTERVAL_using_math(sqrt); + return BOOST_DN( sqrt(x) ); } + T sqrt_up (const T& x) { BOOST_NUMERIC_INTERVAL_using_math(sqrt); + return BOOST_UP( sqrt(x) ); } + T int_down (const T& x) { return BOOST_DN( to_int(x) ); } + T int_up (const T& x) { return BOOST_UP( to_int(x) ); } + #undef BOOST_DN + #undef BOOST_NR + #undef BOOST_UP +}; + +template +struct rounded_arith_opp: Rounding { + void init() { upward(); } + #define BOOST_DN(EXPR) (downward(), force_rounding( EXPR )) + #define BOOST_NR(EXPR) (to_nearest(), force_rounding( EXPR )) + #define BOOST_UP(EXPR) force_rounding( EXPR ) + T add_down (const T& x, const T& y) { return -BOOST_UP( (-x) - y ); } + T add_up (const T& x, const T& y) { return BOOST_UP( x + y ); } + T sub_down (const T& x, const T& y) { return -BOOST_UP( y - x ); } + T sub_up (const T& x, const T& y) { return BOOST_UP( x - y ); } + T mul_down (const T& x, const T& y) { return -BOOST_UP( x * (-y) ); } + T mul_up (const T& x, const T& y) { return BOOST_UP( x * y ); } + T div_down (const T& x, const T& y) { return -BOOST_UP( x / (-y) ); } + T div_up (const T& x, const T& y) { return BOOST_UP( x / y ); } + T median (const T& x, const T& y) { T r = BOOST_NR( (x + y) / 2 ); + upward(); + return r; + } + T sqrt_down(const T& x) { BOOST_NUMERIC_INTERVAL_using_math(sqrt); + T r = BOOST_DN ( sqrt(x) ); + upward(); + return r; + } + T sqrt_up (const T& x) { BOOST_NUMERIC_INTERVAL_using_math(sqrt); + return BOOST_UP( sqrt(x) ); } + T int_down (const T& x) { return -to_int(-x); } + T int_up (const T& x) { return to_int(x); } + #undef BOOST_DN + #undef BOOST_NR + #undef BOOST_UP +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_ROUNDED_ARITH_HPP diff --git a/include/boost/numeric/interval/rounded_transc.hpp b/include/boost/numeric/interval/rounded_transc.hpp new file mode 100644 index 0000000..bfd20e8 --- /dev/null +++ b/include/boost/numeric/interval/rounded_transc.hpp @@ -0,0 +1,139 @@ +/* Boost interval/rounded_transc.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_ROUNDED_TRANSC_HPP +#define BOOST_NUMERIC_INTERVAL_ROUNDED_TRANSC_HPP + +#include +#include +#include + +namespace boost { +namespace numeric { +namespace interval_lib { + +template +struct rounded_transc_exact: Rounding +{ +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) { BOOST_NUMERIC_INTERVAL_using_math(f); return f(x); } \ + T f##_up (const T& x) { BOOST_NUMERIC_INTERVAL_using_math(f); return f(x); } + BOOST_NUMERIC_INTERVAL_new_func(exp) + BOOST_NUMERIC_INTERVAL_new_func(log) + BOOST_NUMERIC_INTERVAL_new_func(sin) + BOOST_NUMERIC_INTERVAL_new_func(cos) + BOOST_NUMERIC_INTERVAL_new_func(tan) + BOOST_NUMERIC_INTERVAL_new_func(asin) + BOOST_NUMERIC_INTERVAL_new_func(acos) + BOOST_NUMERIC_INTERVAL_new_func(atan) + BOOST_NUMERIC_INTERVAL_new_func(sinh) + BOOST_NUMERIC_INTERVAL_new_func(cosh) + BOOST_NUMERIC_INTERVAL_new_func(tanh) +# undef BOOST_NUMERIC_INTERVAL_new_func +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) { BOOST_NUMERIC_INTERVAL_using_ahyp(f); return f(x); } \ + T f##_up (const T& x) { BOOST_NUMERIC_INTERVAL_using_ahyp(f); return f(x); } + BOOST_NUMERIC_INTERVAL_new_func(asinh) + BOOST_NUMERIC_INTERVAL_new_func(acosh) + BOOST_NUMERIC_INTERVAL_new_func(atanh) +# undef BOOST_NUMERIC_INTERVAL_new_func +}; + +template +struct rounded_transc_std: Rounding +{ +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_math(f); downward(); return force_rounding(f(x)); } \ + T f##_up (const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_math(f); upward(); return force_rounding(f(x)); } + BOOST_NUMERIC_INTERVAL_new_func(exp) + BOOST_NUMERIC_INTERVAL_new_func(log) + BOOST_NUMERIC_INTERVAL_new_func(sin) + BOOST_NUMERIC_INTERVAL_new_func(cos) + BOOST_NUMERIC_INTERVAL_new_func(tan) + BOOST_NUMERIC_INTERVAL_new_func(asin) + BOOST_NUMERIC_INTERVAL_new_func(acos) + BOOST_NUMERIC_INTERVAL_new_func(atan) + BOOST_NUMERIC_INTERVAL_new_func(sinh) + BOOST_NUMERIC_INTERVAL_new_func(cosh) + BOOST_NUMERIC_INTERVAL_new_func(tanh) +# undef BOOST_NUMERIC_INTERVAL_new_func +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_ahyp(f); downward(); return force_rounding(f(x)); } \ + T f##_up (const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_ahyp(f); upward(); return force_rounding(f(x)); } + BOOST_NUMERIC_INTERVAL_new_func(asinh) + BOOST_NUMERIC_INTERVAL_new_func(acosh) + BOOST_NUMERIC_INTERVAL_new_func(atanh) +# undef BOOST_NUMERIC_INTERVAL_new_func +}; + +template +struct rounded_transc_opp: Rounding +{ +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_math(f); \ + downward(); T y = force_rounding(f(x)); upward(); return y; } \ + T f##_up (const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_math(f); return force_rounding(f(x)); } + BOOST_NUMERIC_INTERVAL_new_func(exp) + BOOST_NUMERIC_INTERVAL_new_func(log) + BOOST_NUMERIC_INTERVAL_new_func(cos) + BOOST_NUMERIC_INTERVAL_new_func(acos) + BOOST_NUMERIC_INTERVAL_new_func(cosh) +# undef BOOST_NUMERIC_INTERVAL_new_func +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_math(f); return -force_rounding(-f(x)); } \ + T f##_up (const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_math(f); return force_rounding(f(x)); } + BOOST_NUMERIC_INTERVAL_new_func(sin) + BOOST_NUMERIC_INTERVAL_new_func(tan) + BOOST_NUMERIC_INTERVAL_new_func(asin) + BOOST_NUMERIC_INTERVAL_new_func(atan) + BOOST_NUMERIC_INTERVAL_new_func(sinh) + BOOST_NUMERIC_INTERVAL_new_func(tanh) +# undef BOOST_NUMERIC_INTERVAL_new_func +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_ahyp(f); \ + downward(); T y = force_rounding(f(x)); upward(); return y; } \ + T f##_up (const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_ahyp(f); return force_rounding(f(x)); } + BOOST_NUMERIC_INTERVAL_new_func(asinh) + BOOST_NUMERIC_INTERVAL_new_func(atanh) +# undef BOOST_NUMERIC_INTERVAL_new_func +# define BOOST_NUMERIC_INTERVAL_new_func(f) \ + T f##_down(const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_ahyp(f); return -force_rounding(-f(x)); } \ + T f##_up (const T& x) \ + { BOOST_NUMERIC_INTERVAL_using_ahyp(f); return force_rounding(f(x)); } + BOOST_NUMERIC_INTERVAL_new_func(acosh) +# undef BOOST_NUMERIC_INTERVAL_new_func +}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_ROUNDED_TRANSC_HPP diff --git a/include/boost/numeric/interval/rounding.hpp b/include/boost/numeric/interval/rounding.hpp new file mode 100644 index 0000000..bc4e51b --- /dev/null +++ b/include/boost/numeric/interval/rounding.hpp @@ -0,0 +1,112 @@ +/* Boost interval/rounding.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_ROUNDING_HPP +#define BOOST_NUMERIC_INTERVAL_ROUNDING_HPP + +namespace boost { +namespace numeric { +namespace interval_lib { + +/* + * Default rounding_control class (does nothing) + */ + +template +struct rounding_control +{ + typedef int rounding_mode; + static void get_rounding_mode(rounding_mode&) {} + static void set_rounding_mode(rounding_mode) {} + static void upward() {} + static void downward() {} + static void to_nearest() {} + static const T& to_int(const T& x) { return x; } + static const T& force_rounding(const T& x) { return x; } +}; + +/* + * A few rounding control classes (exact/std/opp: see documentation) + * rounded_arith_* control the rounding of the arithmetic operators + * rounded_transc_* control the rounding of the transcendental functions + */ + +template > +struct rounded_arith_exact; + +template > +struct rounded_arith_std; + +template > +struct rounded_arith_opp; + +template +struct rounded_transc_dummy; + +template > +struct rounded_transc_exact; + +template > +struct rounded_transc_std; + +template > +struct rounded_transc_opp; + +/* + * State-saving classes: allow to set and reset rounding control + */ + +namespace detail { + +template +struct save_state_unprotected: Rounding +{ + typedef save_state_unprotected unprotected_rounding; +}; + +} // namespace detail + +template +struct save_state: Rounding +{ + typename Rounding::rounding_mode mode; + save_state() { + get_rounding_mode(mode); + init(); + } + ~save_state() { set_rounding_mode(mode); } + typedef detail::save_state_unprotected unprotected_rounding; +}; + +template +struct save_state_nothing: Rounding +{ + typedef save_state_nothing unprotected_rounding; +}; + +template +struct rounded_math: save_state_nothing > +{}; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_ROUNDING_HPP diff --git a/include/boost/numeric/interval/transc.hpp b/include/boost/numeric/interval/transc.hpp new file mode 100644 index 0000000..ff5eeab --- /dev/null +++ b/include/boost/numeric/interval/transc.hpp @@ -0,0 +1,241 @@ +/* Boost interval/transc.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_TRANSC_HPP +#define BOOST_NUMERIC_INTERVAL_TRANSC_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace numeric { + +template inline +interval exp(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + return I(rnd.exp_down(x.lower()), rnd.exp_up (x.upper()), true); +} + +template inline +interval log(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x) || x.upper() <= static_cast(0)) + return I::empty(); + typename Policies::rounding rnd; + typedef typename Policies::checking checking; + T l = (x.lower() <= static_cast(0)) + ? -checking::inf() : rnd.log_down(x.lower()); + return I(l, rnd.log_up(x.upper()), true); +} + +template inline +interval cos(const interval& x) +{ + if (interval_lib::detail::test_input(x)) + return interval::empty(); + typename Policies::rounding rnd; + typedef interval I; + typedef typename interval_lib::unprotect::type R; + + // get lower bound within [0, pi] + const R pi2 = interval_lib::pi_twice(); + R tmp = fmod((const R&)x, pi2); + if (width(tmp) >= pi2.lower()) + return I(-1, 1, true); // we are covering a full period + if (tmp.lower() >= interval_lib::constants::pi_upper()) + return -cos(tmp - interval_lib::pi()); + T l = tmp.lower(); + T u = tmp.upper(); + + BOOST_NUMERIC_INTERVAL_using_max(min); + + // separate into monotone subintervals + if (u <= interval_lib::constants::pi_lower()) + return I(rnd.cos_down(u), rnd.cos_up(l), true); + else if (u <= pi2.lower()) + return I(-1, rnd.cos_up(min(rnd.sub_down(pi2.lower(), u), l)), true); + else + return I(-1, 1, true); +} + +template inline +interval sin(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect::type R; + I r = cos((const R&)x - interval_lib::pi_half()); + (void)&rnd; + return r; +} + +template inline +interval tan(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + typedef typename interval_lib::unprotect::type R; + + // get lower bound within [-pi/2, pi/2] + const R pi = interval_lib::pi(); + R tmp = fmod((const I&)x, pi); + const T pi_half_d = interval_lib::constants::pi_half_lower(); + if (tmp.lower() >= pi_half_d) + tmp -= pi; + if (tmp.lower() <= -pi_half_d || tmp.upper() >= pi_half_d) + return I::whole(); + return I(rnd.tan_down(tmp.lower()), rnd.tan_up (tmp.upper()), true); +} + +template inline +interval asin(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x) + || x.upper() < static_cast(-1) || x.lower() > static_cast(1)) + return I::empty(); + typename Policies::rounding rnd; + T l = (x.lower() <= static_cast(-1)) + ? -interval_lib::constants::pi_half_upper() + : rnd.asin_down(x.lower()); + T u = (x.upper() >= static_cast(1) ) + ? interval_lib::constants::pi_half_upper() + : rnd.asin_up (x.upper()); + return I(l, u, true); +} + +template inline +interval acos(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x) + || x.upper() < static_cast(-1) || x.lower() > static_cast(1)) + return I::empty(); + typename Policies::rounding rnd; + T l = (x.upper() >= static_cast(1) ) + ? 0 + : rnd.acos_down(x.upper()); + T u = (x.lower() <= static_cast(-1)) + ? interval_lib::constants::pi_upper() + : rnd.acos_up (x.lower()); + return I(l, u, true); +} + +template inline +interval atan(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + return I(rnd.atan_down(x.lower()), rnd.atan_up(x.upper()), true); +} + +template inline +interval sinh(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + return I(rnd.sinh_down(x.lower()), rnd.sinh_up(x.upper()), true); +} + +template inline +interval cosh(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + if (interval_lib::detail::is_neg(x.upper())) + return I(rnd.cosh_down(x.upper()), rnd.cosh_up(x.lower()), true); + else if (!interval_lib::detail::is_neg(x.lower())) + return I(rnd.cosh_down(x.lower()), rnd.cosh_up(x.upper()), true); + else + return I(0, rnd.cosh_up(-x.lower() > x.upper() ? x.lower() : x.upper()), true); +} + +template inline +interval tanh(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + return I(rnd.tanh_down(x.lower()), rnd.tanh_up(x.upper()), true); +} + +template inline +interval asinh(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + typename Policies::rounding rnd; + return I(rnd.asinh_down(x.lower()), rnd.asinh_up(x.upper()), true); +} + +template inline +interval acosh(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x) || x.upper() < static_cast(1)) + return I::empty(); + typename Policies::rounding rnd; + T l = x.lower() <= static_cast(1) ? 0 : rnd.acosh_down(x.lower()); + return I(l, rnd.acosh_up(x.upper()), true); +} + +template inline +interval atanh(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x) + || x.upper() < static_cast(-1) || x.lower() > static_cast(1)) + return I::empty(); + typename Policies::rounding rnd; + typedef typename Policies::checking checking; + T l = (x.lower() <= static_cast(-1)) + ? -checking::inf() : rnd.atanh_down(x.lower()); + T u = (x.upper() >= static_cast(1) ) + ? checking::inf() : rnd.atanh_up (x.upper()); + return I(l, u, true); +} + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_TRANSC_HPP diff --git a/include/boost/numeric/interval/utility.hpp b/include/boost/numeric/interval/utility.hpp new file mode 100644 index 0000000..a01aaae --- /dev/null +++ b/include/boost/numeric/interval/utility.hpp @@ -0,0 +1,326 @@ +/* Boost interval/utility.hpp template implementation file + * + * Copyright Jens Maurer 2000 + * Copyright Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion 2002 + * Permission to use, copy, modify, sell, and distribute this software + * is hereby granted without fee provided that the above copyright notice + * appears in all copies and that both that copyright notice and this + * permission notice appear in supporting documentation, + * + * None of the above authors nor Polytechnic University make any + * representation about the suitability of this software for any + * purpose. It is provided "as is" without express or implied warranty. + * + * $Id$ + * + * Revision history: + * 2002-08-31 Prepared for boost formal review + * 2000-09-24 Separated from interval.hpp + */ + +#ifndef BOOST_NUMERIC_INTERVAL_UTILITY_HPP +#define BOOST_NUMERIC_INTERVAL_UTILITY_HPP + +#include +#include +#include +#include + +/* + * Implementation of simple functions + */ + +namespace boost { +namespace numeric { + +/* + * Utility Functions + */ + +template inline +const T& lower(const interval& x) +{ + return x.lower(); +} + +template inline +const T& upper(const interval& x) +{ + return x.upper(); +} + +template inline +T checked_lower(const interval& x) +{ + if (empty(x)) { + typedef typename Policies::checking checking; + return checking::nan(); + } + return low; +} + +template inline +T checked_upper(const interval& x) +{ + if (empty(x)) { + typedef typename Policies::checking checking; + return checking::nan(); + } + return up; +} + +template inline +T width(const interval& x) +{ + if (interval_lib::detail::test_input(x)) return 0; + typename Policies::rounding rnd; + return rnd.sub_up(x.upper(), x.lower()); +} + +template inline +T median(const interval& x) +{ + if (interval_lib::detail::test_input(x)) { + typedef typename Policies::checking checking; + return checking::nan(); + } + typename Policies::rounding rnd; + return rnd.median(x.lower(), x.upper()); +} + +template inline +interval widen(const interval& x, const T& v) +{ + if (interval_lib::detail::test_input(x)) + return interval::empty(); + typename Policies::rounding rnd; + return interval(rnd.sub_down(x.lower(), v), + rnd.add_up (x.upper(), v), true); +} + +/* + * Set-like operations + */ + +template inline +bool empty(const interval& x) +{ + return interval_lib::detail::test_input(x); +} + +template inline +bool in_zero(const interval& x) +{ + if (interval_lib::detail::test_input(x)) return false; + return x.lower() <= T(0) && T(0) <= x.upper(); +} + +template inline +bool in(const T& x, const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) return false; + return y.lower() <= x && x <= y.upper(); +} + +template inline +bool subset(const interval& x, + const interval& y) +{ + if (empty(x)) return true; + return !empty(y) && y.lower() <= x.lower() && x.upper() <= y.upper(); +} + +template inline +bool proper_subset(const interval& x, + const interval& y) +{ + if (empty(y)) return false; + if (empty(x)) return true; + return y.lower() <= x.lower() && x.upper() <= y.upper() && + (y.lower() != x.lower() || x.upper() != y.upper()); +} + +template inline +bool overlap(const interval& x, + const interval& y) +{ + if (interval_lib::detail::test_input(x, y)) return false; + return x.lower() <= y.lower() && y.lower() <= x.upper() || + y.lower() <= x.lower() && x.lower() <= y.upper(); +} + +template inline +bool singleton(const interval& x) +{ + return !empty(x) && x.lower() == x.upper(); +} + +template inline +bool equal(const interval& x, const interval& y) +{ + if (empty(x)) return empty(y); + return !empty(y) && x.lower() == y.lower() && x.upper() == y.upper(); +} + +template inline +interval intersect(const interval& x, + const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + BOOST_NUMERIC_INTERVAL_using_max(max); + if (interval_lib::detail::test_input(x, y)) + return interval::empty(); + const T& l = max(x.lower(), y.lower()); + const T& u = min(x.upper(), y.upper()); + if (l <= u) return interval(l, u, true); + else return interval::empty(); +} + +template inline +interval hull(const interval& x, + const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + BOOST_NUMERIC_INTERVAL_using_max(max); + bool bad_x = interval_lib::detail::test_input(x); + bool bad_y = interval_lib::detail::test_input(y); + if (bad_x) + if (bad_y) return interval::empty(); + else return y; + else + if (bad_y) return x; + return interval(min(x.lower(), y.lower()), + max(x.upper(), y.upper()), true); +} + +template inline +interval hull(const interval& x, const T& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + BOOST_NUMERIC_INTERVAL_using_max(max); + bool bad_x = interval_lib::detail::test_input(x); + bool bad_y = interval_lib::detail::test_input(y); + if (bad_y) + if (bad_x) return interval::empty(); + else return x; + else + if (bad_x) return interval(y, y, true); + return interval(min(x.lower(), y), + max(x.upper(), y), true); +} + +template inline +interval hull(const T& x, const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + BOOST_NUMERIC_INTERVAL_using_max(max); + bool bad_x = interval_lib::detail::test_input(x); + bool bad_y = interval_lib::detail::test_input(y); + if (bad_x) + if (bad_y) return interval::empty(); + else return y; + else + if (bad_y) return interval(x, x, true); + return interval(min(x, y.lower()), + max(x, y.upper()), true); +} + +template inline +interval hull(const T& x, const T& y) +{ + return interval::hull(x, y); +} + +template inline +std::pair, interval > +bisect(const interval& x) +{ + typedef interval I; + if (interval_lib::detail::test_input(x)) + return std::pair(I::empty(), I::empty()); + const T m = median(x); + return std::pair(I(x.lower(), m, true), I(m, x.upper(), true)); +} + +/* + * Elementary functions + */ + +template inline +interval abs(const interval& x) +{ + BOOST_NUMERIC_INTERVAL_using_max(max); + typedef interval I; + if (interval_lib::detail::test_input(x)) + return I::empty(); + if (!interval_lib::detail::is_neg(x.lower())) return x; + if (interval_lib::detail::is_neg(x.upper())) return -x; + return I(0, max(-x.lower(), x.upper()), true); +} + +template inline +interval max(const interval& x, + const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(max); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + return I(max(x.lower(), y.lower()), max(x.upper(), y.upper()), true); +} + +template inline +interval max(const interval& x, const T& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(max); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + return I(max(x.lower(), y), max(x.upper(), y), true); +} + +template inline +interval max(const T& x, const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(max); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + return I(max(x, y.lower()), max(x, y.upper()), true); +} + +template inline +interval min(const interval& x, + const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + return I(min(x.lower(), y.lower()), min(x.upper(), y.upper()), true); +} + +template inline +interval min(const interval& x, const T& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + return I(min(x.lower(), y), min(x.upper(), y), true); +} + +template inline +interval min(const T& x, const interval& y) +{ + BOOST_NUMERIC_INTERVAL_using_max(min); + typedef interval I; + if (interval_lib::detail::test_input(x, y)) + return I::empty(); + return I(min(x, y.lower()), min(x, y.upper()), true); +} + +} // namespace numeric +} // namespace boost + +#endif // BOOST_NUMERIC_INTERVAL_UTILITY_HPP diff --git a/test/Jamfile b/test/Jamfile new file mode 100644 index 0000000..cd7d0e1 --- /dev/null +++ b/test/Jamfile @@ -0,0 +1,67 @@ +subproject libs/numeric/interval/test ; + +unit-test add + : add.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test cmp + : cmp.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test cmp_exp + : cmp_exp.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test cmp_lex + : cmp_lex.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test det + : det.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test fmod + : fmod.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test mul + : mul.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test overflow + : overflow.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test pi + : pi.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test pow + : pow.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; + +unit-test test_float + : test_float.cpp + ../../../test/build/boost_test_exec_monitor + : $(BOOST_ROOT) + ; diff --git a/test/add.cpp b/test/add.cpp new file mode 100644 index 0000000..2082367 --- /dev/null +++ b/test/add.cpp @@ -0,0 +1,231 @@ +#include +#include +#include +#include +#include +#include +#include + +typedef enum { EXPR_VAR, EXPR_NEG, EXPR_UP, EXPR_DOWN, EXPR_ADD, EXPR_SUB } e_type; + +struct expr; +struct pexpr { + expr *ptr; + expr* operator->() { return ptr; } + pexpr(expr *p = NULL): ptr(p) { } +}; + +struct expr { + e_type type; + int var; + pexpr e; + pexpr e1, e2; +}; + +pexpr var(int v) { + pexpr e = new expr; + e->type = EXPR_VAR; + e->var = v; + return e; +} + +pexpr operator+(pexpr, pexpr); +pexpr operator-(pexpr, pexpr); +pexpr operator-(pexpr); + +pexpr operator+(pexpr a, pexpr b) { + if (a->type == EXPR_NEG) return b - a->e; + if (b->type == EXPR_NEG) return a - b->e; + if (a->type == EXPR_VAR && b->type == EXPR_VAR && a->var > b->var) return b + a; + pexpr c = new expr; + c->type = EXPR_ADD; + c->e1 = a; + c->e2 = b; + return c; +} + +pexpr operator-(pexpr a, pexpr b) { + if (b->type == EXPR_NEG) return a + b->e; + pexpr c = new expr; + c->type = EXPR_SUB; + c->e1 = a; + c->e2 = b; + return c; +} + +pexpr down(pexpr a) { + pexpr e = new expr; + e->type = EXPR_DOWN; + e->e = a; + return e; +} + +pexpr up(pexpr a) { + pexpr e = new expr; + e->type = EXPR_UP; + e->e = a; + return e; +} + +pexpr operator-(pexpr a) { + if (a->type == EXPR_NEG) return a->e; + if (a->type == EXPR_UP) return down(-a->e); + if (a->type == EXPR_DOWN) return up(-a->e); + if (a->type == EXPR_SUB) return a->e2 - a->e1; + if (a->type == EXPR_ADD) return -a->e1 - a->e2; + pexpr e = new expr; + e->type = EXPR_NEG; + e->e = a; + return e; +} + +bool operator==(pexpr a, pexpr b) { + if (a->type != b->type) return false; + if (a->type == EXPR_VAR) return a->var == b->var; + if (a->type == EXPR_DOWN || a->type == EXPR_UP || a->type == EXPR_NEG) + return a->e == b->e; + return a->e1 == b->e1 && a->e2 == b->e2; +} + +bool operator<=(pexpr, pexpr) { return true; } + +namespace boost { +namespace numeric { +namespace interval_lib { + +template<> +struct rounding_control { + typedef enum { RND_U, RND_M, RND_D } rounding_mode; + static rounding_mode mode; + rounding_control() { mode = RND_M; } + void get_rounding_mode(rounding_mode& m) { m = mode; } + void set_rounding_mode(rounding_mode m) { mode = m; } + void upward() { mode = RND_U; } + void downward() { mode = RND_D; } + pexpr force_rounding(pexpr a) { + switch (mode) { + case RND_U: return up(a); + case RND_D: return down(a); + default: throw "Unset rounding mode"; + } + } +}; + +rounding_control::rounding_mode rounding_control::mode = RND_M; + +} // namespace interval_lib +} // namespace numeric +} // namespace boost + +template +bool test_neg() { + I a(var(0), var(1)); + return equal(-a, I(-var(1), -var(0))); +} + +template +bool test_add() { + I a(var(0), var(1)), b(var(2), var(3)); + return equal(a + b, I(down(var(0) + var(2)), up(var(1) + var(3)))); +} + +template +bool test_add1() { + I a(var(0), var(1)); + return equal(a + var(2), I(down(var(0) + var(2)), up(var(1) + var(2)))); +} + +template +bool test_add2() { + I a(var(0), var(1)); + return equal(var(2) + a, I(down(var(0) + var(2)), up(var(1) + var(2)))); +} + +template +bool test_sub() { + I a(var(0), var(1)), b(var(2), var(3)); + return equal(a - b, I(down(var(0) - var(3)), up(var(1) - var(2)))); +} + +template +bool test_sub1() { + I a(var(0), var(1)); + return equal(a - var(2), I(down(var(0) - var(2)), up(var(1) - var(2)))); +} + +template +bool test_sub2() { + I a(var(0), var(1)); + return equal(var(2) - a, I(down(var(2) - var(1)), up(var(2) - var(0)))); +} + +template +bool test_addeq() { + I a(var(0), var(1)), b(var(2), var(3)); + return equal(a += b, I(down(var(0) + var(2)), up(var(1) + var(3)))); +} + +template +bool test_addeq1() { + I a(var(0), var(1)); + return equal(a += var(2), I(down(var(0) + var(2)), up(var(1) + var(2)))); +} + +template +bool test_subeq() { + I a(var(0), var(1)), b(var(2), var(3)); + return equal(a -= b, I(down(var(0) - var(3)), up(var(1) - var(2)))); +} + +template +bool test_subeq1() { + I a(var(0), var(1)); + return equal(a -= var(2), I(down(var(0) - var(2)), up(var(1) - var(2)))); +} + +struct my_checking +{ + static pexpr inf() { throw; } + static pexpr nan() { throw; } + static bool is_nan(const pexpr&) { return false; } + static pexpr empty_lower() { throw; } + static pexpr empty_upper() { throw; } + static bool is_empty(const pexpr&, const pexpr&) { return false; } +}; + +template +struct my_interval { +private: + typedef boost::numeric::interval_lib::save_state my_rounding; + typedef boost::numeric::interval_lib::policies my_policies; +public: + typedef boost::numeric::interval type; +}; + +int test_main(int, char *[]) { + typedef my_interval >::type I1; + typedef my_interval >::type I2; + BOOST_TEST((test_neg())); + BOOST_TEST((test_neg())); + BOOST_TEST((test_add())); + BOOST_TEST((test_add())); + BOOST_TEST((test_add1())); + BOOST_TEST((test_add1())); + BOOST_TEST((test_add2())); + BOOST_TEST((test_add2())); + BOOST_TEST((test_sub())); + BOOST_TEST((test_sub())); + BOOST_TEST((test_sub1())); + BOOST_TEST((test_sub1())); + BOOST_TEST((test_sub2())); + BOOST_TEST((test_sub2())); + BOOST_TEST((test_addeq())); + BOOST_TEST((test_addeq())); + BOOST_TEST((test_addeq1())); + BOOST_TEST((test_addeq1())); + BOOST_TEST((test_subeq())); + BOOST_TEST((test_subeq())); + BOOST_TEST((test_subeq1())); + BOOST_TEST((test_subeq1())); + return 0; +} diff --git a/test/cmp.cpp b/test/cmp.cpp new file mode 100644 index 0000000..df891d1 --- /dev/null +++ b/test/cmp.cpp @@ -0,0 +1,170 @@ +#include "cmp_header.hpp" + +// comparisons between [1,2] and [3,4] + +static void test_12_34() { + const I a(1,2), b(3,4); + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(b > a); + BOOST_CHECK(b >= a); + BOOST_CHECK(!(b < a)); + BOOST_CHECK(!(b <= a)); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,3] and [2,4] + +static void test_13_24() { + const I a(1,3), b(2,4); + + BOOST_C_EXN(a < b); + BOOST_C_EXN(a <= b); + BOOST_C_EXN(a > b); + BOOST_C_EXN(a >= b); + + BOOST_C_EXN(b < a); + BOOST_C_EXN(b <= a); + BOOST_C_EXN(b > a); + BOOST_C_EXN(b >= a); + + BOOST_C_EXN(a == b); + BOOST_C_EXN(a != b); +} + +// comparisons between [1,2] and [2,3] + +static void test_12_23() { + const I a(1,2), b(2,3); + + BOOST_C_EXN(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_C_EXN(a >= b); + + BOOST_CHECK(!(b < a)); + BOOST_C_EXN(b <= a); + BOOST_C_EXN(b > a); + BOOST_CHECK(b >= a); + + BOOST_C_EXN(a == b); + BOOST_C_EXN(a != b); +} + +static void test_12_E() { + I a(1, 2), b(I::empty()); + + BOOST_C_EXN(a < b); + BOOST_C_EXN(a <= b); + BOOST_C_EXN(a > b); + BOOST_C_EXN(a >= b); + + BOOST_C_EXN(b < a); + BOOST_C_EXN(b <= a); + BOOST_C_EXN(b > a); + BOOST_C_EXN(b >= a); + + BOOST_C_EXN(a == b); + BOOST_C_EXN(a != b); +} + +// comparisons between [1,2] and 0 + +static void test_12_0() { + const I a(1,2); + const int b = 0; + + BOOST_CHECK(!(a < b)); + BOOST_CHECK(!(a <= b)); + BOOST_CHECK(a > b); + BOOST_CHECK(a >= b); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,2] and 1 + +static void test_12_1() { + const I a(1,2); + const int b = 1; + + BOOST_CHECK(!(a < b)); + BOOST_C_EXN(a <= b); + BOOST_C_EXN(a > b); + BOOST_CHECK(a >= b); + + BOOST_C_EXN(a == b); + BOOST_C_EXN(a != b); +} + +// comparisons between [1,2] and 2 + +static void test_12_2() { + const I a(1,2); + const int b = 2; + + BOOST_C_EXN(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_C_EXN(a >= b); + + BOOST_C_EXN(a == b); + BOOST_C_EXN(a != b); +} + +// comparisons between [1,2] and 3 + +static void test_12_3() { + const I a(1,2); + const int b = 3; + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +static void test_12_12() { + const I a(1,2), b(1,2); + BOOST_C_EXN(a == b); + BOOST_C_EXN(a != b); +} + +static void test_11_11() { + const I a(1,1), b(1,1); + BOOST_CHECK(a == b); + BOOST_CHECK(!(a != b)); +} + +static void test_11_1() { + const I a(1,1); + const int b = 1; + BOOST_CHECK(a == b); + BOOST_CHECK(!(a != b)); +} + +int test_main(int, char *[]) { + test_12_34(); + test_13_24(); + test_12_23(); + test_12_E(); + test_12_0(); + test_12_1(); + test_12_2(); + test_12_3(); + test_12_12(); + test_11_11(); + test_11_1(); + + return 0; +} diff --git a/test/cmp_exp.cpp b/test/cmp_exp.cpp new file mode 100644 index 0000000..e1d1fdd --- /dev/null +++ b/test/cmp_exp.cpp @@ -0,0 +1,309 @@ +#include "cmp_header.hpp" + +// comparisons between [1,2] and [3,4] + +static void test_12_34() { + const I a(1,2), b(3,4); + + BOOST_CHECK(cerlt(a, b)); + BOOST_CHECK(cerle(a, b)); + BOOST_CHECK(!cergt(a, b)); + BOOST_CHECK(!cerge(a, b)); + + BOOST_CHECK(!cerlt(b, a)); + BOOST_CHECK(!cerle(b, a)); + BOOST_CHECK(cergt(b, a)); + BOOST_CHECK(cerge(b, a)); + + BOOST_CHECK(poslt(a, b)); + BOOST_CHECK(posle(a, b)); + BOOST_CHECK(!posgt(a, b)); + BOOST_CHECK(!posge(a, b)); + + BOOST_CHECK(!poslt(b, a)); + BOOST_CHECK(!posle(b, a)); + BOOST_CHECK(posgt(b, a)); + BOOST_CHECK(posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(cerne(a, b)); + BOOST_CHECK(!poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(cerne(b, a)); + BOOST_CHECK(!poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +// comparisons between [1,3] and [2,4] + +static void test_13_24() { + const I a(1,3), b(2,4); + + BOOST_CHECK(!cerlt(a, b)); + BOOST_CHECK(!cerle(a, b)); + BOOST_CHECK(!cergt(a, b)); + BOOST_CHECK(!cerge(a, b)); + + BOOST_CHECK(!cerlt(b, a)); + BOOST_CHECK(!cerle(b, a)); + BOOST_CHECK(!cergt(b, a)); + BOOST_CHECK(!cerge(b, a)); + + BOOST_CHECK(poslt(a, b)); + BOOST_CHECK(posle(a, b)); + BOOST_CHECK(posgt(a, b)); + BOOST_CHECK(posge(a, b)); + + BOOST_CHECK(poslt(b, a)); + BOOST_CHECK(posle(b, a)); + BOOST_CHECK(posgt(b, a)); + BOOST_CHECK(posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +// comparisons between [1,2] and [2,3] + +static void test_12_23() { + const I a(1,2), b(2,3); + + BOOST_CHECK(!cerlt(a, b)); + BOOST_CHECK(cerle(a, b)); + BOOST_CHECK(!cergt(a, b)); + BOOST_CHECK(!cerge(a, b)); + + BOOST_CHECK(!cerlt(b, a)); + BOOST_CHECK(!cerle(b, a)); + BOOST_CHECK(!cergt(b, a)); + BOOST_CHECK(cerge(b, a)); + + BOOST_CHECK(poslt(a, b)); + BOOST_CHECK(posle(a, b)); + BOOST_CHECK(!posgt(a, b)); + BOOST_CHECK(posge(a, b)); + + BOOST_CHECK(!poslt(b, a)); + BOOST_CHECK(posle(b, a)); + BOOST_CHECK(posgt(b, a)); + BOOST_CHECK(posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +// comparisons between [1,2] and 0 + +static void test_12_0() { + const I a(1,2); + const int b = 0; + + BOOST_CHECK(!cerlt(a, b)); + BOOST_CHECK(!cerle(a, b)); + BOOST_CHECK(cergt(a, b)); + BOOST_CHECK(cerge(a, b)); + + BOOST_CHECK(cerlt(b, a)); + BOOST_CHECK(cerle(b, a)); + BOOST_CHECK(!cergt(b, a)); + BOOST_CHECK(!cerge(b, a)); + + BOOST_CHECK(!poslt(a, b)); + BOOST_CHECK(!posle(a, b)); + BOOST_CHECK(posgt(a, b)); + BOOST_CHECK(posge(a, b)); + + BOOST_CHECK(poslt(b, a)); + BOOST_CHECK(posle(b, a)); + BOOST_CHECK(!posgt(b, a)); + BOOST_CHECK(!posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(cerne(a, b)); + BOOST_CHECK(!poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(cerne(b, a)); + BOOST_CHECK(!poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +// comparisons between [1,2] and 1 + +static void test_12_1() { + const I a(1,2); + const int b = 1; + + BOOST_CHECK(!cerlt(a, b)); + BOOST_CHECK(!cerle(a, b)); + BOOST_CHECK(!cergt(a, b)); + BOOST_CHECK(cerge(a, b)); + + BOOST_CHECK(!cerlt(b, a)); + BOOST_CHECK(cerle(b, a)); + BOOST_CHECK(!cergt(b, a)); + BOOST_CHECK(!cerge(b, a)); + + BOOST_CHECK(!poslt(a, b)); + BOOST_CHECK(posle(a, b)); + BOOST_CHECK(posgt(a, b)); + BOOST_CHECK(posge(a, b)); + + BOOST_CHECK(poslt(b, a)); + BOOST_CHECK(posle(b, a)); + BOOST_CHECK(!posgt(b, a)); + BOOST_CHECK(posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +// comparisons between [1,2] and 2 + +static void test_12_2() { + const I a(1,2); + const int b = 2; + + BOOST_CHECK(!cerlt(a, b)); + BOOST_CHECK(cerle(a, b)); + BOOST_CHECK(!cergt(a, b)); + BOOST_CHECK(!cerge(a, b)); + + BOOST_CHECK(!cerlt(b, a)); + BOOST_CHECK(!cerle(b, a)); + BOOST_CHECK(!cergt(b, a)); + BOOST_CHECK(cerge(b, a)); + + BOOST_CHECK(poslt(a, b)); + BOOST_CHECK(posle(a, b)); + BOOST_CHECK(!posgt(a, b)); + BOOST_CHECK(posge(a, b)); + + BOOST_CHECK(!poslt(b, a)); + BOOST_CHECK(posle(b, a)); + BOOST_CHECK(posgt(b, a)); + BOOST_CHECK(posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +// comparisons between [1,2] and 3 + +static void test_12_3() { + const I a(1,2); + const int b = 3; + + BOOST_CHECK(cerlt(a, b)); + BOOST_CHECK(cerle(a, b)); + BOOST_CHECK(!cergt(a, b)); + BOOST_CHECK(!cerge(a, b)); + + BOOST_CHECK(!cerlt(b, a)); + BOOST_CHECK(!cerle(b, a)); + BOOST_CHECK(cergt(b, a)); + BOOST_CHECK(cerge(b, a)); + + BOOST_CHECK(poslt(a, b)); + BOOST_CHECK(posle(a, b)); + BOOST_CHECK(!posgt(a, b)); + BOOST_CHECK(!posge(a, b)); + + BOOST_CHECK(!poslt(b, a)); + BOOST_CHECK(!posle(b, a)); + BOOST_CHECK(posgt(b, a)); + BOOST_CHECK(posge(b, a)); + + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(cerne(a, b)); + BOOST_CHECK(!poseq(a, b)); + BOOST_CHECK(posne(a, b)); + + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(cerne(b, a)); + BOOST_CHECK(!poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +static void test_12_12() { + const I a(1,2), b(1,2); + BOOST_CHECK(!cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(posne(a, b)); + BOOST_CHECK(!cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(posne(b, a)); +} + +static void test_11_11() { + const I a(1,1), b(1,1); + BOOST_CHECK(cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(!posne(a, b)); + BOOST_CHECK(cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(!posne(b, a)); +} + +static void test_11_1() { + const I a(1,1); + const int b = 1; + BOOST_CHECK(cereq(a, b)); + BOOST_CHECK(!cerne(a, b)); + BOOST_CHECK(poseq(a, b)); + BOOST_CHECK(!posne(a, b)); + BOOST_CHECK(cereq(b, a)); + BOOST_CHECK(!cerne(b, a)); + BOOST_CHECK(poseq(b, a)); + BOOST_CHECK(!posne(b, a)); +} + +int test_main(int, char *[]) { + test_12_34(); + test_13_24(); + test_12_23(); + test_12_0(); + test_12_1(); + test_12_2(); + test_12_3(); + test_12_12(); + test_11_11(); + test_11_1(); + + return 0; +} diff --git a/test/cmp_header.hpp b/test/cmp_header.hpp new file mode 100644 index 0000000..f80f857 --- /dev/null +++ b/test/cmp_header.hpp @@ -0,0 +1,16 @@ +#include +#include +#include +#include +#include + +struct empty_class {}; + +typedef boost::numeric::interval_lib::policies + > + my_policies; + +typedef boost::numeric::interval I; + +#define BOOST_C_EXN(e) \ + BOOST_CHECK_THROW(e, boost::numeric::interval_lib::comparison_error) diff --git a/test/cmp_lex.cpp b/test/cmp_lex.cpp new file mode 100644 index 0000000..5e5da9e --- /dev/null +++ b/test/cmp_lex.cpp @@ -0,0 +1,157 @@ +#include "cmp_header.hpp" + +using namespace boost::numeric::interval_lib::compare::lexicographic; + +// comparisons between [1,2] and [3,4] + +static void test_12_34() { + const I a(1,2), b(3,4); + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(b > a); + BOOST_CHECK(b >= a); + BOOST_CHECK(!(b < a)); + BOOST_CHECK(!(b <= a)); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,3] and [2,4] + +static void test_13_24() { + const I a(1,3), b(2,4); + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(!(b < a)); + BOOST_CHECK(!(b <= a)); + BOOST_CHECK(b > a); + BOOST_CHECK(b >= a); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,2] and [2,3] + +static void test_12_23() { + const I a(1,2), b(2,3); + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(!(b < a)); + BOOST_CHECK(!(b <= a)); + BOOST_CHECK(b > a); + BOOST_CHECK(b >= a); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,2] and 0 + +static void test_12_0() { + const I a(1,2); + const int b = 0; + + BOOST_CHECK(!(a < b)); + BOOST_CHECK(!(a <= b)); + BOOST_CHECK(a > b); + BOOST_CHECK(a >= b); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,2] and 1 + +static void test_12_1() { + const I a(1,2); + const int b = 1; + + BOOST_CHECK(!(a < b)); + BOOST_CHECK(!(a <= b)); + BOOST_CHECK(a > b); + BOOST_CHECK(a >= b); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,2] and 2 + +static void test_12_2() { + const I a(1,2); + const int b = 2; + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +// comparisons between [1,2] and 3 + +static void test_12_3() { + const I a(1,2); + const int b = 3; + + BOOST_CHECK(a < b); + BOOST_CHECK(a <= b); + BOOST_CHECK(!(a > b)); + BOOST_CHECK(!(a >= b)); + + BOOST_CHECK(!(a == b)); + BOOST_CHECK(a != b); +} + +static void test_12_12() { + const I a(1,2), b(1,2); + + BOOST_CHECK(a == b); + BOOST_CHECK(!(a != b)); +} + +static void test_11_11() { + const I a(1,1), b(1,1); + + BOOST_CHECK(a == b); + BOOST_CHECK(!(a != b)); +} + +static void test_11_1() { + const I a(1,1); + const int b = 1; + + BOOST_CHECK(a == b); + BOOST_CHECK(!(a != b)); +} + +int test_main(int, char *[]) { + test_12_34(); + test_13_24(); + test_12_23(); + test_12_0(); + test_12_1(); + test_12_2(); + test_12_3(); + test_12_12(); + test_11_11(); + test_11_1(); + + return 0; +} diff --git a/test/det.cpp b/test/det.cpp new file mode 100644 index 0000000..551b480 --- /dev/null +++ b/test/det.cpp @@ -0,0 +1,89 @@ +#include +#include + +#define size 8 + +template +void det(I (&mat)[size][size]) { + for(int i = 0; i < size; i++) + for(int j = 0; j < size; j++) + mat[i][j] = I(1) / I(i + j + 1); + + for(int i = 0; i < size - 1; i++) { + int m = i, n = i; + typename I::base_type v = 0; + for(int a = i; a < size; a++) + for(int b = i; b < size; b++) { + typename I::base_type w = abs(mat[a][b]).lower(); + if (w > v) { m = a; n = b; v = w; } + } + if (n != i) + for(int a = 0; a < size; a++) { + I t = mat[a][n]; + mat[a][n] = mat[a][i]; + mat[a][i] = t; + } + if (m != i) + for(int b = i; b < size; b++) { + I t = mat[m][b]; + mat[m][b] = mat[m][i]; + mat[m][i] = t; + } + if (((m + n) & 1) == 1) { }; + I c = mat[i][i]; + for(int j = i + 1; j < size; j++) { + I f = mat[j][i] / c; + for(int k = i; k < size; k++) + mat[j][k] -= f * mat[i][k]; + } + if (in_zero(c)) return; + } +} + +namespace my_namespace { + +using namespace boost; +using namespace numeric; +using namespace interval_lib; + +template +struct variants { + typedef interval I_op; + typedef typename change_rounding > >::type I_sp; + typedef typename unprotect::type I_ou; + typedef typename unprotect::type I_su; + typedef T type; +}; + +} + +template +bool test() { + typedef my_namespace::variants types; + types::I_op mat_op[size][size]; + types::I_sp mat_sp[size][size]; + types::I_ou mat_ou[size][size]; + types::I_su mat_su[size][size]; + det(mat_op); + det(mat_sp); + { types::I_op::traits_type::rounding rnd; det(mat_ou); } + { types::I_sp::traits_type::rounding rnd; det(mat_su); } + for(int i = 0; i < size; i++) + for(int j = 0; j < size; j++) { + typedef types::I_op I; + I d_op = mat_op[i][j]; + I d_sp = mat_sp[i][j]; + I d_ou = mat_ou[i][j]; + I d_su = mat_su[i][j]; + if (!(equal(d_op, d_sp) && equal(d_sp, d_ou) && equal(d_ou, d_su))) + return false; + } + return true; +} + +int test_main(int, char *[]) { + BOOST_TEST(test()); + BOOST_TEST(test()); + BOOST_TEST(test()); + return 0; +} diff --git a/test/fmod.cpp b/test/fmod.cpp new file mode 100644 index 0000000..b869618 --- /dev/null +++ b/test/fmod.cpp @@ -0,0 +1,42 @@ +#include +#include +#include +#include +#include +#include +#include + +struct my_rounded_arith { + int sub_down(int x, int y) { return x - y; } + int sub_up (int x, int y) { return x - y; } + int mul_down(int x, int y) { return x * y; } + int mul_up (int x, int y) { return x * y; } + int div_down(int x, int y) { + int q = x / y; + return (x % y < 0) ? (q - 1) : q; + } + int int_down(int x) { return x; } +}; + +using namespace boost; +using namespace numeric; +using namespace interval_lib; + +typedef change_rounding, save_state_nothing >::type I; + +int test_main(int, char *[]) { + + BOOST_CHECK(equal(fmod(I(6,9), 7), I(6,9))); + BOOST_CHECK(equal(fmod(6, I(7,8)), I(6,6))); + BOOST_CHECK(equal(fmod(I(6,9), I(7,8)), I(6,9))); + + BOOST_CHECK(equal(fmod(I(13,17), 7), I(6,10))); + BOOST_CHECK(equal(fmod(13, I(7,8)), I(5,6))); + BOOST_CHECK(equal(fmod(I(13,17), I(7,8)), I(5,10))); + + BOOST_CHECK(equal(fmod(I(-17,-13), 7), I(4,8))); + BOOST_CHECK(equal(fmod(-17, I(7,8)), I(4,7))); + BOOST_CHECK(equal(fmod(I(-17,-13), I(7,8)), I(4,11))); + + return 0; +} diff --git a/test/io.cpp b/test/io.cpp new file mode 100644 index 0000000..4199da8 --- /dev/null +++ b/test/io.cpp @@ -0,0 +1,95 @@ +#include +#include +#include + +#include +#include +#include +#include + +using namespace boost; +using namespace interval_lib; + + +namespace my_namespace { + +template +bool test(int precision) { + interval pi_i( 3.1415926, 3.1415927 ); + // write pi_i to a string + std::ostringstream oss; + oss.precision(precision); + oss << pi_i; + std::string pi_string = oss.str(); +#ifdef BOOST_NUMERIC_INTERVAL_IO_VERBOSE + std::cerr << "FYI: printing [3.1415926,3.1415927] with precision " + << precision << " as " + << std::setprecision(precision) << pi_i << std::endl; +#endif + // read pi_o from same string + interval pi_o; + std::istringstream iss(pi_string); + iss >> pi_o; + // the result should be at least as imprecise +#ifdef BOOST_NUMERIC_INTERVAL_IO_VERBOSE + if (!subset(pi_i, pi_o)) + std::cerr << "Error: getting [3.1415926,3.1415927] with precision " + << precision << " as " << std::setprecision(precision) + << "[" << pi_o.lower() << std::setprecision(precision) + << "," << pi_o.upper() << "]" << std::endl; + else + std::cerr << "Good: getting [3.1415926,3.1415927] with precision " + << precision << " as " << std::setprecision(precision) + << "[" << pi_o.lower() << std::setprecision(precision) + << "," << pi_o.upper() << "]" << std::endl; +#endif + return subset(pi_i, pi_o); +} + +} // namespace my_namespace + +int test_main(int, char *[]) { + using my_namespace::test; + BOOST_TEST(test(0)); + BOOST_TEST(test(1)); + BOOST_TEST(test(2)); + BOOST_TEST(test(3)); + BOOST_TEST(test(4)); + BOOST_TEST(test(5)); + BOOST_TEST(test(6)); + BOOST_TEST(test(7)); + BOOST_TEST(test(8)); + BOOST_TEST(test(0)); + BOOST_TEST(test(1)); + BOOST_TEST(test(2)); + BOOST_TEST(test(3)); + BOOST_TEST(test(4)); + BOOST_TEST(test(5)); + BOOST_TEST(test(6)); + BOOST_TEST(test(7)); + BOOST_TEST(test(8)); + BOOST_TEST(test(9)); + BOOST_TEST(test(10)); + BOOST_TEST(test(11)); + BOOST_TEST(test(12)); + BOOST_TEST(test(13)); + BOOST_TEST(test(14)); + BOOST_TEST(test(15)); + BOOST_TEST(test(0)); + BOOST_TEST(test(1)); + BOOST_TEST(test(2)); + BOOST_TEST(test(3)); + BOOST_TEST(test(4)); + BOOST_TEST(test(5)); + BOOST_TEST(test(6)); + BOOST_TEST(test(7)); + BOOST_TEST(test(8)); + BOOST_TEST(test(9)); + BOOST_TEST(test(10)); + BOOST_TEST(test(11)); + BOOST_TEST(test(12)); + BOOST_TEST(test(13)); + BOOST_TEST(test(14)); + BOOST_TEST(test(15)); + return 0; +} diff --git a/test/mul.cpp b/test/mul.cpp new file mode 100644 index 0000000..7dc5ae7 --- /dev/null +++ b/test/mul.cpp @@ -0,0 +1,111 @@ +#include +#include +#include + +typedef boost::numeric::interval I; + +static double min(double a, double b, double c, double d) { + return std::min(std::min(a, b), std::min(c, d)); +} + +static double max(double a, double b, double c, double d) { + return std::max(std::max(a, b), std::max(c, d)); +} + +static bool test_mul(double al, double au, double bl, double bu) { + I a(al, au), b(bl, bu); + I c = a * b; + return c.lower() == min(al*bl, al*bu, au*bl, au*bu) + && c.upper() == max(al*bl, al*bu, au*bl, au*bu); +} + +static bool test_mul1(double ac, double bl, double bu) { + I a(ac), b(bl, bu); + I c = ac * b; + I d = b * ac; + I e = a * b; + return equal(c, d) && equal(d, e); +} + +static bool test_div(double al, double au, double bl, double bu) { + I a(al, au), b(bl, bu); + I c = a / b; + return c.lower() == min(al/bl, al/bu, au/bl, au/bu) + && c.upper() == max(al/bl, al/bu, au/bl, au/bu); +} + +static bool test_div1(double al, double au, double bc) { + I a(al, au), b(bc); + I c = a / bc; + I d = a / b; + return equal(c, d); +} + +static bool test_div2(double ac, double bl, double bu) { + I a(ac), b(bl, bu); + I c = ac / b; + I d = a / b; + return equal(c, d); +} + +static bool test_square(double al, double au) { + I a(al, au); + I b = square(a); + I c = a * a; + return b.upper() == c.upper() && + (b.lower() == c.lower() || (c.lower() <= 0 && b.lower() == 0)); +} + +static bool test_sqrt(double al, double au) { + I a(al, au); + I b = square(sqrt(a)); + return subset(abs(a), b); +} + +int test_main(int, char*[]) { + BOOST_CHECK(test_mul(2, 3, 5, 7)); + BOOST_CHECK(test_mul(2, 3, -5, 7)); + BOOST_CHECK(test_mul(2, 3, -7, -5)); + BOOST_CHECK(test_mul(-2, 3, 5, 7)); + BOOST_CHECK(test_mul(-2, 3, -5, 7)); + BOOST_CHECK(test_mul(-2, 3, -7, -5)); + BOOST_CHECK(test_mul(-3, -2, 5, 7)); + BOOST_CHECK(test_mul(-3, -2, -5, 7)); + BOOST_CHECK(test_mul(-3, -2, -7, -5)); + + BOOST_CHECK(test_mul1(3, 5, 7)); + BOOST_CHECK(test_mul1(3, -5, 7)); + BOOST_CHECK(test_mul1(3, -7, -5)); + BOOST_CHECK(test_mul1(-3, 5, 7)); + BOOST_CHECK(test_mul1(-3, -5, 7)); + BOOST_CHECK(test_mul1(-3, -7, -5)); + + BOOST_CHECK(test_div(30, 42, 2, 3)); + BOOST_CHECK(test_div(30, 42, -3, -2)); + BOOST_CHECK(test_div(-30, 42, 2, 3)); + BOOST_CHECK(test_div(-30, 42, -3, -2)); + BOOST_CHECK(test_div(-42, -30, 2, 3)); + BOOST_CHECK(test_div(-42, -30, -3, -2)); + + BOOST_CHECK(test_div1(30, 42, 3)); + BOOST_CHECK(test_div1(30, 42, -3)); + BOOST_CHECK(test_div1(-30, 42, 3)); + BOOST_CHECK(test_div1(-30, 42, -3)); + BOOST_CHECK(test_div1(-42, -30, 3)); + BOOST_CHECK(test_div1(-42, -30, -3)); + + BOOST_CHECK(test_div2(30, 2, 3)); + BOOST_CHECK(test_div2(30, -3, -2)); + BOOST_CHECK(test_div2(-30, 2, 3)); + BOOST_CHECK(test_div2(-30, -3, -2)); + + BOOST_CHECK(test_square(2, 3)); + BOOST_CHECK(test_square(-2, 3)); + BOOST_CHECK(test_square(-3, 2)); + + BOOST_CHECK(test_sqrt(2, 3)); + BOOST_CHECK(test_sqrt(5, 7)); + BOOST_CHECK(test_sqrt(-1, 2)); + + return 0; +} diff --git a/test/overflow.cpp b/test/overflow.cpp new file mode 100644 index 0000000..1bc0f18 --- /dev/null +++ b/test/overflow.cpp @@ -0,0 +1,27 @@ +#include +#include + +template +void test_one(typename I::base_type x, typename I::base_type f) { + I y = x; + typename I::base_type g = 1 / f; + const int nb = 10000; + for(int i = 0; i < nb; i++) y *= f; + for(int i = 0; i < nb; i++) y *= g; + BOOST_TEST(in(x, y)); +} + +template +void test() { + test_one(1., 25.); + test_one(1., 0.04); + test_one(-1., 25.); + test_one(-1., 0.04); +} + +int test_main(int, char *[]) { + test >(); + test >(); + test >(); + return 0; +} diff --git a/test/pi.cpp b/test/pi.cpp new file mode 100644 index 0000000..b04183f --- /dev/null +++ b/test/pi.cpp @@ -0,0 +1,48 @@ +#include +#include +#include + +#define PI 3.14159265358979323846 + +typedef boost::numeric::interval I_i; +typedef boost::numeric::interval I_f; +typedef boost::numeric::interval I_d; +typedef boost::numeric::interval I_ld; + +using boost::numeric::interval_lib::pi; +using boost::numeric::interval_lib::pi_half; +using boost::numeric::interval_lib::pi_twice; + +int test_main(int, char *[]) { + I_i pi_i = pi(); + I_f pi_f = pi(); + I_d pi_d = pi(); + I_ld pi_ld = pi(); + + BOOST_TEST(in((int) PI, pi_i)); + BOOST_TEST(in((float) PI, pi_f)); + BOOST_TEST(in((double)PI, pi_d)); + BOOST_TEST(subset(pi_i, widen(I_i((int) PI), 1))); + BOOST_TEST(subset(pi_f, widen(I_f((float) PI), std::numeric_limits ::min()))); + BOOST_TEST(subset(pi_d, widen(I_d((double)PI), std::numeric_limits::min()))); + + // We can't test the following equalities for interval. + I_f pi_f_half = pi_half(); + I_f pi_f_twice = pi_twice(); + + I_d pi_d_half = pi_half(); + I_d pi_d_twice = pi_twice(); + + I_ld pi_ld_half = pi_half(); + I_ld pi_ld_twice = pi_twice(); + + BOOST_TEST(equal(2.0f * pi_f_half, pi_f)); + BOOST_TEST(equal(2.0 * pi_d_half, pi_d)); + BOOST_TEST(equal(2.0l * pi_ld_half, pi_ld)); + + BOOST_TEST(equal(2.0f * pi_f, pi_f_twice)); + BOOST_TEST(equal(2.0 * pi_d, pi_d_twice)); + BOOST_TEST(equal(2.0l * pi_ld, pi_ld_twice)); + + return 0; +} diff --git a/test/pow.cpp b/test/pow.cpp new file mode 100644 index 0000000..e01bd53 --- /dev/null +++ b/test/pow.cpp @@ -0,0 +1,28 @@ +#include +#include + +bool test_pow(double al, double au, double bl, double bu, int p) { + typedef boost::numeric::interval I; + I b = pow(I(al, au), p); + return b.lower() == bl && b.upper() == bu; +} + +int test_main(int, char *[]) { + BOOST_TEST(test_pow(2, 3, 8, 27, 3)); + BOOST_TEST(test_pow(2, 3, 16, 81, 4)); + BOOST_TEST(test_pow(-3, 2, -27, 8, 3)); + BOOST_TEST(test_pow(-3, 2, 0, 81, 4)); + BOOST_TEST(test_pow(-3, -2, -27, -8, 3)); + BOOST_TEST(test_pow(-3, -2, 16, 81, 4)); + + BOOST_TEST(test_pow(2, 4, 1./64, 1./8, -3)); + BOOST_TEST(test_pow(2, 4, 1./256, 1./16, -4)); + BOOST_TEST(test_pow(-4, -2, -1./8, -1./64, -3)); + BOOST_TEST(test_pow(-4, -2, 1./256, 1./16, -4)); + + BOOST_TEST(test_pow(2, 3, 1, 1, 0)); + BOOST_TEST(test_pow(-3, 2, 1, 1, 0)); + BOOST_TEST(test_pow(-3, -2, 1, 1, 0)); + + return 0; +} diff --git a/test/test_float.cpp b/test/test_float.cpp new file mode 100644 index 0000000..584bf09 --- /dev/null +++ b/test/test_float.cpp @@ -0,0 +1,103 @@ +#include +#include + +template +void test_unary() { + typedef typename F::I I; + for(I a(-10., -9.91); a.lower() <= 10.; a += 0.3) { + if (!F::validate(a)) continue; + I rI = F::f_I(a); + T rT1 = F::f_T(a.lower()), rT2 = F::f_T(a.upper()), + rT3 = F::f_T(median(a)); + BOOST_CHECK(in(rT1, rI)); + BOOST_CHECK(in(rT2, rI)); + BOOST_CHECK(in(rT3, rI)); + } +} + +template +void test_binary() { + typedef typename F::I I; + for(I a(-10., -9.91); a.lower() <= 10.; a += 0.3) { + for(I b(-10., -9.91); b.lower() <= 10.; b += 0.3) { + if (!F::validate(a, b)) continue; + T al = a.lower(), au = a.upper(), bl = b.lower(), bu = b.upper(); + I rII = F::f_II(a, b); + I rIT1 = F::f_IT(a, bl), rIT2 = F::f_IT(a, bu); + I rTI1 = F::f_TI(al, b), rTI2 = F::f_TI(au, b); + T rTT1 = F::f_TT(al, bl), rTT2 = F::f_TT(al, bu), + rTT3 = F::f_TT(au, bl), rTT4 = F::f_TT(au, bu); + BOOST_CHECK(in(rTT1, rIT1)); + BOOST_CHECK(in(rTT3, rIT1)); + BOOST_CHECK(in(rTT2, rIT2)); + BOOST_CHECK(in(rTT4, rIT2)); + BOOST_CHECK(in(rTT1, rTI1)); + BOOST_CHECK(in(rTT2, rTI1)); + BOOST_CHECK(in(rTT3, rTI2)); + BOOST_CHECK(in(rTT4, rTI2)); + BOOST_CHECK(subset(rIT1, rII)); + BOOST_CHECK(subset(rIT2, rII)); + BOOST_CHECK(subset(rTI1, rII)); + BOOST_CHECK(subset(rTI2, rII)); + } + } +} + +#define new_unary_bunch(name, op, val) \ + template \ + struct name { \ + typedef boost::numeric::interval I; \ + static I f_I(const I& a) { return op(a); } \ + static T f_T(const T& a) { return op(a); } \ + static bool validate(const I& a) { return val; } \ + } + +new_unary_bunch(bunch_pos, +, true); +new_unary_bunch(bunch_neg, -, true); +//new_unary_bunch(bunch_sqrt, sqrt, a.lower() >= 0.); +//new_unary_bunch(bunch_abs, abs, true); + +template +void test_all_unaries() { + BOOST_CHECKPOINT("pos"); test_unary >(); + BOOST_CHECKPOINT("neg"); test_unary >(); + // BOOST_CHECKPOINT("sqrt"); test_unary >(); + // BOOST_CHECKPOINT("abs"); test_unary >(); +} + +#define new_binary_bunch(name, op, val) \ + template \ + struct name { \ + typedef boost::numeric::interval I; \ + static I f_II(const I& a, const I& b) { return a op b; } \ + static I f_IT(const I& a, const T& b) { return a op b; } \ + static I f_TI(const T& a, const I& b) { return a op b; } \ + static T f_TT(const T& a, const T& b) { return a op b; } \ + static bool validate(const I& a, const I& b) { return val; } \ + } + +new_binary_bunch(bunch_add, +, true); +new_binary_bunch(bunch_sub, -, true); +new_binary_bunch(bunch_mul, *, true); +new_binary_bunch(bunch_div, /, !in_zero(b)); + +template +void test_all_binaries() { + BOOST_CHECKPOINT("add"); test_binary >(); + BOOST_CHECKPOINT("sub"); test_binary >(); + BOOST_CHECKPOINT("mul"); test_binary >(); + BOOST_CHECKPOINT("div"); test_binary >(); +} + +int test_main(int, char *[]) { + BOOST_CHECKPOINT("float tests"); + test_all_unaries (); + test_all_binaries (); + BOOST_CHECKPOINT("double tests"); + test_all_unaries(); + test_all_binaries(); + BOOST_CHECKPOINT("long double tests"); + test_all_unaries(); + test_all_binaries(); + return 0; +}