RE: wrong results from generic rint function (s_rint.c)

["Howland Craig D (Craig)" <howland at LGSInnovations dot com>]

I have reproduced the problem.  I'll take a look this weekend as to
what might be wrong (assuming that nobody else jumps in and takes care
of it sooner).
Craig 

-----Original Message-----
From: newlib-owner@sourceware.org [mailto:newlib-owner@sourceware.org]
On Behalf Of Ruppert
Sent: Friday, March 05, 2010 5:33 AM
To: newlib@sources.redhat.com
Cc: dieter_ruppert@siemens.com
Subject: wrong results from generic rint function (s_rint.c)


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