Re: modules 'round', 'roundf', 'roundl' for review

[Bruno Haible]

Hello Ben,

> I'm appending a program that can be used to test the correctness
> of each of the three implementations, by evaluating them for
> every possible 32-bit value of float and checking the results
> against the system roundf's result.

Well done!

This test program can be made faster, by varying only the
highest and the lowest bits; the bits in the middle can be set to
00....00 or 11...11 without weakening the check. This allows to
generalize the test to 'double' (64-bit).

When applying this test program to 'floor', 'ceil', 'trunc', it
uncovered bugs there. I've fixed them now:
http://git.sv.gnu.org/gitweb/?p=gnulib.git;a=commitdiff;h=3db9b0d802e34ec28fc924076ae4c196bd7d23a9
http://git.sv.gnu.org/gitweb/?p=gnulib.git;a=commitdiff;h=33d6a0e5ca06c6d0a8b0fddba8c7f41bbd6794ae
http://git.sv.gnu.org/gitweb/?p=gnulib.git;a=commitdiff;h=c7f33ef75bc531bc875e34019324811962b001d1

> I took the liberty of
> correcting a bug that had crept into round_bruno.

Thanks, this was a mistake.

> (By the way, is it intentional that TWO_MANT_DIG is always of type
> "double"?  Should it be type DOUBLE instead?)

A mistake as well.

> round_bruno does not: it fails in many cases, such as
> 0.5-2*epsilon:
>         round 0.500000(0x1.fffffep-2) = 0.000000(0x0p+0) or 1.000000(0x1p+0)?

Here's a better version of both.

==============================================================================

DOUBLE
round_bruno1 (DOUBLE x)
{
  /* The use of 'volatile' guarantees that excess precision bits are dropped
     at each addition step and before the following comparison at the caller's
     site.  It is necessary on x86 systems where double-floats are not IEEE
     compliant by default, to avoid that the results become platform and 
compiler
     option dependent.  'volatile' is a portable alternative to gcc's
     -ffloat-store option.  */
  volatile DOUBLE y = x;
  volatile DOUBLE z = y;

  if (z > L_(0.0))
    {
      /* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1).  */
      if (z < L_(0.5))
        z = L_(0.0);
      /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1.  */
      else if (z < TWO_MANT_DIG)
        {
          /* Add 0.5 to the absolute value.  */
          y = z += L_(0.5);
          /* Round to the next integer (nearest or up or down, doesn't matter). 
 */
          z += TWO_MANT_DIG;
          z -= TWO_MANT_DIG;
          /* Enforce rounding down.  */
          if (z > y)
            z -= L_(1.0);
        }
    }
  else if (z < L_(0.0))
    {
      /* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)).  */
      if (z > - L_(0.5))
        z = L_(0.0);
      /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1.  */
      else if (z > -TWO_MANT_DIG)
        {
          /* Add 0.5 to the absolute value.  */
          y = z -= L_(0.5);
          /* Round to the next integer (nearest or up or down, doesn't matter). 
 */
          z -= TWO_MANT_DIG;
          z += TWO_MANT_DIG;
          /* Enforce rounding up.  */
          if (z < y)
            z += L_(1.0);
        }
    }
  return z;
}

DOUBLE
round_bruno2 (DOUBLE x)
{
  if (x >= L_(0.0))
    {
      /* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1).  */
      if (x < L_(0.5))
        return L_(0.0);
      if (x < TWO_MANT_DIG)
        return FLOOR (x + L_(0.5));
    }
  else
    {
      /* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)).  */
      if (x > - L_(0.5))
        return L_(0.0);
      if (x > - TWO_MANT_DIG)
        return CEIL (x - L_(0.5));
    }
  return x;
}

==============================================================================

What to do with these two functions? I would suggest to
1) test them both as reference implementations, like done in test-truncf2.c
   and test-trunc2.c.
2) use round_bruno1 as the implementation in lib/round.c if the corresponding
   functions floor, ceil (or floorf, ceilf, or floorl, ceill) are not
   declared. If they are declared, one can assume that they are efficiently
   implemented, and then round_blp looks like quite efficient.

Bruno