Re: [PATCH][4.3] Expand lround inline for x86_64/i?86 SSE math

[Brooks Moses <bmoses at stanford dot edu>]

(Added a cc: to fortran@, so as to increase the population of people who 
will complain if I say something that's incorrect!)

Richard Guenther wrote:
> On 10/19/06, Brooks Moses <bmoses@stanford.edu> wrote:
> > Toon Moene wrote:
> > > Richard Guenther wrote:
> > > > I wonder if fortran specifies round differently, as the frontend
> > > > explicitly converts NINT(x) = INT(x + ((x > 0) ? 0.5 : -0.5)).
> > >
> > > Yep - sorry, don't have the reference handy.
> >
> > For what it's worth, since I do have the reference handy, that's
> > essentially a direct translation of how the Fortran 95 standard defines
> > the NINT intrinsic.
>
> So, does it define how the x + 0.5 is carried out with respect to
> intermediate rounding before converting to INT?  Literaly writing
> the above yields 1.0 for NINT ( 0.5 - epsilon ) assuming the hardware
> rounds to nearest even for the addition.

It doesn't define this, so far as I am aware -- and, really, I'm not at 
all sure whether the standard authors intended the description to be 
interpreted as exact math, or as numerically-approximate math.

On the other hand, I just checked the Fortran 2003 standard to see if it 
was any clearer, and it looks like they decided to sidestep the issue 
entirely -- the result of NINT(A) is now defined as "the integer nearest 
A, or if there are two integers equally near A, the result is whichever 
such integer has the greater magnitude."  Since this isn't listed as an 
incompatibility between F2003 and F95, it's presumably fair to consider 
this a clarification of the intent of the F95 standard as well.

Thus, this is probably a bug in Gfortran; the following program should 
demonstrate it:

    program nint_bug
      real :: x
      x = 0.5 - epsilon(x)/4.0
      write(*,*) x - 0.5, nint(x)
      if (nint(x) /= 0) call abort
    end program

Interestingly, on i686 machines, this program fails to show the bug; 
even int(x+0.5) in this example produces 0.  This is presumably due to 
the use of extra-precision registers; if one stores x+0.5 in another 
variable (thereby pulling it out of the registers) and then calls int() 
on it, one does get an "erroneous" value of 1.

Unfortunately, I don't have a non-Intel machine handy to test this on; 
perhaps someone who has one could test it and (if appropriate) file a PR?

- Brooks