Re: Problem with compl() in gawk-3.1.3

[Paul Eggert]

Thanks for reporting that interesting bug in Gawk.
Here is a proposed patch, relative to gawk 3.1.3.

2003-09-27  Paul Eggert  <address@hidden>

        Fix a bug reported by Mike Romaniw <address@hidden>
        to bug-gnu-utils today: compl(compl(0xf0f)) returned 0xfff on hosts
        with 64-bit uintmax_t and 64-bit IEEE-764 double, due to rounding
        errors.
        * doc/gawk.texi (Bitwise Functions): Leading nonzero bits are
        removed in order to fit the result into a C 'double' without rounding
        error.
        * builtin.c: Include <float.h> if available.
        (FLT_RADIX, FLT_MANT_DIG, DBL_MANT_DIG): Define if not already defined.
        (AWKSMALL_MANT_DIG, AWKNUM_MANT_DIG, AWKNUM_FRACTION_BITS): New macros.
        (tmp_integer): New function.
        (do_lshift, do_rshift, do_and, do_or, do_xor, do_compl): Use them.

===================================================================
RCS file: doc/gawk.texi,v
retrieving revision 3.1.3.0
retrieving revision 3.1.3.1
diff -pu -r3.1.3.0 -r3.1.3.1
--- doc/gawk.texi       2003/07/04 17:40:47     3.1.3.0
+++ doc/gawk.texi       2003/09/28 06:14:51     3.1.3.1
@@ -13906,7 +13906,9 @@ Return the value of @var{val}, shifted r
 
 For all of these functions, first the double-precision floating-point value is
 converted to the widest C unsigned integer type, then the bitwise operation is
-performed and then the result is converted back into a C @code{double}. (If
+performed.  If the result cannot be represented exactly as a C @code{double},
+leading nonzero bits are removed one by one until it can be represented
+exactly.  The result is then converted back into a C @code{double}.  (If
 you don't understand this paragraph, don't worry about it.)
 
 Here is a user-defined function
===================================================================
RCS file: builtin.c,v
retrieving revision 3.1.3.0
retrieving revision 3.1.3.1
diff -pu -r3.1.3.0 -r3.1.3.1
--- builtin.c   2003/07/06 22:08:08     3.1.3.0
+++ builtin.c   2003/09/28 06:14:51     3.1.3.1
@@ -75,9 +75,22 @@ extern int output_is_tty;
 
 static NODE *sub_common P((NODE *tree, long how_many, int backdigs));
 
+#ifdef STDC_HEADERS
+#include <float.h>
+#endif
+/* Assume IEEE-754 arithmetic on pre-C89 hosts.  */
+#ifndef FLT_RADIX
+#define FLT_RADIX 2
+#endif
+#ifndef FLT_MANT_DIG
+#define FLT_MANT_DIG 24
+#endif
+#ifndef DBL_MANT_DIG
+#define DBL_MANT_DIG 53
+#endif
+
 #ifdef _CRAY
 /* Work around a problem in conversion of doubles to exact integers. */
-#include <float.h>
 #define Floor(n) floor((n) * (1.0 + DBL_EPSILON))
 #define Ceil(n) ceil((n) * (1.0 + DBL_EPSILON))
 
@@ -2397,6 +2410,56 @@ sgfmt(char *buf, /* return buffer; assum
 }
 #endif /* GFMT_WORKAROUND */
 
+/*
+ * The number of base-FLT_RADIX digits in an AWKNUM fraction, assuming
+ * that AWKNUM is not long double.
+ */
+#define AWKSMALL_MANT_DIG \
+  (sizeof (AWKNUM) == sizeof (double) ? DBL_MANT_DIG : FLT_MANT_DIG)
+
+/*
+ * The number of base-FLT_DIGIT digits in an AWKNUM fraction, even if
+ * AWKNUM is long double.  Don't mention 'long double' unless
+ * LDBL_MANT_DIG is defined, for the sake of ancient compilers that
+ * lack 'long double'.
+ */
+#ifdef LDBL_MANT_DIG
+#define AWKNUM_MANT_DIG \
+  (sizeof (AWKNUM) == sizeof (long double) ? LDBL_MANT_DIG : AWKSMALL_MANT_DIG)
+#else
+#define AWKNUM_MANT_DIG AWKSMALL_MANT_DIG
+#endif
+
+/*
+ * The number of bits in an AWKNUM fraction, assuming FLT_RADIX is
+ * either 2 or 16.  IEEE and VAX formats use radix 2, and IBM
+ * mainframe format uses radix 16; we know of no other radices in
+ * practical use.
+ */
+#if FLT_RADIX != 2 && FLT_RADIX != 16
+Please port the following code to your weird host;
+#endif
+#define AWKNUM_FRACTION_BITS (AWKNUM_MANT_DIG * (FLT_RADIX == 2 ? 1 : 4))
+ 
+/* tmp_integer - Convert an integer to a temporary number node.  */
+
+static NODE *
+tmp_integer(uintmax_t n)
+{
+  /*
+   * If uintmax_t is so wide that AWKNUM cannot represent all its
+   * values, strip leading nonzero bits of integers that are so large
+   * that they cannot be represented exactly as AWKNUMs, so that their
+   * low order bits are represented exactly, without rounding errors.
+   * This is more desirable in practice, since it means the user sees
+   * integers that are the same width as the AWKNUM fractions.
+   */
+  if (AWKNUM_FRACTION_BITS < CHAR_BIT * sizeof n)
+    n &= ((uintmax_t) 1 << AWKNUM_FRACTION_BITS) - 1;
+
+  return tmp_number((AWKNUM) n);
+}
+
 /* do_lshift --- perform a << operation */
 
 NODE *
@@ -2431,7 +2494,7 @@ do_lshift(NODE *tree)
        ushift = (uintmax_t) shift;
 
        res = uval << ushift;
-       return tmp_number((AWKNUM) res);
+       return tmp_integer(res);
 }
 
 /* do_rshift --- perform a >> operation */
@@ -2468,7 +2531,7 @@ do_rshift(NODE *tree)
        ushift = (uintmax_t) shift;
 
        res = uval >> ushift;
-       return tmp_number((AWKNUM) res);
+       return tmp_integer(res);
 }
 
 /* do_and --- perform an & operation */
@@ -2503,7 +2566,7 @@ do_and(NODE *tree)
        uright = (uintmax_t) right;
 
        res = uleft & uright;
-       return tmp_number((AWKNUM) res);
+       return tmp_integer(res);
 }
 
 /* do_or --- perform an | operation */
@@ -2538,7 +2601,7 @@ do_or(NODE *tree)
        uright = (uintmax_t) right;
 
        res = uleft | uright;
-       return tmp_number((AWKNUM) res);
+       return tmp_integer(res);
 }
 
 /* do_xor --- perform an ^ operation */
@@ -2573,7 +2636,7 @@ do_xor(NODE *tree)
        uright = (uintmax_t) right;
 
        res = uleft ^ uright;
-       return tmp_number((AWKNUM) res);
+       return tmp_integer(res);
 }
 
 /* do_compl --- perform a ~ operation */
@@ -2600,7 +2663,7 @@ do_compl(NODE *tree)
 
        uval = (uintmax_t) d;
        uval = ~ uval;
-       return tmp_number((AWKNUM) uval);
+       return tmp_integer(uval);
 }
 
 /* do_strtonum --- the strtonum function */