Re: "seq .1 .1" would mistakenly generate no output on FreeBSD 6.1

[Paul Eggert]

Jim Meyering wrote:

>  (gdb) p last
>  $4 = 0.10000000000000000000135525271560688
>  (gdb) p x
>  $5 = 0.1000000000000000055511151231257827

That's odd, since X was set via this assignment:

  x = first + i * step;

and I is zero, which means X should equal FIRST, and (with "seq .1 .1")
FIRST should equal LAST.  All the values in question have the type
"long double" so there shouldn't be any rounding error here.

I suspect the problem came up because FreeBSD has what I like to call
"varying width" floating point.  That is, "long double" sometimes has
a 53-bit fraction, and sometimes a 64-bit fraction, depending on which
temporary happens to hold the "long double" value.  This would explain
the problem, since the closest approximation to 0.1 with a 53-bit
fraction (IEEE double) is exactly
0.1000000000000000055511151231257827021181583404541015625, and with a
64-bit fraction (Intel extended) it's exactly
0.1000000000000000000013552527156068805425093160010874271392822265625,
and these numbers match the values shown above (modulo the finite
print width of GDB).

"Varying width" floating point conforms to C89, but it's a real hassle
to code portably to.  For example, if you have code like this:

        long double x = (something);
        long double y = x;

you cannot assume that y == x.  Here's another, more-drastic example:

        long double x = (something);
        long double y = (something_else);
        assert (x == y || x != y);

The assertion might fail, since the first use of X might have a 64-bit
fraction, but the second might have a 53-bit fraction.

I looked at the revised print_numbers function and found what I think
is one or two other instances of similar problems.  I hope the
following code will have a better chance of surviving similar problems
in the future.  (The proposed code is a tad shorter and avoids some
code duplication and IF_LINT stuff; that's a good sign...)

2007-11-18  Paul Eggert  <address@hidden>

        * src/seq.c (print_numbers): Rewrite in an attempt to avoid the
        more-general rounding issues exposed by the previous patch.

diff --git a/src/seq.c b/src/seq.c
index eec5ed5..261a44b 100644
--- a/src/seq.c
+++ b/src/seq.c
@@ -238,24 +238,32 @@ static void
 print_numbers (char const *fmt, struct layout layout,
               long double first, long double step, long double last)
 {
-  long double i;
-  long double x0 IF_LINT (= 0);
+  bool out_of_range = (step < 0 ? first < last : last < first);

-  for (i = 0; /* empty */; i++)
+  if (! out_of_range)
     {
-      long double x = first + i * step;
+      long double x = first;
+      long double i;

-      if (step < 0 ? x < last : last < x)
+      for (i = 1; ; i++)
        {
-         /* If we go one past the end, but that number prints as a
-            value equal to "last", and prints differently from the
-            previous number, then print "last".  This avoids problems
-            with rounding.  For example, with the x86 it causes "seq
-            0 0.000001 0.000003" to print 0.000003 instead of
-            stopping at 0.000002.  */
-
-         if (i)
+         long double x0 = x;
+         printf (fmt, x);
+         if (out_of_range)
+           break;
+         x = first + i * step;
+         out_of_range = (step < 0 ? x < last : last < x);
+
+         if (out_of_range)
            {
+             /* If the number just past LAST prints as a value equal
+                to LAST, and prints differently from the previous
+                number, then print the number.  This avoids problems
+                with rounding.  For example, with the x86 it causes
+                "seq 0 0.000001 0.000003" to print 0.000003 instead
+                of stopping at 0.000002.  */
+
+             bool print_extra_number = false;
              long double x_val;
              char *x_str;
              int x_strlen = asprintf (&x_str, fmt, x);
@@ -269,34 +277,20 @@ print_numbers (char const *fmt, struct layout layout,
                  char *x0_str = NULL;
                  if (asprintf (&x0_str, fmt, x0) < 0)
                    xalloc_die ();
-                 if (!STREQ (x0_str, x_str))
-                   {
-                     fputs (separator, stdout);
-                     fputs (x_str, stdout);
-                   }
+                 print_extra_number = !STREQ (x0_str, x_str);
                  free (x0_str);
                }

              free (x_str);
+             if (! print_extra_number)
+               break;
            }

-         /* With floating point arithmetic, we may well reach this point
-            with i == 0 and first == last.  E.g., ./seq .1 .1 on FreeBSD 6.1.
-            Hence the first conjunct: don't break out of this loop when
-            i == 0.  *unless* first and last themselves are out of order,
-            in which case we must print nothing, e.g. for ./seq -1  */
-         if (i || (0 < step && last < first) || (step < 0 && first < last))
-           break;
+         fputs (separator, stdout);
        }

-      if (i)
-       fputs (separator, stdout);
-      printf (fmt, x);
-      x0 = x;
+      fputs (terminator, stdout);
     }
-
-  if (i)
-    fputs (terminator, stdout);
 }

 /* Return the default format given FIRST, STEP, and LAST.  */