Topic: Floating point arithmetic


Author: Greg Herlihy <greghe@mac.com>
Date: Sat, 21 Mar 2009 12:53:22 CST
Raw View
On Mar 5, 7:06    am, Kai-Uwe Bux <jkherci...@gmx.net> wrote:
>
> Clause [5/10] allows excess precision for expressions involving floating
> point numbers. The question is, at which points the standard _requires_
> truncation of that excess precision. Put:
>
>          double x = ...
>          double y = ...

Essentially, an implementation may take advantage of extended
precision floating point calculations but only as long as the
semantics of the program itself are not affected. So, although interim
floating point calculations might use extended precisions, assignment
or conversions ("load and store" operations) of floating point values
to floating point types - may not.

>          bool equal_a ( double lhs, double rhs ) {
>               return ( lhs == rhs );
>          }
>
>          bool equal_b ( double const & lhs, double const & rhs ) {
>               return ( lhs == rhs );
>          }
>
> With regard to initializing parameters, consider:
>
> a) The expression x*y == x*y seems to be allwed to yield false since either
> side is allowed but neither side is required to be evaluated with excess
> precision and the built-in comparison that is part of the expression can
> compare _different_ floating point numbers.

No, multiplying x and y yields a single value - not a range of values.
And although implementations may differ on how best to represent the
product of x and y - they must all agree that that product has in fact
only one value. Or, as the C99 Standard explains:

"Implementations employing wide registers have to take care to honor
appropriate semantics. Values are independent of whether they are
represented in a register or in memory. For example, an implicit
"spilling" of a register is not permitted to alter the value. Also,an
explicit "store and load" is required to round to the precision of the
storage type. In particular, casts and assignments are required to
perform their speci   ed conversion." [C99 5.1.2.3]

In other words, if a program calculates "x*y" with extended precision
- then the program must always calculate x*y with extended precision -
because multiplication is a deterministic operation (even if the
precision of the result is left up to the implementation.). Therefore
x*y must yield a single, deterministic value for any given
implementation. C++'s program execution model, moreover, requires that
an implementation produce identical output when executing identical
programs that perform the same, deterministic operations on identical
input. So an implementation that allowed the product of two
identically-valued, identical types to vary - would not be able to
adhere to C++'s program execution model.

> b) The expression equal_b( x*y, x*y ) seems to be required to yield true:
> the initialization of the double objects in passing the parameters involves
> truncating the excess precision. But I am not sure.

Both "x*y" will yield identical values - whether those identical
values compare equal depends on the value itself (in particular, the
equality test will fail if x*y is a NaN, otherwise it will succeed).

> c) About the expression equal_a( x*y, x*y ), I don't know. As opposed to
> equal_b, the function equal_a takes its parameters by value. Therefore,
> initializing the parameters might not involve the construction of a
> floating point object.
>
> d) is std::sin(x) == std::sin(x) required to yield true by the standard?

A sine function (like most mathematical functions) establishes a one-
to-one relation between the value of its argument and the value of its
result. Therefore, calling std::sin() twice, with identical arguments
on the right and left side of an equality comparison - will yield
identical values on both sides. Unless (and a robust implementation is
expected to handle this possibility) the value of pi were to change
between the first and second call to std::sin().

Greg

--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: "Joe Smith" <unknown_kev_cat@hotmail.com>
Date: Sun, 22 Mar 2009 21:27:58 CST
Raw View
"James Kanze" wrote:
>
>     A floating expression may be contracted, that is,
>     evaluated as though it were an atomic operation, thereby
>     omitting rounding errors implied by the source code and
>     the expression evaluation method. The FP_CONTRACT pragma
>     in <math.h> provides a way to disallow contracted
>     expressions. Otherwise, whether and how expressions are
>     contracted is implementation-defined.
>
> In other words, others, in addition to myself, felt the need for
> this extension of the precision to be predictable (for a given
> compiler), and changed the wording to clearly require it.  (I'm
> not sure I like the new wording that much, but it clearly makes
> stricter requirements than the wording in the C++ standard.)
>
You do realize that (in practise at least) FP_CONTRACT is really about
allowing the use of special multi-step instructions, such as a
floating-multiply-and-add instruction which may have a smaller rounding
error than a multiply followed by an add. (Compilers might be allowed to do
some other optimizations when it is active, but I've not heard of any that
do.)

I have seen nothing that indicated it has a bearing on the extended
precision issue.

I do agree that it is not good for C++ and C to be out-of-sync with regards
to extended precision.

I would ideally prefer is the standards required trucantion to double or
single (as per the variable datatype) for any comparison operators
(including the equality operator), but at the least they really should be
conistant with each other.


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Kai-Uwe Bux <jkherciueh@gmx.net>
Date: Thu, 5 Mar 2009 09:06:36 CST
Raw View
Hi,


over in comp.lang.c++

http://groups.google.com/group/comp.lang.c++/browse_frm/thread/67d1aa830ced63bb#

we ran into a problem that requires a language lawer.


Clause [5/10] allows excess precision for expressions involving floating
point numbers. The question is, at which points the standard _requires_
truncation of that excess precision. Put:

   double x = ...
   double y = ...

   bool equal_a ( double lhs, double rhs ) {
     return ( lhs == rhs );
   }

   bool equal_b ( double const & lhs, double const & rhs ) {
     return ( lhs == rhs );
   }

With regard to initializing parameters, consider:

a) The expression x*y == x*y seems to be allwed to yield false since either
side is allowed but neither side is required to be evaluated with excess
precision and the built-in comparison that is part of the expression can
compare _different_ floating point numbers.

b) The expression equal_b( x*y, x*y ) seems to be required to yield true:
the initialization of the double objects in passing the parameters involves
truncating the excess precision. But I am not sure.

c) About the expression equal_a( x*y, x*y ), I don't know. As opposed to
equal_b, the function equal_a takes its parameters by value. Therefore,
initializing the parameters might not involve the construction of a
floating point object.


With regard to return statements, one could argue that each return involves
copy-initialzation and hence the result of a function is not allowed excess
precision. On the other hand, RVO in conjunction with [5/10] may give
license to bypass copy-initialization. So:

d) is std::sin(x) == std::sin(x) required to yield true by the standard?


Best

Kai-Uwe Bux

--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: "Alf P. Steinbach" <alfps@start.no>
Date: Sun, 8 Mar 2009 10:13:47 CST
Raw View
* Kai-Uwe Bux:
>
> over in comp.lang.c++
>
> http://groups.google.com/group/comp.lang.c++/browse_frm/thread/67d1aa830ced63bb#
>
> we ran into a problem that requires a language lawer.

[snipped examples of how fp comparisions *may* yield arbitrary results]


The larger issue is IMHO the non-well-definedness of C++ operations.

Someone thinking the standard can't really allow arbitrary fp comparision
results, or other such apparently insane behavior, may just code away,
happily
unaware of the quicksand underlying the apparently firm surface of C++
operations. After testing they have an app that works on at least the
tested
platform(s) and with the relevant compiler(s) and option(s). Hurray.

Someone else, taking a language lawyer's word that /this/ is what the
standard
guarantees and /that/ is what isn't guaranteed, well, with this kind of
subtlety
the code to avoid the missing guarantees will introduce complexity and
inefficiency, and the code exploiting the guarantees will fail to
achieve the
indended result with actual compilers and environments. This is
guaranteed by
Murphy, and not the least by the obfuscation of the issues that's
effectively
present in the standard due to the IMHO wrong idea of "UB allows
efficiency".

In order to be a useful programming language standard, there should be
nothing
to discuss about the effect of elementary operations. For when the
guarantees or
absence of guarantees for fundamental operations can be discussed, then the
effect of any program is very much in question. Which then means that
it's only
in name and "officially" that the document is a standard.

So I think it's really high time that all that fundamental op UB stuff and
exquisite multi-layered indirect subtlety is replaced with less
interesting but
far more practically useful well-definedness and simplity and obviousness.


Cheers,

- Alf

--
Due to hosting requirements I need visits to [http://alfps.izfree.com/].
No ads, and there is some C++ stuff! :-) Just going there is good. Linking
to it is even better! Thanks in advance!

[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: James Kanze <james.kanze@gmail.com>
Date: Sun, 8 Mar 2009 17:15:59 CST
Raw View
On Mar 8, 5:13 pm, "Alf P. Steinbach" <al...@start.no> wrote:
> * Kai-Uwe Bux:
> > over in comp.lang.c++

> >http://groups.google.com/group/comp.lang.c++/browse_frm/thread/67d1aa...

> > we ran into a problem that requires a language lawer.

> [snipped examples of how fp comparisions *may* yield arbitrary results]

> The larger issue is IMHO the non-well-definedness of C++
> operations.

> Someone thinking the standard can't really allow arbitrary fp
> comparision results, or other such apparently insane behavior,
> may just code away, happily unaware of the quicksand
> underlying the apparently firm surface of C++ operations.
> After testing they have an app that works on at least the
> tested platform(s) and with the relevant compiler(s) and
> option(s).

For the values they actually tested.  Testing doesn't guarantee
correctness, especially where floating point is concerned.

> Someone else, taking a language lawyer's word that /this/ is
> what the standard guarantees and /that/ is what isn't
> guaranteed, well, with this kind of subtlety the code to avoid
> the missing guarantees will introduce complexity and
> inefficiency, and the code exploiting the guarantees will fail
> to achieve the indended result with actual compilers and
> environments. This is guaranteed by Murphy, and not the least
> by the obfuscation of the issues that's effectively present in
> the standard due to the IMHO wrong idea of "UB allows
> efficiency".

It's not a question of just language lawyering.  You have to
know what you can and cannot count on.

> In order to be a useful programming language standard, there
> should be nothing to discuss about the effect of elementary
> operations. For when the guarantees or absence of guarantees
> for fundamental operations can be discussed, then the effect
> of any program is very much in question. Which then means that
> it's only in name and "officially" that the document is a
> standard.

> So I think it's really high time that all that fundamental op
> UB stuff and exquisite multi-layered indirect subtlety is
> replaced with less interesting but far more practically useful
> well-definedness and simplity and obviousness.

And how does that affect floating point operations?  Java
limited itself to IEEE (which C++ intentionally does not), and
attempted to define the language so that the results were
guaranteed to be the same, everywhere.  In the end, they had to
back off, because the runtime cost was too high on some
platforms.

--
James Kanze (GABI Software)             email:james.kanze@gmail.com
Conseils en informatique orient   e objet/
                    Beratung in objektorientierter Datenverarbeitung
9 place S   mard, 78210 St.-Cyr-l'   cole, France, +33 (0)1 30 23 00 34


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: "Alf P. Steinbach" <alfps@start.no>
Date: Mon, 9 Mar 2009 02:43:40 CST
Raw View
* James Kanze:
>
> On Mar 8, 5:13 pm, "Alf P. Steinbach" <al...@start.no> wrote:
>
>> So I think it's really high time that all that fundamental op
>> UB stuff and exquisite multi-layered indirect subtlety is
>> replaced with less interesting but far more practically useful
>> well-definedness and simplity and obviousness.
>
> And how does that affect floating point operations?

Among other things, the issues Kai-Uwe raised about comparisions. Fair
enough that we don't know what the result of 1.2+1.3 is. But we should
know, as an example, that 1.2+1.3 == 1.2+1.3  --  without consulting a
language lawyer.


>  Java
> limited itself to IEEE (which C++ intentionally does not), and
> attempted to define the language so that the results were
> guaranteed to be the same, everywhere.  In the end, they had to
> back off, because the runtime cost was too high on some
> platforms.

Good point.

Showing one way not to do it.

But I'm not sure it's the complete story, because the (previous) IEEE
standard wasn't and isn't all that confining, except possibly for the
issue of NaN comparisions (which as I understand it is where g++ had
to "back down", via an option ;-)). Perhaps it shows that one needs at
least at two floating point types, safe versus fast, with implicit
conversion only one way. But no matter solution, which should better
be suggested by someone with very solid experience (in various
environments) in numerical work, having solid foundation is key.


Cheers,

- Alf

--
Due to hosting requirements I need visits to <url: http://alfps.izfree.com/>.
No ads, and there is some C++ stuff! :-) Just going there is good. Linking
to it is even better! Thanks in advance!

[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: James Kanze <james.kanze@gmail.com>
Date: Mon, 9 Mar 2009 11:13:05 CST
Raw View
On Mar 9, 9:43 am, "Alf P. Steinbach" <al...@start.no> wrote:
> * James Kanze:
> > On Mar 8, 5:13 pm, "Alf P. Steinbach" <al...@start.no> wrote:

> >> So I think it's really high time that all that fundamental
> >> op UB stuff and exquisite multi-layered indirect subtlety
> >> is replaced with less interesting but far more practically
> >> useful well-definedness and simplity and obviousness.

> > And how does that affect floating point operations?

> Among other things, the issues Kai-Uwe raised about
> comparisions. Fair enough that we don't know what the result
> of 1.2+1.3 is. But we should know, as an example, that 1.2+1.3
> == 1.2+1.3  --  without consulting a language lawyer.

Na   vely, I'd tend to agree with you.  Practically, Kahan and
others introduced extended precision for a reason, and they're a
lot better expert in numerics than I am, so I'll assume that
there are legitimate arguments in favor of using extended
precision for the intermediate results.  I know as well that
some numeric specialists disagree.  But I'm not really competent
to judge the arguments, one way or another.

> > Java limited itself to IEEE (which C++ intentionally does
> > not), and attempted to define the language so that the
> > results were guaranteed to be the same, everywhere.  In the
> > end, they had to back off, because the runtime cost was too
> > high on some platforms.

> Good point.

> Showing one way not to do it.

> But I'm not sure it's the complete story, because the
> (previous) IEEE standard wasn't and isn't all that confining,
> except possibly for the issue of NaN comparisions (which as I
> understand it is where g++ had to "back down", via an option
> ;-)). Perhaps it shows that one needs at least at two floating
> point types, safe versus fast, with implicit conversion only
> one way. But no matter solution, which should better be
> suggested by someone with very solid experience (in various
> environments) in numerical work, having solid foundation is
> key.

At the very least, compilers on platforms where it is an issue
need to provide options to control it.  On an Intel, for
example, I can imagine three variants:

  -- every intermediate result is rounded to double---this is
     what Java originally required---;

  -- all casts, assignments and initializations forced to double,
     but intermediate values use extended precision---this is
     conformant to the standard, but from a QoI point of view, a
     compiler which does this should do it systematically,
     including splills to memory (and part of the question here
     is whether return values and function arguments count as
     "initializations"); and

  -- the compiler goes all out for speed, even when this leads to
     counter-intuitive or even non-conformant behavior (caveat
     emptor).

Personally, I've favor the first option being the default, but I
can understand arguments for the second, *provided* the extended
precision is predictable.  (Thus, for example, "1.2+1.3 ===
1.2+.1.3" is guaranteed, because everything is extended
precision.  "double x = 1.2 + 1.3 ; if ( x == 1.2 + 1.3 )..."
isn't.)

The issue here is, I think, really where the standard allows
such extended precision; there are cases where it isn't really
very clear.  I would also argue that in such cases, the standard
should require compilers to be coherent---they shouldn't be
allowed to change the precision just because a value spills from
register.  (All MHO, of course.)

--
James Kanze (GABI Software)             email:james.kanze@gmail.com
Conseils en informatique orient   e objet/
                    Beratung in objektorientierter Datenverarbeitung
9 place S   mard, 78210 St.-Cyr-l'   cole, France, +33 (0)1 30 23 00 34


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Jean-Marc Bourguet <jm@bourguet.org>
Date: Mon, 9 Mar 2009 11:13:45 CST
Raw View
"Alf P. Steinbach" <alfps@start.no> writes:

> But I'm not sure it's the complete story, because the (previous) IEEE
> standard wasn't and isn't all that confining, except possibly for the
> issue of NaN comparisions (which as I understand it is where g++ had
> to "back down", via an option ;-)).

I don't think there is much problem with NaN.

I know there are performance problem with subnormals (this seems the new
favourite denonimation of what was called previously denormal) whose use
trap and backup to software handling in some processor.  (ISTR that some
processors needs software handling even in their absence if one want
subnormals to be handled correctly).

Another area which could be problematic is the one of rounding (IEEE 754
has four rounding modes)

Yours,

--
Jean-Marc

[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Bart van Ingen Schenau <bart@ingen.ddns.info>
Date: Tue, 10 Mar 2009 18:19:41 CST
Raw View
James Kanze wrote:

> On Mar 9, 9:43 am, "Alf P. Steinbach" <al...@start.no> wrote:
>> * James Kanze:
>> > On Mar 8, 5:13 pm, "Alf P. Steinbach" <al...@start.no> wrote:
>
>> >> So I think it's really high time that all that fundamental
>> >> op UB stuff and exquisite multi-layered indirect subtlety
>> >> is replaced with less interesting but far more practically
>> >> useful well-definedness and simplity and obviousness.
>
>> > And how does that affect floating point operations?
>
>> Among other things, the issues Kai-Uwe raised about
>> comparisions. Fair enough that we don't know what the result
>> of 1.2+1.3 is. But we should know, as an example, that 1.2+1.3
>> == 1.2+1.3  --  without consulting a language lawyer.
>
> Na   vely, I'd tend to agree with you.  Practically, Kahan and
> others introduced extended precision for a reason, and they're a
> lot better expert in numerics than I am, so I'll assume that
> there are legitimate arguments in favor of using extended
> precision for the intermediate results.  I know as well that
> some numeric specialists disagree.  But I'm not really competent
> to judge the arguments, one way or another.

I will add my agreement here.
As a layman, my impression is that what trips people up most, with
regard to the extended precision of floating point intermediaries, is
that two expressions that should calculate the same result do not
compare equal, even if the inherent limitation of floating point types
are not an issue for the expressions in question.
IOW, it confuses the hell out of people that the equation X == X may
yield false for an arbitrary expression X, even if it is known that the
result of X is not a NaN.

What would break if we limit the use of extended precision a bit, by
disallowing it for the (direct) operands of the equality and relational
operators?

<snip>
> At the very least, compilers on platforms where it is an issue
> need to provide options to control it.  On an Intel, for
> example, I can imagine three variants:
>
>   -- every intermediate result is rounded to double---this is
>      what Java originally required---;
>
>   -- all casts, assignments and initializations forced to double,
>      but intermediate values use extended precision---this is
>      conformant to the standard, but from a QoI point of view, a
>      compiler which does this should do it systematically,
>      including splills to memory (and part of the question here
>      is whether return values and function arguments count as
>      "initializations"); and

I would like to add option 2a:
   -- all casts, assignments, initializations and operands of relational
      and equality operators forced to double, but (other) intermediate
      values use extended precision; and

>
>   -- the compiler goes all out for speed, even when this leads to
>      counter-intuitive or even non-conformant behavior (caveat
>      emptor).
>
> Personally, I've favor the first option being the default, but I
> can understand arguments for the second, *provided* the extended
> precision is predictable.  (Thus, for example, "1.2+1.3 ===
> 1.2+.1.3" is guaranteed, because everything is extended
> precision.  "double x = 1.2 + 1.3 ; if ( x == 1.2 + 1.3 )..."
> isn't.)
>
> The issue here is, I think, really where the standard allows
> such extended precision; there are cases where it isn't really
> very clear.  I would also argue that in such cases, the standard
> should require compilers to be coherent---they shouldn't be
> allowed to change the precision just because a value spills from
> register.  (All MHO, of course.)
>
> --
> James Kanze (GABI Software)             email:james.kanze@gmail.com
>

Bart v Ingen Schenau
--
a.c.l.l.c-c++ FAQ: http://www.comeaucomputing.com/learn/faq
c.l.c FAQ: http://c-faq.com/
c.l.c++ FAQ: http://www.parashift.com/c++-faq-lite/

[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Kai-Uwe Bux <jkherciueh@gmx.net>
Date: Wed, 11 Mar 2009 16:49:42 CST
Raw View
Bart van Ingen Schenau wrote:

> James Kanze wrote:
>
>> On Mar 9, 9:43 am, "Alf P. Steinbach" <al...@start.no> wrote:
>>> * James Kanze:
>>> > On Mar 8, 5:13 pm, "Alf P. Steinbach" <al...@start.no> wrote:
>>
>>> >> So I think it's really high time that all that fundamental
>>> >> op UB stuff and exquisite multi-layered indirect subtlety
>>> >> is replaced with less interesting but far more practically
>>> >> useful well-definedness and simplity and obviousness.
>>
>>> > And how does that affect floating point operations?
>>
>>> Among other things, the issues Kai-Uwe raised about
>>> comparisions. Fair enough that we don't know what the result
>>> of 1.2+1.3 is. But we should know, as an example, that 1.2+1.3
>>> == 1.2+1.3  --  without consulting a language lawyer.
>>
>> Na   vely, I'd tend to agree with you.  Practically, Kahan and
>> others introduced extended precision for a reason, and they're a
>> lot better expert in numerics than I am, so I'll assume that
>> there are legitimate arguments in favor of using extended
>> precision for the intermediate results.  I know as well that
>> some numeric specialists disagree.  But I'm not really competent
>> to judge the arguments, one way or another.
>
> I will add my agreement here.
> As a layman, my impression is that what trips people up most, with
> regard to the extended precision of floating point intermediaries, is
> that two expressions that should calculate the same result do not
> compare equal, even if the inherent limitation of floating point types
> are not an issue for the expressions in question.
> IOW, it confuses the hell out of people that the equation X == X may
> yield false for an arbitrary expression X, even if it is known that the
> result of X is not a NaN.
>
> What would break if we limit the use of extended precision a bit, by
> disallowing it for the (direct) operands of the equality and relational
> operators?

I would like to stress (not add, since it's already in your question) that
it is important to require truncation for _all_ relational operators.
The discussion on comp.lang.c++, which triggered these questions was about a
comparison predicate that crashed sort. Consider something like:

 struct Point {
   double x;
   double y;
 };

 bool l1_is_shorter ( Point const & lhs, Point const & rhs ) {
   return ( abs( lhs.x ) + abs( lhs.y ) < abs( rhs.x ) + abs( rhs.y ) );
 }

Such a comparison predicate might yield inconsistent answers and fail to be
a strict weak order due to floating point issues.


Best

Kai-Uwe Bux
[snip]


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Geert-Jan Giezeman <geert@cs.uu.nl>
Date: Wed, 11 Mar 2009 16:49:42 CST
Raw View
Bart van Ingen Schenau wrote:

> I will add my agreement here.
> As a layman, my impression is that what trips people up most, with
> regard to the extended precision of floating point intermediaries, is
> that two expressions that should calculate the same result do not
> compare equal, even if the inherent limitation of floating point types
> are not an issue for the expressions in question.
> IOW, it confuses the hell out of people that the equation X == X may
> yield false for an arbitrary expression X, even if it is known that the
> result of X is not a NaN.
>
> What would break if we limit the use of extended precision a bit, by
> disallowing it for the (direct) operands of the equality and relational
> operators?
>

>
> I would like to add option 2a:
>   -- all casts, assignments, initializations and operands of relational
>      and equality operators forced to double, but (other) intermediate
>      values use extended precision; and


I agree that this is most annoying. Let me give a use case that
occurred in my code recently. It is a minimisation algorithm that uses
an evolutionary algorithm. I want to stop when no progress is made for
some number of generations. Here is the relevant part of the code:

struct evolutionary_pool
{
   virtual void next_generation()=0;
   virtual double current_minimum() const=0;
};

void optimise(evolutionary_pool &pool)
{
   const int Max = 10000;
   double best_till_now = pool.current_minimum();
   int no_progress = 0;
   while (no_progress < Max) {
       pool.next_generation();
       if (pool.current_minimum() < best_till_now) { // Not safe now
           best_till_now = pool.current_minimum();
           no_progress = 0;
       } else {
           ++no_progress;
       }
   }
}


Because in this case there are a finite number of internal states and
each state has a fixed value, I thought that this program is
guaranteed to terminate. After the recent discussions about floating
point issues, I now realise that this is not the case. Even if the
internal state doesn't change, and the value returned by
'current_minimum' stays the same, the comparison marked 'Not safe now'
may yield true each time.
Highly surprising to me.


Geert-Jan Giezeman

--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: James Kanze <james.kanze@gmail.com>
Date: Wed, 11 Mar 2009 23:27:52 CST
Raw View
On Mar 11, 1:19 am, Bart van Ingen Schenau <b...@ingen.ddns.info>
wrote:
> James Kanze wrote:
> > On Mar 9, 9:43 am, "Alf P. Steinbach" <al...@start.no> wrote:

      [...]
> As a layman, my impression is that what trips people up most,
> with regard to the extended precision of floating point
> intermediaries, is that two expressions that should calculate
> the same result do not compare equal, even if the inherent
> limitation of floating point types are not an issue for the
> expressions in question.  IOW, it confuses the hell out of
> people that the equation X == X may yield false for an
> arbitrary expression X, even if it is known that the result of
> X is not a NaN.

> What would break if we limit the use of extended precision a
> bit, by disallowing it for the (direct) operands of the
> equality and relational operators?

I don't think it would change much.  You'd still have the
problem that something like a*b+c == a*b+c might evaluate false.

What would help, IMHO, is requiring it to be an all or nothing
thing.  The compiler can use extended precision where ever it
wants, but if it does, it must do so consistently; it can't
drop back to the standard precision just because it has to spill
some intermediate values to memory.  (This could cause problems
for a processor which uses extended precision in its FPU, but
has no instructions for reading or writing values with this
precision.  I don't know of such a processor, however.)

> <snip>
> > At the very least, compilers on platforms where it is an
> > issue need to provide options to control it.  On an Intel,
> > for example, I can imagine three variants:

> >   -- every intermediate result is rounded to double---this is
> >      what Java originally required---;

> >   -- all casts, assignments and initializations forced to double,
> >      but intermediate values use extended precision---this is
> >      conformant to the standard, but from a QoI point of view, a
> >      compiler which does this should do it systematically,
> >      including splills to memory (and part of the question here
> >      is whether return values and function arguments count as
> >      "initializations"); and

> I would like to add option 2a:
>    -- all casts, assignments, initializations and operands of relational
>       and equality operators forced to double, but (other) intermediate
>       values use extended precision; and

I think this is an added comlexity, for no real benefit.  The
rule would basically be that a compiler is allowed to use
extended precision for "temporaries", but that it is required to
1) document this (it should be implementation defined, not
unspecified), and 2) use the same precision for all
"temporaries".

Given that, it remains for the standard to define exactly what
is meant by "temporaries" here---are function return values
"temporaries", in this sense, or not, for example.

As the standard currently reads, they are temporaries once they
have been constructed, but the rules for initializing them seem
to imply that they behave as concrete objects during
initialization.  The simplest rule, IMHO, would be that if it
has a name, it isn't a temporary, and if it doesn't, it is.  So
return values are temporaries (and so would have extended
precision), but function parameters aren't.  (I don't really
like the idea of return values and parameters being treated
differently, but I suspect you can't have everything.)  What is
thrown in an exception is also a temporary, according to this,
but it doesn't matter much, because the temporary has to be used
to initialize an object before it is used (and I doubt the
behavior of floating point types used as exceptions really
affects very many programs anyway).

--
James Kanze (GABI Software)             email:james.kanze@gmail.com
Conseils en informatique orient   e objet/
                    Beratung in objektorientierter Datenverarbeitung
9 place S   mard, 78210 St.-Cyr-l'   cole, France, +33 (0)1 30 23 00 34


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Greg Herlihy <greghe@mac.com>
Date: Thu, 12 Mar 2009 10:59:24 CST
Raw View
On Mar 5, 7:06    am, Kai-Uwe Bux <jkherci...@gmx.net> wrote:
>
> Clause [5/10] allows excess precision for expressions involving floating
> point numbers. The question is, at which points the standard _requires_
> truncation of that excess precision. Put:
>
>          double x = ...
>          double y = ...

Essentially, an implementation may take advantage of extended
precision floating point calculations - but only as long as the
semantics of the program itself are not affected. For example, if a
floating point expression like (a*b) == (c*d) is true when calculated
without extended precision, then (a*b) == (c*d) must still be true
when the identical multiplication is performed with extended
precision. A particular value is - in a way - independent of the
precision of its representation.

>          bool equal_a ( double lhs, double rhs ) {
>               return ( lhs == rhs );
>          }
>
>          bool equal_b ( double const & lhs, double const & rhs ) {
>               return ( lhs == rhs );
>          }
>
> With regard to initializing parameters, consider:
>
> a) The expression x*y == x*y seems to be allwed to yield false since either
> side is allowed but neither side is required to be evaluated with excess
> precision and the built-in comparison that is part of the expression can
> compare _different_ floating point numbers.

No, the x*y comparison must be made between identical values. The
explanation is simple: multiplying x and y yields a single value - not
a range of values. And although implementations may differ on how
accurately to represent the product of x and y - they must all agree
that that product has in fact only one value (that is, that the
equation has only one solution). As the C99 Standard explains:

"Implementations employing wide registers have to take care to honor
appropriate semantics. Values are independent of whether they are
represented in a register or in memory. For example, an implicit
"spilling" of a register is not permitted to alter the value. Also,an
explicit "store and load" is required to round to the precision of the
storage type. In particular, casts and assignments are required to
perform their speci   ed conversion." [C99   5.1.2.3]

In other words, if a program calculates "x*y" with extended precision
- then the program must always calculate x*y with extended precision
(at least in those situations where the extended precision would make
a difference) - because multiplication is a deterministic operation
that yields a single value (even if the precision of its
representation is unspecified.)

C++'s program execution model, moreover, requires that an
implementation produce identical output when executing identical
programs that perform defined operations on identical input. An
implementation that allowed the product of two identically-valued,
identical types to differ - would therefore not be able to adhere to C+
+'s program execution model, because such an implementation could not
guarantee consistent results as required.

> b) The expression equal_b( x*y, x*y ) seems to be required to yield true:
> the initialization of the double objects in passing the parameters involves
> truncating the excess precision. But I am not sure.

The "x*y" expression will yield identical values in all three cases.
Now, whether those identical values compare equal - depends on the
value itself (Specifically, the equality test will fail if x*y is a
NaN, otherwise it will succeed).

> c) About the expression equal_a( x*y, x*y ), I don't know. As opposed to
> equal_b, the function equal_a takes its parameters by value. Therefore,
> initializing the parameters might not involve the construction of a
> floating point object.
>
> d) is std::sin(x) == std::sin(x) required to yield true by the standard?

A sine function (like most mathematical functions) establishes a one-
to-one relation between the value of its input and the value of its
output (or result). Therefore, calling std::sin() with identical
arguments on the right and left side of an equality comparison - will
yield identical values on both sides. Unless (and a robust
implementation will handle this possibility) the value of pi were to
change between the first and second call to std::sin().

Greg



--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: Kai-Uwe Bux <jkherciueh@gmx.net>
Date: Thu, 12 Mar 2009 16:59:15 CST
Raw View
Greg Herlihy wrote:

> On Mar 5, 7:06?? am, Kai-Uwe Bux <jkherci...@gmx.net> wrote:
>>
>> Clause [5/10] allows excess precision for expressions involving floating
>> point numbers. The question is, at which points the standard _requires_
>> truncation of that excess precision. Put:
>>
>> ??  ?? double x = ...
>> ??  ?? double y = ...
>
> Essentially, an implementation may take advantage of extended
> precision floating point calculations - but only as long as the
> semantics of the program itself are not affected. For example, if a
> floating point expression like (a*b) == (c*d) is true when calculated
> without extended precision, then (a*b) == (c*d) must still be true
> when the identical multiplication is performed with extended
> precision. A particular value is - in a way - independent of the
> precision of its representation.


I don't see that this is supported by the standard. [5/10] seems to say
otherwise.


>> ??  ?? bool equal_a ( double lhs, double rhs ) {
>> ??  ??  ?? return ( lhs == rhs );
>> ??  ?? }
>>
>> ??  ?? bool equal_b ( double const & lhs, double const & rhs ) {
>> ??  ??  ?? return ( lhs == rhs );
>> ??  ?? }
>>
>> With regard to initializing parameters, consider:
>>
>> a) The expression x*y == x*y seems to be allwed to yield false since
>> either side is allowed but neither side is required to be evaluated with
>> excess precision and the built-in comparison that is part of the
>> expression can compare _different_ floating point numbers.
>
> No, the x*y comparison must be made between identical values. The
> explanation is simple: multiplying x and y yields a single value - not
> a range of values. And although implementations may differ on how
> accurately to represent the product of x and y - they must all agree
> that that product has in fact only one value (that is, that the
> equation has only one solution). As the C99 Standard explains:
>
> "Implementations employing wide registers have to take care to honor
> appropriate semantics. Values are independent of whether they are
> represented in a register or in memory. For example, an implicit
> "spilling" of a register is not permitted to alter the value. Also,an
> explicit "store and load" is required to round to the precision of the
> storage type. In particular, casts and assignments are required to
> perform their speci?ed conversion." [C99    5.1.2.3]
>
> In other words, if a program calculates "x*y" with extended precision
> - then the program must always calculate x*y with extended precision
> (at least in those situations where the extended precision would make
> a difference) - because multiplication is a deterministic operation
> that yields a single value (even if the precision of its
> representation is unspecified.)

I don't see how C99 is relevant.

> C++'s program execution model, moreover, requires that an
> implementation produce identical output when executing identical
> programs that perform defined operations on identical input. An
> implementation that allowed the product of two identically-valued,
> identical types to differ - would therefore not be able to adhere to C+
> +'s program execution model, because such an implementation could not
> guarantee consistent results as required.

Even assuming that the C++ execution model required a program to be
consistent between runs, a single expression like

   x*y == x*y

can still evaluate to false. It just has to do so consistently. Nothing in
the execution model prevents that.

[snip arguments based on the above]

>> d) is std::sin(x) == std::sin(x) required to yield true by the standard?
>
> A sine function (like most mathematical functions) establishes a one-
> to-one relation between the value of its input and the value of its
> output (or result).

Nit: it's not one-to-one but many-to-one.

> Therefore, calling std::sin() with identical
> arguments on the right and left side of an equality comparison - will
> yield identical values on both sides. Unless (and a robust
> implementation will handle this possibility) the value of pi were to
> change between the first and second call to std::sin().

The question is not about sin(), but about return values and whether the
standard requires them to be truncated.



Best

Kai-Uwe Bux


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]





Author: James Kanze <james.kanze@gmail.com>
Date: Fri, 13 Mar 2009 10:46:09 CST
Raw View
On Mar 12, 5:59 pm, Greg Herlihy <gre...@mac.com> wrote:
> On Mar 5, 7:06    am, Kai-Uwe Bux <jkherci...@gmx.net> wrote:
> > Clause [5/10] allows excess precision for expressions
> > involving floating point numbers. The question is, at which
> > points the standard _requires_ truncation of that excess
> > precision. Put:

> >          double x = ...
> >          double y = ...

> Essentially, an implementation may take advantage of extended
> precision floating point calculations - but only as long as
> the semantics of the program itself are not affected.

I'm not sure I understand what you mean.  A compiler can do
literally anything, as long as the observable behavior of the
program is not affected.  The statement which triggered this
discussion is: "The values of the floating operands and the
results of floating expressions may be represented in greater
precision and range than that required by the type; the types
are not changed thereby.", which seems to be an explicit licence
to allow using extended precision even if it changes the
observable behavior.

Consider the following:

     #include <iostream>

     int
     main()
     {
         double a = 1e-18 ;
         double b = 1.0 ;
         double c = -1.0 ;
         double d = (a + b) + c ;
         std::cout << d << std::endl ;
         return 0 ;
     }

And assume that the machine in question uses IEEE floating
point.  If no extended precision is used (i.e. all of the
calculations are in double), the results are 0 (which is what I
get on a Sun Sparc, which doesn't have extended precision).  If
I run it on an Intel, however (at least with g++ under Linux),
the output is "9.75782e-19".  The observable behavior has
changed because extended precision was used; my interpretation
of the sentence quoted above is that it is there to explicitly
allow this.  (In the case of an Intel, not allowing it would
slow the program down considerably.  And in most cases, like
here, result in less precise results.)

> For example, if a floating point expression like (a*b) ==
> (c*d) is true when calculated without extended precision, then
> (a*b) == (c*d) must still be true when the identical
> multiplication is performed with extended precision. A
> particular value is - in a way - independent of the precision
> of its representation.

The problem occurs when the precision starts being mixed.  If we
force the precision on one side to exactly double, without
forcing it on the other side, then the results may be false.

> >          bool equal_a ( double lhs, double rhs ) {
> >               return ( lhs == rhs );
> >          }

> >          bool equal_b ( double const & lhs, double const & rhs ) {
> >               return ( lhs == rhs );
> >          }

> > With regard to initializing parameters, consider:

> > a) The expression x*y == x*y seems to be allwed to yield
> > false since either side is allowed but neither side is
> > required to be evaluated with excess precision and the
> > built-in comparison that is part of the expression can
> > compare _different_ floating point numbers.

> No, the x*y comparison must be made between identical values.
> The explanation is simple: multiplying x and y yields a single
> value - not a range of values.

Any given multiplication yields a single value.  That value
depends on the precision, however.

> And although implementations may differ on how accurately to
> represent the product of x and y - they must all agree that
> that product has in fact only one value (that is, that the
> equation has only one solution). As the C99 Standard explains:

As the C99 standard explains.  Yes (and I was unaware of this
previously).  The C90 standard didn't state this; nor does C++
(including in the latest draft).

> "Implementations employing wide registers have to take care to
> honor appropriate semantics. Values are independent of whether
> they are represented in a register or in memory. For example,
> an implicit "spilling" of a register is not permitted to alter
> the value. Also,an explicit "store and load" is required to
> round to the precision of the storage type. In particular,
> casts and assignments are required to perform their
> specified conversion." [C99    5.1.2.3]

Of course, that particular paragraph is in a non-normative
example.  It is, however, a good indication of the intent.  More
importantly, C99 has modified the sentence I quoted above (which
I believe is taken directly from C90, although I don't have my
copy of the C90 standard handy to verify) to:

     A floating expression may be contracted, that is,
     evaluated as though it were an atomic operation, thereby
     omitting rounding errors implied by the source code and
     the expression evaluation method. The FP_CONTRACT pragma
     in <math.h> provides a way to disallow contracted
     expressions. Otherwise, whether and how expressions are
     contracted is implementation-defined.

In other words, others, in addition to myself, felt the need for
this extension of the precision to be predictable (for a given
compiler), and changed the wording to clearly require it.  (I'm
not sure I like the new wording that much, but it clearly makes
stricter requirements than the wording in the C++ standard.)

Regretfully, these changes in wording and the additional
functionality (the FP_CONTRACT pragma) haven't been adopted by
C++0x.

> In other words, if a program calculates "x*y" with extended
> precision - then the program must always calculate x*y with
> extended precision (at least in those situations where the
> extended precision would make a difference) - because
> multiplication is a deterministic operation that yields a
> single value (even if the precision of its representation is
> unspecified.)

According to C99.  C++ specifically doesn't give this guarantee.

> C++'s program execution model, moreover, requires that an
> implementation produce identical output when executing
> identical programs that perform defined operations on
> identical input.

No it doesn't.  It's easy to come up with counter examples:

     int global = 0 ;

     int
     f()
     {
         return global ;
     }

     int
     g()
     {
         ++ global ;
         return global ;
     }

     int
     main()
     {
         std::cout << f() + g() << std::endl ;
         return 0 ;
     }

A conforming implementation may output 1 or 2.  Theoretically,
at least, the value output may change even from one run of the
program to the next.

> An implementation that allowed the product of two
> identically-valued, identical types to differ - would
> therefore not be able to adhere to C+ +'s program execution
> model, because such an implementation could not guarantee
> consistent results as required.

A false premise leads to a false conclusion.  The standard
defines a parameterized set of possible observable behavior for
any given program and input; there's nothing in the standard
requiring consistency, only that the observable behavior is in
the set of allowed values.

> > b) The expression equal_b( x*y, x*y ) seems to be required
> > to yield true: the initialization of the double objects in
> > passing the parameters involves truncating the excess
> > precision. But I am not sure.

> The "x*y" expression will yield identical values in all three
> cases.  Now, whether those identical values compare equal -
> depends on the value itself (Specifically, the equality test
> will fail if x*y is a NaN, otherwise it will succeed).

> > c) About the expression equal_a( x*y, x*y ), I don't know.
> > As opposed to equal_b, the function equal_a takes its
> > parameters by value. Therefore, initializing the parameters
> > might not involve the construction of a floating point
> > object.

> > d) is std::sin(x) == std::sin(x) required to yield true by
> > the standard?

> A sine function (like most mathematical functions) establishes
> a one- to-one relation between the value of its input and the
> value of its output (or result).  Therefore, calling
> std::sin() with identical arguments on the right and left side
> of an equality comparison - will yield identical values on
> both sides. Unless (and a robust implementation will handle
> this possibility) the value of pi were to change between the
> first and second call to std::sin().

The actual question here was not so much whether the values
returned by std::sin would be identical; in practice, they will
be, whether the standard requires it or not.  The question is
whether they could have extended precision, rather than just
double precision.  If they can have extended precision, then the
compiler could spill one to memory (and in C++, truncating the
precision), while maintaining the other in extended precision.

You might want to try the following:

     for ( double e = 0.0 ; e < 1.0 ; e += 0.1 ) {
         if ( std::sin( e ) != std::sin( e ) ) {
             std::cout << "Oops! (" << e << ')' << std::endl ;
         }
     }

With g++ (4.1.0) under Linux, it outputs Oops for all of the
values in the loop (but not with VC++, although they seem to use
the same code in the loop---maybe, VC++ forces double precision
of the return value.  My own interpretation of this is that it
is illegal, but I think it is far from clear.  My interpretation
is based on reading a quite a bit into one paragraph (   8.5/12)
which is only distantly concerned.

--
James Kanze (GABI Software)             email:james.kanze@gmail.com
Conseils en informatique orient   e objet/
                    Beratung in objektorientierter Datenverarbeitung
9 place S   mard, 78210 St.-Cyr-l'   cole, France, +33 (0)1 30 23 00 34


--
[ comp.std.c++ is moderated.  To submit articles, try just posting with ]
[ your news-reader.  If that fails, use mailto:std-c++@netlab.cs.rpi.edu]
[              --- Please see the FAQ before posting. ---               ]
[ FAQ: http://www.comeaucomputing.com/csc/faq.html                      ]