Avoid floating-point compares on Stem:head_positions; issue 4310

Avoid floating-point compares on Stem:head_positions; issue 4310

On the cross-compiled windows executable, head_positions().is_empty()
seems to sometimes return true for a stem with a single note.
head_positions() returns something like [3.00... , 3.00...]

[lemzwerg, 2015-03-11 05:19:34]

LGTM.

[dak, 2015-03-11 09:10:29]

_How_ do we arrive at different versions of 3.0 here, or at GCC considering 3.0
> 3.0 ?  While doing floating point comparisons is fraught with its own
problems, small integers _are_ exactly representable in IEEE floating point
numbers regardless of their precision.

The fix certainly makes sense in that it makes the code more straightforward. 
However, the condition triggering the bad behavior is likely to bite us
elsewhere with random symptoms.

So we should find the cause and fix it, even though in this particular case the
bad effect can be made to go away by sanitizing the code.

[lemzwerg, 2015-03-11 09:40:11]

Any idea how to find it?

[Keith, 2015-03-11 16:36:44]

On Wed, 11 Mar 2015 02:10:29 -0700, <dak@gnu.org> wrote:

> _How_ do we arrive at different versions of 3.0 here, or at GCC
> considering 3.0 > 3.0 ?

Probably we get different numbers that are close to 3.00... or 2.99...

The head positions are individually multiplied by staff-space in
Grob::get_offset() and then divided by staff-space
Staff_symbol_referencer::internal_get_position() and at magstep=2 staff-space is
2^(1/3)

I think a standard-conforming implementation should give bit-identical
floating-point results for these parallel computations, but suspect some
compiler option is set to allow non-conformance (retaining the most-recent value
in register with its guard bits intact for the comparison, maybe).

In general, I think it wisest for us to code in a way that allows floating point
to be inexact. The broader change that I considered was to count only only [+inf
-inf] as the 'empty' interval in is_empty() as we considered doing before.

> https://codereview.appspot.com/217720043/

[dak, 2015-03-11 18:14:37]

On 2015/03/11 16:36:44, Keith wrote:
> On Wed, 11 Mar 2015 02:10:29 -0700, <mailto:dak@gnu.org> wrote:
> 
> > _How_ do we arrive at different versions of 3.0 here, or at GCC
> > considering 3.0 > 3.0 ?
> 
> Probably we get different numbers that are close to 3.00... or 2.99...
> 
> The head positions are individually multiplied by staff-space in
> Grob::get_offset() and then divided by staff-space
> Staff_symbol_referencer::internal_get_position() and at magstep=2 staff-space
is
> 2^(1/3)
> 
> I think a standard-conforming implementation should give bit-identical
> floating-point results for these parallel computations,

No, not at all.  In C++, any computation may be carried out at any precision at
least as high as that of the operands.  Only when the results are stored in a
variable (or possibly explicitly cast to some fp type) is it guaranteed that the
precision is not higher than specified.

If we indeed do the equivalent of n*x/x, we are no longer guaranteed that
starting from an integer n we arrive at an integer, and when we do this in two
separate code paths, one calculation may end up as an integer while another
doesn't.

If we want to be safe about preserving integers, we must not scale back and
forth with some non-integer factor.

That's a not really fun part of the C standard.  We might use some gcc option
like -fstore-floats or something similar to avoid that if we can't clean up our
code sufficiently.  But the scaling back-and-forth seems bad in its own right.

[Keith, 2015-03-12 03:25:20]

On 2015/03/11 18:14:37, dak wrote:
> On 2015/03/11 16:36:44, Keith wrote:
> > I think a standard-conforming implementation should give bit-identical
> > floating-point results for these parallel computations,
> 
> No, not at all.  In C++, any computation may be carried out at any precision
at
> least as high as that of the operands.  Only when the results are stored in a
> variable (or possibly explicitly cast to some fp type) is it guaranteed that
the
> precision is not higher than specified.
> 
> If we indeed do the equivalent of n*x/x, we are no longer guaranteed that
> starting from an integer n we arrive at an integer, and when we do this in two
> separate code paths, one calculation may end up as an integer while another
> doesn't.
> 

Oh.  That's probably it, then.