wrong results from generic rint function (s_rint.c)

[Jeff Johnston]

Thanks Craig.  Patch checked in.  Would your test be of any benefit in 
the testsuite directory or is it just the problem scenario?

-- Jeff J.

On 03/08/2010 11:25 AM, Howland Craig D (Craig) wrote:
> Patch attached.  Tested on i386 (i686, to be more precise).
>
> The source of the problem is due to an optional step that limits how
> many
> bits are set in the fractional part before the primary rounding step is
> done.  When the integer-part/fraction-part boundary crossed the MSB/LSB
> boundary in the bit-fiddling code, the error occurred when 18 bits were
> in the integral part.  The value of the fraction that could cause a
> problem was basically any in which the .5 bit and any bit .25 or smaller
> was also set.  (e.g. .5625, 1.75, 1.9375).
>
> The apparent intent of this step of limiting the number of bits in the
> fraction is to avoid double-rounding errors in implementations which
> use guard bits.  (Although I can't figure out how it actually would
> make any difference.  I kept it, though, in case I am missing
> something.)
>
> This error was simply due to a missing "else i1 = 0;".
>
> However, the code also contained an inconsistency, as there are really
> 2 places where the number of bits in the fraction are limited.  The
> first is when the int/frac boundary is near the MSB/LSB boundary, where
> the error was.  It was setting 0b.001.  The second place was when it
> falls only in the MSB, where it sets 0b.01.  So in addition to fixing
> the error, I also changed the two fraction-limit steps to be the same,
> going with the non-errored MSB-only step of 0b.01.  Not surprisingly,
> sf_rint.c contained an inconsistency in this step when compared against
> s_rint.c, so I have adjusted that also.  (Likely no bug in rintf(), but
> it's better to have the rint() and rintf() algorithms exactly match.)
>
> While fixing the problems, I also put in a few (although not complete)
> comments in the code that help explain what is being done.
>
> By the way, the FreeBSD code made an attempt to fix this, but I think
> missed a case--with errors possibly happening with a 17-bit integer
> part.
> Since I'm not on the FreeBSD mailing list, if someone on this list is on
> the FreeBSD list, it might be good to pass this along.
>
> Craig
>
> -----Original Message-----
> From: Howland Craig D (Craig)
> Sent: Friday, March 05, 2010 10:40 AM
> To: 'Ruppert'
> Cc: newlib@sources.redhat.com
> Subject: RE: wrong results from generic rint function (s_rint.c)
>
> 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.
>
> ...
>
> Regards
> Dieter Ruppert
>
>
>