[boost] Re: Numeric Conversion Formal Review

From: Guillaume Melquiond (guillaume.melquiond_at_[hidden])
Date: 2003-12-12 06:24:49

Le ven 12/12/2003 à 00:37, Fernando Cacciola a écrit :
> Guillaume Melquiond wrote:

[...]

> > When doing float to int conversion, you use this formula: floor(f +
> > 0.5) Unfortunately, the formula won't work for some conversions. Here
> > is a
> > small snippet of a C program (the "volatile" keyword is used all over
> > the place so that you can test it on x86 and PPC).
> >
> > #include <math.h>
> > #include <stdio.h>
> > #include <assert.h>
> > int main() {
> >   int a = (1 << 24) - 1;
> >   volatile float b = a;
> >   volatile float c = b + 0.5f;
> >   volatile float d = floorf(c);
> >   int e = d;
> >   assert((double)a == (double)b);
> >   printf("%d %d\n", a, e);
> >   return 0;
> > }
> >
> > The integer is the biggest odd integer representable as simple
> > precision floating-point number. So the conversion from float to int
> > should not
> > modify the value (since it already is an integer). However, as you can
> > see, the final value is different from the first value. The same will
> > happen with any odd integer between 2^23 and 2^24 (and their negative
> > counterpart). So it's not a few cases, there is unfortunately a bunch
> > of them.
> >
> Hmmm...
> OK, I see the problem....
> Seems difficult to deal with without compromising efficiency for the safe
> cases (<2^23).
> Anyway, it is a problem rather specific of the fixed -known size- floating
> point formats (such as the builtin floats), so I figure that I'll only need
> to provide custom specialization of the rounder policies for these types.

Yes it is rather specific. It only happens when converting a
floating-point format to an integer format with a size bigger than the
size of the mantissa of the floating-point. So float to int (or bigger)
and double to (long) long int.

> > Another problem with rounding to nearest (that's because it's a
> > rounding
> > to "even"). This range may be too big: ( s >= S(LowestT)-S(0.5) ) &&
> > ( s <= S(HighestT)+S(0.5) ). It's probably better if it's only < for
> > the
> > positive part because the biggest positive number of an integer type
> > is usually odd.
> 
> Yes I see... essentially the same problem as before right? I think the
> change you suggest is appropriate, thanks!

No, it's a different problem. The previous problem was caused by the
fact the mantissa was not wide enough to handle the +0.5 without
changing the integer value. Now there is enough precision.

However, if the floating-point value is HighestT+0.5 with HighestT being
odd, the closest integer is HighestT+1 when rounding to nearest, ties to
even. Consequently you get an integer outside of your range. That's for
the theory (it doesn't depend on the implementation).

Now that I'm explaining the theory, I realize that there is a more
general problem when converting a floating-point number halfway between
two integers. Let's suppose k is an integer. If the floating-point value
is f=k+0.5, then the implementation always compute floor(f+0.5)=k+1 as
the rounded value. So the applied rounding is rounding to nearest, ties
to plus infinity.

So the problem is: this implemented rounding is not really a standard
one. The standard one is "ties to even" and in the future there may be
"ties away from zero" (in the revised ieee 754 that includes decimal
floating-point numbers).

So I'm not sure what the best thing to do is. At least, the
documentation should mention the default rounding to nearest is a
strange one. But it may be worthwhile to use a more thorough
implementation.

And rather than only relying only on floor and ceil, maybe a third
function could be used (it's just a quick thought, there probably are
some drawbacks).

> > Concerning the documentation, there are a few places where the < and >
> > are badly written. By rgrep-ing with '[^&][lg]t;' and '&[lg]t[^;]' you
> > should find all the occurrences.
> 
> Oh, Ok. I'll fix them.
> 
> >
> > In the definition part of the documentation, the notions of "range"
> > of a
> > type and "overflow" don't match the later explanation of float to int
> > rounding.
> I can't find what 'later explanation' are you referring to... or I don't see
> the error :-)
> Can you be more specific?

In the definition part, you explain a value is out of range if it is
smaller than the smallest target value or bigger than the biggest target
value. But, in the rounding to integer part, a value is out of range if
it is outside of the destination range plus a fuzz (-1/0, -0.5/+0.5,
0/+1). Maybe I'm misunderstanding the definitions.

> Anyway, this part of the documentation is intended to be as consistent and
> precise as possible, so that is a good place to be pedantic :-)
> The more and finer corrections I get the better the result.
> 
> > You should also be a bit more cautious when speaking of
> > "correct rounding". Sometimes what you describe is only a "faithful
> > rounding" (a faithful rounding becomes correct only when the direction
> > of rounding is precisely known). But I'm splitting hairs.
> 
> I never quite understood the difference :-)
> So if you can help me out with the right descriptions I'd much appreciate
> it.

The rounding is faithful when the user can't foretell the direction of
rounding. When the result is only guaranteed to be less than 1 ulp away
to the real value, the rounding is faithful.

The rounding is directed when a direction is given. Rounding towards
-inf, towards +inf, towards zero, away from zero, etc are directed
rounding.

And finally there is the rounding to nearest. There are lot of them
depending on the way the ties are handled (random, to even, to odd,
toward zero, away from zero, etc).

Depending on the people, correct rounding may mean two things: directed
rounding and to nearest rounding, or only to nearest rounding.

In fact, the only problem may be this sentence: "if v==prev or v==next,
the representation is correctly rounded". You should remove "correctly"
and put the "correctly rounded" of the next sentence in bold. Later you
can add that in case of rounding to nearest, the precision is 1/2 ulp
and not just 1 ulp.

Finally in the "Standard (numeric) Conversions", you may want to use
"faithfully" rather than "correctly" since the user can't know the
direction of rounding unless she takes a look at the implementation
(both software and hardware).

> > Finally what is the rational for returning 0 in bounds::smallest when
> > it's an integer type?
> 
> The decision was rather arbitrary (way too arbitrary :-)
> Now after thinking about it I realize that it should be 1, not 0; which is
> the proper parallel to the smallest float value (which happens not to be 0)

Yes, I think so.

> > And why return the smallest normalized
> > floating-point number? Why not the smallest floating-point number?
> 
> This is because the ISO98 standard doesn't know about subnormals, so the
> smallest value that can be portably represented as far as C++ is concerned
> must be normalized.
> (this is why numeric_limits<a_float_type>::min() returns a normalized
> number)
> 
> > Or even better the two of them?
> 
> Because subnormals don't legally exist in C++.

Yes, and they also are not forbidden, are they? It's why I suggested to
put both a smallest and a really_smallest (hum...) values.

> > And what about ulp(1) or the smallest
> > floating-point number following 1?
> >
> This is way too big, isn't it?
> smallest() substitutes numeric_limits<a_float_type>::min().

No, I was just suggesting adding a new field, not replacing smallest. I
rather had in mind the necessary values in order to fully define a
classical floating-point type.

> Thanks for your mini-review....
> I'm looking forward for a full review from you!
> 
> Best regards,
> 
> Fernando Cacciola
> SciSoft

Regards,

Guillaume