bug-gnulib
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Ben Pfaff
Subject: Re: modules 'round', 'roundf', 'roundl' for review
Date: Sun, 07 Oct 2007 14:14:40 -0700
User-agent: Gnus/5.11 (Gnus v5.11) Emacs/22.1 (gnu/linux)

Bruno Haible <address@hidden> writes:

> 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.

OK.

[...Ben's implementation: call it round_blp...]
> It's less arithmetic operations if you use CEIL, knowing that
> FLOOR (-x) = - CEIL(x).

True.  I was optimizing for fewer dependencies, rather than fewer
arithmetic operations.

> I would have written the function like this:
[...Bruno's first implementation: call it round_bruno...]

> PS: A third possible implementation - if you want to rely on FLOOR and CEIL -
> is like this:
[...Bruno's second implementation: call it round_bruno2...]

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

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.  I took the liberty of
correcting a bug that had crept into round_bruno.

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

Results: round_blp and round_bruno2 match roundf in all cases.
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)?

That leaves round_blp and round_bruno2.  I like round_bruno2
better, because it is easier to understand.  

>> +  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.

You are very observant!  I will adjust the test.

>> +/* The Compaq (ex-DEC) C 6.4 compiler chokes on the expression 0.0 / 0.0.  
>> */
[...]
> This workaround is not needed for 'long double'.

Thanks, I'll drop it.

>> +int
>> +main ()
>> +{
>
> Will not work on NetBSD without the use of BEGIN_LONG_DOUBLE_ROUNDING ();
> (uses module 'fpucw').

Thanks, I'll fix that.

> And another nit:
>
>> +    HAVE_DECL_ROUND=0
>> +    AC_LIBOBJ([round])
>> +    ROUND_LIBM=
>
> If the lib/round.c substitute relies on floor, then it is wrong to set
> ROUND_LIBM to empty. In this case you need to set ROUND_LIBM="$FLOOR_LIBM"
> (and ensure that FLOOR_LIBM is already computed at this point in the
> configure file).

OK.

Here is my test program, followed by the incremental diff versus
the first version of the module that I posted:

/* Tested with glibc 2.6.1, GCC 4.1.3, on x86.  Not written for
   portability. */
#define _GNU_SOURCE 1
#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>

float
round_blp (float x)
{
  if (x >= 0.0f) 
    {
      float y = floorf (x);
      if (x - y >= 0.5f)
        y += 1.0f;
      return y;
    }
  else
    {
      float y = floorf (-x);
      if (-x - y >= 0.5f)
        y += 1.0f;
      return -y;
    }
}

/* 2^(MANT_DIG-1).  */
static const double TWO_MANT_DIG =
  /* Assume MANT_DIG <= 5 * 31.
     Use the identity
       n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5).  */
  (float) (1U << ((FLT_MANT_DIG - 1) / 5))
  * (float) (1U << ((FLT_MANT_DIG - 1 + 1) / 5))
  * (float) (1U << ((FLT_MANT_DIG - 1 + 2) / 5))
  * (float) (1U << ((FLT_MANT_DIG - 1 + 3) / 5))
  * (float) (1U << ((FLT_MANT_DIG - 1 + 4) / 5));

float
round_bruno (float 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 float z = x;
  volatile float y;

  if (z > 0.0f)
    {
      /* Add 0.5 to the absolute value.  */
      y = z += 0.5f;
      /* 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 -= 1.0f;
    }
  else if (z < 0.0f)
    {
      /* Add 0.5 to the absolute value.  */
      y = z -= 0.5f;
      /* 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 += 1.0f;
    }
  return z;
}

float
round_bruno2 (float x) 
{
  if (x >= 0.0f)
    return floorf (x + 0.5f);
  else
    return ceilf (x - 0.5f);
}

int
main (void) 
{
  uint32_t x;

  for (x = 0; x < UINT32_MAX; x++) 
    {
      float f, g, h;

      if (!(x << 7))
        {
          putchar ('.');
          fflush (stdout);
        }

      assert (sizeof x == sizeof f);
      memcpy (&f, &x, sizeof f);
      g = roundf (f);
      h = round_blp (f);
      if ((isnan (g) && isnan (h)) || g == h)
        continue;

      printf ("round %f(%a) = %f(%a) or %f(%a)?\n", f, f, g, g, h, h);
      break;
    }
  return 0;
}

diff --git a/lib/round.c b/lib/round.c
index 8baf119..3f5ec9d 100644
--- a/lib/round.c
+++ b/lib/round.c
@@ -26,16 +26,19 @@
 #ifdef USE_LONG_DOUBLE
 # define ROUND roundl
 # define FLOOR floorl
+# define CEIL ceill
 # define DOUBLE long double
 # define L_(literal) literal##L
 #elif ! defined USE_FLOAT
 # define ROUND round
 # define FLOOR floor
+# define CEIL ceil
 # define DOUBLE double
 # define L_(literal) literal
 #else /* defined USE_FLOAT */
 # define ROUND roundf
 # define FLOOR floorf
+# define CEIL ceilf
 # define DOUBLE float
 # define L_(literal) literal##f
 #endif
@@ -43,18 +46,8 @@
 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;
-    }
+  if (x >= L_(0.0))
+    return FLOOR (x + L_(0.5));
   else
-    {
-      DOUBLE y = FLOOR (-x);
-      if (-x - y >= L_(0.5))
-        y += L_(1.0);
-      return -y;
-    }
+    return CEIL (x - L_(0.5));
 }
diff --git a/m4/round.m4 b/m4/round.m4
index 76ab809..86ca9fa 100644
--- a/m4/round.m4
+++ b/m4/round.m4
@@ -39,9 +39,38 @@ AC_DEFUN([gl_FUNC_ROUND],
       ROUND_LIBM=
     fi
   else
+    dnl We need to use our substitute round, which depends on floor
+    dnl and ceil.  Test whether floor can be used without libm.
+    dnl (We assume that floor and ceil are in the same library.)
+    FLOOR_LIBM=?
+    AC_TRY_LINK([
+       #ifndef __NO_MATH_INLINES
+       # define __NO_MATH_INLINES 1 /* for glibc */
+       #endif
+       #include <math.h>
+       double x;],
+      [x = floor(x);],
+      [FLOOR_LIBM=])
+    if test "$FLOOR_LIBM" = "?"; then
+      save_LIBS="$LIBS"
+      LIBS="$LIBS -lm"
+      AC_TRY_LINK([
+         #ifndef __NO_MATH_INLINES
+         # define __NO_MATH_INLINES 1 /* for glibc */
+         #endif
+         #include <math.h>
+         double x;],
+        [x = floor(x);],
+        [FLOOR_LIBM="-lm"])
+      LIBS="$save_LIBS"
+    fi
+    if test "$FLOOR_LIBM" = "?"; then
+      FLOOR_LIBM=
+    fi
+
     HAVE_DECL_ROUND=0
     AC_LIBOBJ([round])
-    ROUND_LIBM=
+    ROUND_LIBM=$FLOOR_LIBM
   fi
   AC_SUBST([HAVE_DECL_ROUND])
   AC_SUBST([ROUND_LIBM])
diff --git a/modules/roundl-tests b/modules/roundl-tests
index 9bed031..957c934 100644
--- a/modules/roundl-tests
+++ b/modules/roundl-tests
@@ -2,6 +2,8 @@ Files:
 tests/test-roundl.c
 
 Depends-on:
+fpucw
+isnanl-nolibm
 
 configure.ac:
 
diff --git a/tests/test-roundf.c b/tests/test-roundf.c
index 56fc1db..6311e3a 100644
--- a/tests/test-roundf.c
+++ b/tests/test-roundf.c
@@ -63,9 +63,9 @@ main ()
   ASSERT (roundf (2.5f) == 3.0f);
   ASSERT (roundf (1.999f) == 2.0f);
   ASSERT (roundf (2.0f) == 2.0f);
-  ASSERT (roundf (65535.999f) == 65536.0f);
+  ASSERT (roundf (65535.99f) == 65536.0f);
   ASSERT (roundf (65536.0f) == 65536.0f);
-  ASSERT (roundf (65536.001f) == 65536.0f);
+  ASSERT (roundf (65536.01f) == 65536.0f);
   ASSERT (roundf (2.341e31f) == 2.341e31f);
   /* Negative numbers.  */
   ASSERT (roundf (-0.3f) == 0.0f);
@@ -76,9 +76,9 @@ main ()
   ASSERT (roundf (-2.5f) == -3.0f);
   ASSERT (roundf (-1.999f) == -2.0f);
   ASSERT (roundf (-2.0f) == -2.0f);
-  ASSERT (roundf (-65535.999f) == -65536.0f);
+  ASSERT (roundf (-65535.99f) == -65536.0f);
   ASSERT (roundf (-65536.0f) == -65536.0f);
-  ASSERT (roundf (-65536.001f) == -65536.0f);
+  ASSERT (roundf (-65536.01f) == -65536.0f);
   ASSERT (roundf (-2.341e31f) == -2.341e31f);
   /* Infinite numbers.  */
   ASSERT (roundf (1.0 / 0.0f) == 1.0 / 0.0f);
diff --git a/tests/test-roundl.c b/tests/test-roundl.c
index 4e0ef4e..1495eae 100644
--- a/tests/test-roundl.c
+++ b/tests/test-roundl.c
@@ -25,6 +25,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 
+#include "fpucw.h"
+#include "isnanl-nolibm.h"
+
 #define ASSERT(expr) \
   do                                                                        \
     {                                                                       \
@@ -36,21 +39,13 @@
     }                                                                       \
   while (0)
 
-/* 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
-
 int
 main ()
 {
+  DECL_LONG_DOUBLE_ROUNDING
+
+  BEGIN_LONG_DOUBLE_ROUNDING ();
+
   /* Zero.  */
   ASSERT (roundl (0.0L) == 0.0L);
   ASSERT (roundl (-0.0L) == 0.0L);
@@ -84,7 +79,7 @@ main ()
   ASSERT (roundl (1.0 / 0.0L) == 1.0 / 0.0L);
   ASSERT (roundl (-1.0 / 0.0L) == -1.0 / 0.0L);
   /* NaNs.  */
-  ASSERT (isnan (roundl (NaN ())));
+  ASSERT (isnanl (roundl (0.0L / 0.0L)));
 
   return 0;
 }

-- 
Ben Pfaff 
http://benpfaff.org




reply via email to

[Prev in Thread] Current Thread [Next in Thread]