--- a/liboctave/Range.cc	Fri Oct 06 06:14:45 1995 +0000
+++ b/liboctave/Range.cc	Sat Oct 07 03:39:26 1995 +0000
@@ -29,7 +29,9 @@
 #include <config.h>
 #endif
 
+#include <cfloat>
 #include <climits>
+#include <cmath>
 
 #include <iostream.h>
 
@@ -136,6 +138,81 @@
   return is;
 }
 
+// C  See Knuth, Art Of Computer Programming, Vol. 1, Problem 1.2.4-5.
+// C
+// C===Tolerant FLOOR function.
+// C
+// C    X  -  is given as a Double Precision argument to be operated on.
+// C          It is assumed that X is represented with M mantissa bits.
+// C    CT -  is   given   as   a   Comparison   Tolerance   such   that
+// C          0.LT.CT.LE.3-SQRT(5)/2. If the relative difference between
+// C          X and A whole number is  less  than  CT,  then  TFLOOR  is
+// C          returned   as   this   whole   number.   By  treating  the
+// C          floating-point numbers as a finite ordered set  note  that
+// C          the  heuristic  EPS=2.**(-(M-1))   and   CT=3*EPS   causes
+// C          arguments  of  TFLOOR/TCEIL to be treated as whole numbers
+// C          if they are  exactly  whole  numbers  or  are  immediately
+// C          adjacent to whole number representations.  Since EPS,  the
+// C          "distance"  between  floating-point  numbers  on  the unit
+// C          interval, and M, the number of bits in X'S mantissa, exist
+// C          on  every  floating-point   computer,   TFLOOR/TCEIL   are
+// C          consistently definable on every floating-point computer.
+// C
+// C          For more information see the following references:
+// C    (1) P. E. Hagerty, "More On Fuzzy Floor And Ceiling," APL  QUOTE
+// C        QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5.
+// C    (2) L. M. Breed, "Definitions For Fuzzy Floor And Ceiling",  APL
+// C        QUOTE QUAD 8(3):16-23, March 1978. This paper cites FL1 through
+// C        FL5, the history of five years of evolutionary development of
+// C        FL5 - the seven lines of code below - by open collaboration
+// C        and corroboration of the mathematical-computing community.
+// C
+// C  Penn State University Center for Academic Computing
+// C  H. D. Knoble - August, 1978.
+
+static inline double
+tfloor (double x, double ct)
+{
+// C---------FLOOR(X) is the largest integer algebraically less than
+// C         or equal to X; that is, the unfuzzy FLOOR function.
+
+//  DINT (X) = X - DMOD (X, 1.0);
+//  FLOOR (X) = DINT (X) - DMOD (2.0 + DSIGN (1.0, X), 3.0);
+
+// C---------Hagerty's FL5 function follows...
+
+  double q = 1.0;
+
+  if (x < 0.0)
+    q = 1.0 - ct;
+
+  double rmax = q / (2.0 - ct);
+
+  double t1 = 1.0 + floor (x);
+  t1 = (ct / q) * (t1 < 0.0 ? -t1 : t1);
+  t1 = rmax < t1 ? rmax : t1;
+  t1 = ct > t1 ? ct : t1;
+  t1 = floor (x + t1);
+
+  if (x <= 0.0 || t1 - x < rmax)
+    return t1;
+  else
+    return t1 - 1.0;
+}
+
+static inline double
+tceil (double x, double ct)
+{
+  return -tfloor (-x, ct);
+}
+
+static inline double
+round (double x, double ct)
+{
+  return tfloor (x+0.5, ct);
+}
+
+
 // Find an approximate number of intervals, then do the best we can to
 // find the number of intervals that we would get if we had done
 // something like
@@ -152,65 +229,13 @@
 int
 Range::nelem_internal (void) const
 {
-  // We can't have more than INT_MAX elements in the range.
-
-  double d_n_intervals = (rng_limit - rng_base) / rng_inc;
-  int max_intervals = INT_MAX - 1;
-  double d_max_val = (double) max_intervals;
-
-  if (d_n_intervals > d_max_val)
-    return -1;
+  double ct = 3.0 * DBL_EPSILON;
 
-  int n_intervals = (d_n_intervals > 0)
-    ? ((int) (d_n_intervals + 0.5))
-    : ((int) (d_n_intervals - 0.5)); 
-
-  if (rng_limit > rng_base && rng_inc > 0)
-    {
-      // Our approximation may have been too big.
-
-      while (rng_base + n_intervals * rng_inc > rng_limit && n_intervals > 0)
-	n_intervals--;
-
-      // Now that we are close, get the actual number.  Try to avoid
-      // problems with extended precision registers.
+  double tmp = round ((rng_limit - rng_base + rng_inc) / rng_inc, ct);
 
-      for (;;)
-	{
-	  volatile double tmp_inc = (n_intervals + 1) * rng_inc;
-	  volatile double tmp_val = rng_base + tmp_inc;
-	  if (tmp_val <= rng_limit && n_intervals < max_intervals)
-	    n_intervals++;
-	  else
-	    break;
-	}
-    }
-  else if (rng_limit < rng_base && rng_inc < 0)
-    {
-      // Our approximation may have been too big.
-
-      while (rng_base + n_intervals * rng_inc < rng_limit && n_intervals > 0)
-	n_intervals--;
+  int n_intervals = (int) (tmp > 0.0 ? tmp : 0);
 
-      // Now that we are close, get the actual number.  Try to avoid
-      // problems with extended precision registers.
-
-      for (;;)
-	{
-	  volatile double tmp_inc = (n_intervals + 1) * rng_inc;
-	  volatile double tmp_val = rng_base + tmp_inc;
-	  if (tmp_val >= rng_limit && n_intervals < max_intervals)
-	    n_intervals++;
-	  else
-	    break;
-	}
-    }
-  else if (rng_limit == rng_base)
-    n_intervals = 0;
-  else
-    n_intervals = -1;
-
-  return (n_intervals >= max_intervals) ? -1 : n_intervals + 1;
+  return (n_intervals >= INT_MAX - 1) ? -1 : n_intervals;
 }
 
 /*