Fixed rint(x) in the following cases:

(1) In round-to-nearest mode, on all machines, fdlibm rint() never
    worked for |x| = n+0.75 where n is an even integer between 262144
    and 524286 inclusive (2*131072 cases).  To avoid double rounding
    on some machines, we begin by adjusting x to a value with the 0.25
    bit not set, essentially by moving the 0.25 bit to a lower bit
    where it works well enough as a guard, but we botched the adjustment
    when log2(|x|) == 18 (2*2**52 cases) and ended up just clearing
    the 0.25 bit then.  Most subcases still worked accidentally since
    another lower bit serves as a guard.  The case of odd n worked
    accidentally because the rounding goes the right way then.  However,
    for even n, after mangling n+0.75 to 0.5, rounding gives n but the
    correct result is n+1.
(2) In round-towards-minus-infinity mode, on all machines, fdlibm rint()
    never for x = n+0.25 where n is any integer between -524287 and
    -262144 inclusive (262144 cases).  In these cases, after mangling
    n+0.25 to n, rounding gives n but the correct result is n-1.
(3) In round-towards-plus-infinity mode, on all machines, fdlibm rint()
    never for x = n+0.25 where n is any integer between 262144 and
    524287 inclusive (262144 cases).  In these cases, after mangling
    n+0.25 to n, rounding gives n but the correct result is n+1.

A variant of this bug was fixed for the float case in rev.1.9 of s_rintf.c,
but the analysis there is incomplete (it only mentions (1)) and the fix
is buggy.

Example of the problem with double rounding: rint(1.375) on a machine
which evaluates double expressions with just 1 bit of extra precision
and is in round-to-nearest mode.  We evaluate the result using
(double)(2**52 + 1.375) - 2**52.  Evaluating 2**52 + 1.375 in (53+1) bit
prcision gives 2**52 + 1.5 (first rounding).  (Second) rounding of this
to double gives 2**52 + 2.0.  Subtracting 2**52 from this gives 2.0 but
we want 1.0.  Evaluating 2**52 + 1.375 in double precision would have
given the desired intermediate result of 2**52 + 1.0.

The double rounding problem is relatively rare, so the botched adjustment
can be fixed for most machines by removing the entire adjustment.  This
would be a wrong fix (using it is 1 of the bugs in rev.1.9 of s_rintf.c)
since fdlibm is supposed to be generic, but it works in the following cases:
- on all machines that evaluate double expressions in double precision,
  provided either long double has the same precision as double (alpha,
  and i386's with precision forced to double) or my earlier fix to use
  a long double 2**52 is modified to avoid using long double precision.
- on all machines that evaluate double expressions in many more than 11
  bits of extra precision.  The 1 bit of extra precision in the example
  is the worst case.  With N bits of extra precision, it sufices to
  adjust the bit N bits below the 0.5 bit.  For N >= about 52 there is
  no such bit so the adjustment is both impossible and unnecessary.  The
  fix in rev.1.9 of s_rintf.c apparently depends on corresponding magic
  in float precision: on all supported machines N is either 0 or >= 24,
  so double rounding doesn't occur in practice.
- on all machines that don't use fdlibm rint*() (i386's).
So under FreeBSD, the double rounding problem only affects amd64 now, but
should only affect i386 in future (when double expressions are evaluated
in long double precision).

1 files changed, 9 insertions, 0 deletions

diff --git a/lib/msun/src/s_rint.c b/lib/msun/src/s_rint.c
index 9a2b42e..2d9b672 100644
--- a/lib/msun/src/s_rint.c
+++ b/lib/msun/src/s_rint.c
@@ -65,7 +65,16 @@ rint(double x)
 		if(((i0&i)|i1)==0) return x; /* x is integral */
 		i>>=1;
 		if(((i0&i)|i1)!=0) {
+		    /*
+		     * Some bit is set after the 0.5 bit.  To avoid the
+		     * possibility of errors from double rounding in
+		     * w = TWO52[sx]+x, adjust the 0.25 bit to a lower
+		     * guard bit.  We do this for all j0<=51.  The
+		     * adjustment is trickiest for j0==18 and j0==19
+		     * since then it spans the word boundary.
+		     */
 		    if(j0==19) i1 = 0x40000000; else
+		    if(j0==18) i1 = 0x80000000; else
 		    i0 = (i0&(~i))|((0x20000)>>j0);
 		}
 	    }