wrong results from generic rint function (s_rint.c)
Ruppert
dieter_ruppert@siemens.com
Fri Mar 5 10:33:00 GMT 2010
Hi,
there seems to be a bug in the generic C implementation of rint(double)
in s_rint.c which manifests itself in wrong rounding behaviour for certain
input values.
Example:
rint(262144.75) should return 262145.0
(with the normal "round to nearest" behaviour), instead it returns 262144.0,
as if rounding towards zero. Slightly higher or lower input values (262144.749
or 262144.751, for example) yield the correct result 262145.0.
I could reproduce this with a small test program on x86/Linux, Solaris/Sparc,
Solaris/x86 and on a PowerPC Mac, so this is probably independent from
processor or FPU properties. I took s_rint.c from newlib-1.18, but this
implementation is the same for all older versions.
I poked around a bit in this "bit twiddling" implementation of rint, but
could not spot the problem. It is, however, apparent that (at least) the
input values
262144.75, 262146.75, 262148.75, ... 524286.75
result in this false rounding behaviour. These values have
- an even integral part from 2^18 up to and including 2^19 - 2
- a fractional part of 0.75
For all these values, the low order 32 bits of the mantissa are zero, and
the three least significant bits of the upper part of the mantissa are 011.
The exponent (18) is such that the final rounding (by adding and subtracting
from 2^52) shifts these three bits to the right so that the '0' bit becomes
the least significant bit of the mantissa.
I would suggest that somebody familiar with the rationale behind all these
shift and mask operations in the implementation takes a look at this. I assume
that all newlib uses which don't employ a "fastmath" implementation might
experience this strange bug.
Should I file a bug report somewhere and/or provide my test program?
Regards
Dieter Ruppert
More information about the Newlib
mailing list