This is the mail archive of the fortran@gcc.gnu.org mailing list for the GNU Fortran project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


(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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]