questionnable nint() behavior
Dominique Dhumieres
dominiq@lps.ens.fr
Tue Mar 21 09:22:00 GMT 2006
To make things clearer, one can replace
x = 8388609.0
y = 0.4999999701976776123046875
by
x = 1.0 + 1.0/epsilon(1.0)
y = nearest(0.5,-1.0)
floats larger or equal to 1.0/epsilon(1.0) have no fractional
part and are integers, then they are their nint().
As pointed by Steve Kargl the problem comes from
NINT(x) = INT(x + ((x > 0) ? 0.5 : -0.5))
when the rounding mode is towards nearest in which
8388609.5 is rounded to 8388610. One possible solution
would be to save the rounding flags, set the rounding
mode towards 0, then restore the rounding flags.
This could be simpler and faster than using floor().
Dominique
PS note that xlf uses the fact that adding 1.0/epsilon(1.0)
to any positive float, then substracting it returns the
integer part of any float when the rounding mode is
towards 0. Then mod() is much faster on PPC with xlf than
with g77, gfc, and g95.
More information about the Fortran
mailing list