wrong results from generic rint function (s_rint.c)

[Ruppert <dieter_ruppert at siemens dot com>]

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