questionnable nint() behavior
Timothy C Prince
tprince@myrealbox.com
Tue Mar 21 00:52:00 GMT 2006
-----Original Message-----
From: Steve Kargl <sgk@troutmask.apl.washington.edu>
To: Dominique Dhumieres <dominiq@lps.ens.fr>
Date: Mon, 20 Mar 2006 16:35:03 -0800
Subject: Re: questionnable nint() behavior
On Mon, Mar 20, 2006 at 10:24:12PM +0100, Dominique Dhumieres wrote:
> I have found in my file the following program
>
> program main
> real x, y
> x = 8388609.0
> y = 0.4999999701976776123046875
> print '(A6,F9.1,A3,I8)', 'nint (', x, ') =', nint (x)
> print '(A6,F10.8,A3,I2,A18,L1)', 'nint (', y, ') =', nint (y), &
> & ', where y < 0.5 = ', y < 0.5
> end
>
> xlf, ifc, and g95 return
>
> nint (8388609.0) = 8388609
> nint (0.49999997) = 0, where y < 0.5 = T
>
> IMHO the correct result, g77, pgf, and gfc return
>
> nint (8388609.0) = 8388610
> nint (0.49999997) = 1, where y < 0.5 = T
>
> as if nint(a) was implemented as int(a+0.5) without
> taking care of rounding mode (rounding to nearest
> rounds to even on tie). I don't remember where I borrowed
> the code, but I thought the problem was fixed in gfortran.
>
In trans-intrinsics.c, you'll find
/* This is needed because the gcc backend only implements FIX_TRUNC_EXPR
NINT(x) = INT(x + ((x > 0) ? 0.5 : -0.5)). */
static tree
build_round_expr (stmtblock_t * pblock, tree arg, tree type)
A much better algorithm may be found in FreeBSD's round[fl].c. :-)
--
______________________________________
We could get away with this, by promoting float to double, and double to long double, as happens already in gfortran i386 for some -march= values. g77 (libf2c) had a more careful and slow version, using floor().
Tim Prince
More information about the Fortran
mailing list