This is the mail archive of the
mailing list for the GNU Fortran project.
Re: [PATCH][4.3] Expand lround inline for x86_64/i?86 SSE math
- From: Steve Kargl <sgk at troutmask dot apl dot washington dot edu>
- To: Brooks Moses <bmoses at stanford dot edu>
- Cc: fortran at gcc dot gnu dot org, gcc-patches at gcc dot gnu dot org
- Date: Thu, 19 Oct 2006 15:34:23 -0700
- Subject: Re: [PATCH][4.3] Expand lround inline for x86_64/i?86 SSE math
- References: <Pine.LNX.email@example.com> <EEA92755-FDBD-4FF8-AAAB-ACBB9461A6F7@adacore.com> <9645931.1161211060458.JavaMail.firstname.lastname@example.org> <4537C969.email@example.com> <4537DC19.firstname.lastname@example.org> <email@example.com> <4537F091.firstname.lastname@example.org>
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 <email@example.com> 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