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