RE: [boost] Interval Wrinkles

From: Guillaume Melquiond (gmelquio_at_[hidden])
Date: 2003-04-13 10:49:37

On Sun, 13 Apr 2003, Paul A. Bristow wrote:

> > Is it only a typo in your mail or did you really try to convert the result
> > of 'abs(x)' to the type 'T'? The function 'abs' computes:
> >
> >   abs(X) == {abs(x) | x in X}
>
> Not being a mathematician, I naively expected abs of the interval [-1.9, -2.1]
> to return [1.9 , 2.1], but on giving it more a millisecond of thought, I can see
> that [-0.1, 0.1] would not be so useful. Worth documenting for other interval
> novices?

First, [-1.9,-2.1] is not a valid interval, [-2.1,-.19] is. :-)
And, abs([-2.1,-1.9]) == [1.9,2.1]; abs([-0.1,0.1]) == [0,0.1].

All these results are deduced from the inclusion property, which is
already in the documentation. Maybe some more examples of this property
should be added.

> > In the next version of the library, there will be a new function called
> > 'norm' that returns:
> >
> >   norm(X) = max {abs(x) | x in X}
> >
> > Maybe it's this kind of behavior you were expecting?
>
> I suspect so. Neal Becker's version 2 includes an externally defined norm
>
> template<typename T>
> T Norm (T x) { return x * x ; }

What's that? I hope it is not supposed to be what mathematicians call a
'norm'. Because it isn't (it doesn't respect the property: norm(a.x) ==
abs(a).norm(x) when a is a scalar). Maybe it is meant to be the square of
a norm? In this case, the naming should be changed in my opinion.

> If this exists, I'd be grateful for it.

I don't know what it is, so I cannot say if it exists. But if it's useful,
it can be added.

> > > 2    I found several places where expressions like T mean/integer
> > count failed
> > > when T was an interval type.  This is because there is no conversion from
> > > integer, though there is from double etc. writing mean/T(count) works OK.
> > >
> > > I half_x = x/2; // fails
> > > I half_x = x/2.0; // OK
> > > I half_x = x/ I(2); // OK
> > >
> > > This is a very common situation and would be helpful to avoid
> > widespread changes
> > > to existing code - but perhaps this is an exercise for the student :-)
> >
> > Yes, there are only mixed operations between expressions of type
> > 'interval<T,_>' and 'T'. So it is normal for the first expression to fail.
> > The correct way to do it is:
> >
> >    I half_x = x/T(2);
> >
> > We didn't include automatic promotion of 'int' to 'T' or 'interval<T,_>',
> > because the library don't have any information of the validity of the
> > conversion between 'int' and 'T'. If you think it is needed, we can add
> > it.
>
> It feels a reasonable assumption that an integer 2 is 'exact' and so should be
> converted to a 'singleton' interval [2.00000000000000000, 2.000000000000000000]?

Maybe I was not clear. What I wanted to say is that there is no promotion
mechanism in the library. If you manipulate data of type 'T' and
'interval<T>' (let's abbreviate it to 'I') and do the product T * I. The
data of type T is not converted to type I before doing the operation.
There is a function specialized in doing this particular product.

As a side note, because of the use of template operators, C++ compilers
would have not been able to do the conversion from T to I if only the
product I * I was implemented.

So the problem is not just the conversion from int to interval (the
compiler is already able to do the conversion int to interval<double>
since it will promote int to double when it encounters this line:
interval<double> a = 2;). The problem is the lack of mixed operators
between int and interval<T>.

For example, let's suppose you have this kind of code:

  void f(interval<double> a) {
    interval<double> b = 2 * a;
    ...

For this code to compile, there must exist a function whose signature
comes down after template instanciation to:

  interval<double> operator* (int, interval<double>);

And it is this kind of functions that doesn't exist in the library.

> Lots of type templated code is going to assume this (and so require lots of
> tiresome modifications to make it 'interval-type compatible'). So unless there
> is any really good 'interval maths disadvantage' I think this would be a good
> thing.

In order for the user to access integer functions, I can add this
kind of code:

  template<class T, class Policies>
  interval<T, Policies> operator* (int x, interval<T, Policies> y) {
    return static_cast<T>(x) * y;
  }

It will solve the issue. But it may induce some strange behavior when int
is not exactly convertible to T. For example, x is not always equal to
float(x) (for some big integers).

So, when designing the library, we were confronted to a problem. Should we
provide some functions that will improve the ease of use of the library
but could induce wrong results if the user is not cautious enough? It was
decided that the library should not trust the user... Maybe it's time for
a macro THE_USER_ASSUMES_ALL_THE_RESPONSABILITIES? Or to put these helper
functions in a separate header that would not be included by default.

> (In passing, I note the display style might optionally allow a display if the
> number is 'exact', for example [2.] rather than [2., 2.]?  I know you can get at
> this using the singleton function. I find the way the display of
> interval<double. changes with precision slightly surprising - as you change
> cout.precision(2) to precision(16) you get [1.9, 2.] ... [1.999999999999999,
> 2.].

The display of intervals is something we didn't want to do. It's up to the
user to define its own display. The header interval/io.hpp is merely
provided as a debug facilities (it isn't included by default). The main
reason is: << and >> usually don't handle correctly floating-point
numbers.

For example, at first, someone may think it's enough to do:

  out << '[' << x.lower() << ',' << x.upper() << ']';

However, if x is [2.1,2.2] and precision is 0, the output will be [2.,2.].
It isn't satisfying because the display doesn't enclose x (it should have
been [2.,3.]).

The correct way would be to print the exact representation of the bounds
(in a binary format for example). But there is no reason it's what the
user expects.

So you are free to do whatever you want for the display of intervals, it's
outside the scope of the library.

> > > 4   x1 = exp(x); // fails - exp_up and exp_down not a member ...

> > Yes, it's normal. Because we don't expect the mathematical functions to be
> > correctly rounded (and so the results would not be correct), the exp_up
> > and exp_down are not provided by default. In the first version of the
> > library, there were placeholders for these functions. But during the Boost
> > submission, the reviewers explained they prefered the compiler to return a
> > compilation error rather than the functions to return an unusable result
> > ('exp(x)' would have been '<0,+inf>'). So we removed them.
> >
> > > 5  similarly with log_up and log_down.
>
> If you are using an interval<double> x to hold information about uncertainty,
> then is it reasonable to expect exp(x) to give [exp(x.lower(), exp(x.higher()]?

It's a bit more complicated than that. When using interval<double>, the
bounds of your interval are double. However, the exponential of a double
precision floating-point number is not a double (exp is a transcendental
function). So it isn't enough to only compute [exp(lower),exp(upper)], the
function needs to compute [exp_down(lower),exp_up(upper)] in order for the
result to correctly contain all possible results.

However, and here lies the problem, there is no way to ask the
mathematical library (the one which provides the exp function) to compute
rounded down and rounded up versions of the exp function. Because of this
limitation, the user needs to rely on an external mathematical library
which provides them. And it's the reason why the library doesn't do
anything by default.

> So if the user wants this, he has to provide them?
> (I believe this will be a common requirement, and optional code examples would
> be helpful).

I will try to add an example which uses an external library.

> > > 6  Using interval pi I was puzzled that:
> > >
> > > 	using boost::numeric::interval_lib::pi<
> > boost::numeric::interval<double> >; //
> > > triggers compile fail at pi<I>()
> > > 	using boost::numeric::interval_lib::pi; // illegal use of expression
> > > 	I pi_d = pi<I>();
> > >
> > > whereas
> > >
> > > 	using namespace boost::numeric::interval_lib;
> > > 	I pi_d = pi<I>();
> > > and
> > >
> > > 	I pi_d = boost::numeric::interval_lib::pi<
> > boost::numeric::interval<double>
> > > >();
> > >
> > > are fine.
> >
> > I'm not sure to understand what you mean. Is it supposed to be another
> > deficiency of the compiler?
>
> Or is it my (lack of) understanding of C++?

This code:

  using boost::numeric::interval_lib::pi;
  I a = pi<I>();

works with GCC and Intel CC. And this code:

  using boost::numeric::interval_lib::pi<I>;
  I a = pi<I>();

is not allowed in C++.

> Thanks for your prompt assistance.
>
> Paul A Bristow, Prizet Farmhouse, Kendal, Cumbria, LA8 8AB  UK
> [phone]   Mobile [phone]
> Mobile mailto:pabristow_at_[hidden]
> mailto:pbristow_at_[hidden]

Guillaume