cvs commit: src/lib/msun/src s_rint.c

[Bruce Evans]

bde         2005-12-03 07:23:30 UTC

  FreeBSD src repository

  Modified files:
    lib/msun/src         s_rint.c 
  Log:
  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).
  
  Revision  Changes    Path
  1.12      +9 -0      src/lib/msun/src/s_rint.c