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

[Bruno Haible]

Ben Pfaff wrote:
> There is a huge amount of redundancy in the m4 code across all
> six modules (trunc*, round*) and perhaps others.  Perhaps I
> should write a common m4 macro that just takes a parameter or two
> and avoids the redundancy.  Bruno, do you think that this is a
> good idea?

It may be a good idea in a year or so. For now, we have not yet tested
these functions on various platforms. I expect that we will discover
bugs in AIX, HP-UX, BeOS etc. - and that after the workarounds are
in place, the symmetry of the macros is gone.

> +DOUBLE
> +ROUND (DOUBLE x)
> +{
> +  if (x >= L_(0.0)) 
> +    {
> +      DOUBLE y = FLOOR (x);
> +      if (x - y >= L_(0.5))
> +        y += L_(1.0);
> +      return y;
> +    }
> +  else
> +    {
> +      DOUBLE y = FLOOR (-x);
> +      if (-x - y >= L_(0.5))
> +        y += L_(1.0);
> +      return -y;

It's less arithmetic operations if you use CEIL, knowing that
FLOOR (-x) = - CEIL(x).

I would have written the function like this:

DOUBLE
ROUND (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))
    {
      /* Add 0.5 to the absolute value.  */
      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))
    {
      /* Add 0.5 to the absolute value.  */
      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;
}

Can you evaluate the two implementations against each other (regarding
numerical correctness, number of operations, etc.) and then choose the
best parts of each?


> diff --git a/tests/test-roundf.c b/tests/test-roundf.c
> new file mode 100644
> index 0000000..56fc1db
> --- /dev/null
> +++ b/tests/test-roundf.c

> +  ASSERT (roundf (65535.999f) == 65536.0f);
> +  ASSERT (roundf (65536.0f) == 65536.0f);
> +  ASSERT (roundf (65536.001f) == 65536.0f);

In the usual IEEE single-float arithmetic, 65535.999f and 65536.001f
are the same as 65536.0f. You need to remove one of the mantissa digits
to make the test meaningful.

> diff --git a/tests/test-roundl.c b/tests/test-roundl.c
> new file mode 100644
> index 0000000..4e0ef4e
> --- /dev/null
> +++ b/tests/test-roundl.c

> +/* The Compaq (ex-DEC) C 6.4 compiler chokes on the expression 0.0 / 0.0.  */
> +#ifdef __DECC
> +static long double
> +NaN ()
> +{
> +  static long double zero = 0.0L;
> +  return zero / zero;
> +}
> +#else
> +# define NaN() (0.0L / 0.0L)
> +#endif

This workaround is not needed for 'long double'.

> +int
> +main ()
> +{

Will not work on NetBSD without the use of BEGIN_LONG_DOUBLE_ROUNDING ();
(uses module 'fpucw').

Bruno