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

[Steve Kargl]

On Thu, Oct 19, 2006 at 02:39:29PM -0700, Brooks Moses wrote:
> (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.

The standard doesn't defined anything about intermediate rounding.

It states that NINT(X) = INT(X+0.5) for X > 0.  You then
find INT(A) = 0 for |A| < 1.  A quick scan over the section
on numeric binary operators did not reveal the normal 
weasel words, but I suspect that the result of X+0.5 is a
processor-dependent approximation and round-to-nearest
meets that criterium.

OTOH, I think it is a bug because we can choose a better
algorithm to determine if x+0.5 is truly less than 1.
See for example round[f] in intrinsics/c99_functions.c

-- 
Steve