# HG changeset patch
# User Daniel J Sebald <daniel.sebald@ieee.org>
# Date 1412154354 18000
# Node ID 53f156abe4c95e33aec601291aa92b9662c05abd
# Parent  1111d2d5ff95b64fcfa61bf209c60a8f4567ed95
New integerizing range algorithm for non-integer elements using better interp.

* configure.ac: Add nextafter library function to AC_CHECK_FUNCS list.

* liboctave/array/Range.h
(Range): Add variables lspc_fromA, lspc_toB, lspc_Ninc and lspc_scale for
linear spacing algorithm.  Add functions construct_linear_spacing(),
compute_int_elem(), compute_nonint_elem(), compute_elem() and
compute_elements().
(Range::compute_elem) Added.  Inline choice of algorithm based on whether
float parameters are integer.
(Range::Range): Call construct_linear_spacing() in all constructors.
(Range::Range(double,double)): Removed.  Use default input in third argument
of Range::Range(double,double,double).
(Range::set_base, Range::set_limit, Range::set_inc): Removed.  Disallow direct
access to params.

* liboctave/array/Range.cc
(Range::matrix_value): Remove variables b and increment.  Call
compute_elements() rather than use formula.
(Range::checkelem, Range::elem): Call compute_elem() rather than use formula.
(Range::operator <<): Replace variable b with a.base().  Replace variable
increment with a.inc().  Replace interpolation formula with a.compute_elem().
Call a.construct_linear_spacing() before returning.
(Range::operator -): Set linear spacing parameters while reversing polarity
of limits.
(Range::operator -, double and Range): Set linear spacing parameters while
subtracting limits from x divided by scale value.
(Range::operator -, Range and double): Set linear spacing parameters while
subtracting x divided by scale value from limits.
(Range::operator +): Set linear spacing parameters while adding x divided by
scale value.
(Range::operator *): Set linear spacing parameters while multiplying scale
value by x.
(Range::float_is_int): Added.  Declare a float value an integer if it is within
one floating point step from its rounded value.
(Range::construct_linear_spacing): Added.  Use float_is_int to determine which
normalized parameters are integers, then initialize linear spacing parameters
according to all the different combinations of limits and increment being
integer.
(Range::compute_int_elem) Added.  Existing algorithm.
(Range::compute_nonint_elem) Added.  New algorithm.
(Range::compute_elements) Added.  Looping routine utilizing compute_int_elem()
or compute_nonint_elem() accordingly.

diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -2181,7 +2181,7 @@
 AC_CHECK_FUNCS([getpgrp getpid getppid getpwent getpwuid getuid])
 AC_CHECK_FUNCS([isascii kill])
 AC_CHECK_FUNCS([lgamma lgammaf lgamma_r lgammaf_r])
-AC_CHECK_FUNCS([log1p log1pf pipe])
+AC_CHECK_FUNCS([log1p log1pf nextafter pipe])
 AC_CHECK_FUNCS([realpath resolvepath roundl])
 AC_CHECK_FUNCS([select setgrent setpwent siglongjmp strsignal])
 AC_CHECK_FUNCS([tcgetattr tcsetattr tgammaf toascii])
diff --git a/liboctave/array/Range.cc b/liboctave/array/Range.cc
--- a/liboctave/array/Range.cc
+++ b/liboctave/array/Range.cc
@@ -54,15 +54,12 @@
   if (rng_nelem > 0 && cache.nelem () == 0)
     {
       cache.resize (1, rng_nelem);
-      double b = rng_base;
-      double increment = rng_inc;
       if (rng_nelem > 0)
         {
           // The first element must always be *exactly* the base.
           // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment).
-          cache(0) = b;
-          for (octave_idx_type i = 1; i < rng_nelem; i++)
-            cache(i) = b + i * increment;
+          cache(0) = rng_base;
+          compute_elements (cache, 1, rng_nelem - 1);
         }
 
       // On some machines (x86 with extended precision floating point
@@ -89,7 +86,7 @@
   if (i == 0)
     return rng_base;
   else if (i < rng_nelem - 1)
-    return rng_base + i * rng_inc;
+    return compute_elem (i);
   else
     {
       double end = rng_base + i * rng_inc;
@@ -110,7 +107,7 @@
   if (i == 0)
     return rng_base;
   else if (i < rng_nelem - 1)
-    return rng_base + i * rng_inc;
+    return compute_elem (i);
   else
     {
       double end = rng_base + i * rng_inc;
@@ -135,6 +132,7 @@
     if (i == 0)
       *array++ = base;
     else if (i < nmax)
+      // Indeces use the all-integer algorithm.
       *array++ = base + i * inc;
     else
       {
@@ -349,20 +347,18 @@
 std::ostream&
 operator << (std::ostream& os, const Range& a)
 {
-  double b = a.base ();
-  double increment = a.inc ();
   octave_idx_type num_elem = a.nelem ();
 
   if (num_elem > 1)
     {
       // First element must be the base *exactly* (-0).
-      os << b << " ";
+      os << a.base () << " ";
       for (octave_idx_type i = 1; i < num_elem-1; i++)
-        os << b + i * increment << " ";
+        os << a.compute_elem (i) << " ";
     }
 
   // Prevent overshoot.  See comment in the matrix_value method above.
-  os << (increment > 0 ? a.max () : a.min ()) << "\n";
+  os << (a.inc () > 0 ? a.max () : a.min ()) << "\n";
 
   return os;
 }
@@ -381,18 +377,30 @@
         }
     }
 
+  a.construct_linear_spacing ();
+
   return is;
 }
 
 Range
 operator - (const Range& r)
 {
-  return Range (-r.base (), -r.limit (), -r.inc (), r.nelem ());
+  Range result (-r.base (), -r.limit (), -r.inc (), r.nelem ());
+  result.lspc_fromA = -r.lspc_fromA ;
+  result.lspc_toB = -r.lspc_toB ;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = r.lspc_scale;
+
+  return result;
 }
 
 Range operator + (double x, const Range& r)
 {
   Range result (x + r.base (), x + r.limit (), r.inc (), r.nelem ());
+  result.lspc_fromA = x/r.lspc_scale + r.lspc_fromA ;
+  result.lspc_toB = x/r.lspc_scale + r.lspc_toB ;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = r.lspc_scale;
   if (result.rng_nelem < 0)
     result.cache = x + r.matrix_value ();
 
@@ -402,6 +410,10 @@
 Range operator + (const Range& r, double x)
 {
   Range result (r.base () + x, r.limit () + x, r.inc (), r.nelem ());
+  result.lspc_fromA = r.lspc_fromA + x/r.lspc_scale;
+  result.lspc_toB = r.lspc_toB + x/r.lspc_scale;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = r.lspc_scale;
   if (result.rng_nelem < 0)
     result.cache = r.matrix_value () + x;
 
@@ -411,6 +423,10 @@
 Range operator - (double x, const Range& r)
 {
   Range result (x - r.base (), x - r.limit (), -r.inc (), r.nelem ());
+  result.lspc_fromA = x/r.lspc_scale - r.lspc_fromA ;
+  result.lspc_toB = x/r.lspc_scale - r.lspc_toB ;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = r.lspc_scale;
   if (result.rng_nelem < 0)
     result.cache = x - r.matrix_value ();
 
@@ -420,6 +436,10 @@
 Range operator - (const Range& r, double x)
 {
   Range result (r.base () - x, r.limit () - x, r.inc (), r.nelem ());
+  result.lspc_fromA = r.lspc_fromA - x/r.lspc_scale;
+  result.lspc_toB = r.lspc_toB - x/r.lspc_scale;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = r.lspc_scale;
   if (result.rng_nelem < 0)
     result.cache = r.matrix_value () - x;
 
@@ -429,6 +449,10 @@
 Range operator * (double x, const Range& r)
 {
   Range result (x * r.base (), x * r.limit (), x * r.inc (), r.nelem ());
+  result.lspc_fromA = r.lspc_fromA;
+  result.lspc_toB = r.lspc_toB;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = x * r.lspc_scale;
   if (result.rng_nelem < 0)
     result.cache = x * r.matrix_value ();
 
@@ -438,6 +462,10 @@
 Range operator * (const Range& r, double x)
 {
   Range result (r.base () * x, r.limit () * x, r.inc () * x, r.nelem ());
+  result.lspc_fromA = r.lspc_fromA;
+  result.lspc_toB = r.lspc_toB;
+  result.lspc_Ninc = r.lspc_Ninc;
+  result.lspc_scale = r.lspc_scale * x;
   if (result.rng_nelem < 0)
     result.cache = r.matrix_value () * x;
 
@@ -543,6 +571,9 @@
       octave_idx_type n_elt = (tmp > 0.0 ? static_cast<octave_idx_type> (tmp)
                                          : 0);
 
+      // Perhaps the routine float_is_int (double x) located in this
+      // file will help here.
+
       // If the final element that we would compute for the range is
       // equal to the limit of the range, or is an adjacent floating
       // point number, accept it.  Otherwise, try a range with one
@@ -568,3 +599,137 @@
 
   return retval;
 }
+
+bool
+float_is_int (double x)
+{
+  // https://www.gnu.org/software/libc/manual/html_node/Rounding-Functions.html
+  // https://www.gnu.org/software/libc/manual/html_node/FP-Bit-Twiddling.html
+  double absx = fabs (x);
+  double roundabsx = gnulib::round (absx);
+  return (nextafter (absx, roundabsx) == roundabsx);
+}
+
+void
+Range::construct_linear_spacing (void)
+{
+#define SHOW_RANGE_CLASS 0
+  // Attempt to integerize the range as much as possible.
+  double A_norm = rng_base/rng_inc;
+  double B_norm = rng_limit/rng_inc;
+  double N_norm = (rng_limit - rng_base)/rng_inc;
+  bool A_is_int = float_is_int (A_norm);
+  bool B_is_int = float_is_int (B_norm);
+  bool N_is_int = float_is_int (N_norm);
+  if (A_is_int && B_is_int)
+    {
+#if SHOW_RANGE_CLASS
+      fprintf(stderr, "A_is_int && B_is_int\n");
+#endif
+      // Can integerize, Ninc derived.
+      lspc_fromA = gnulib::round (A_norm);
+      lspc_toB = gnulib::round (B_norm);
+      lspc_scale = rng_inc;
+      lspc_Ninc = abs(lspc_fromA - lspc_toB);
+    }
+  else if (A_is_int && N_is_int)
+    {
+#if SHOW_RANGE_CLASS
+      fprintf(stderr, "A_is_int && N_is_int\n");
+#endif
+      // Can integerize, Ninc accepted.
+      lspc_fromA = gnulib::round (A_norm);
+      lspc_toB = B_norm;
+      lspc_scale = rng_inc;
+      lspc_Ninc = abs (gnulib::round (N_norm));
+    }
+  else if (A_is_int)
+    {
+#if SHOW_RANGE_CLASS
+      fprintf(stderr, "A_is_int\n");
+#endif
+      // Can integerize, but Ninc and toB derived.
+      lspc_fromA = gnulib::round (A_norm);
+      lspc_scale = rng_inc;
+      lspc_Ninc = gnulib::floor (abs (N_norm));
+      lspc_toB = lspc_fromA + lspc_Ninc;
+    }
+  else if (N_is_int)
+    {
+#if SHOW_RANGE_CLASS
+      fprintf(stderr, "N_is_int\n");
+#endif
+      // Cannot integerize, but limits accepted.
+      lspc_fromA = rng_base;
+      lspc_toB = rng_limit;
+      lspc_scale = 1.0;
+      lspc_Ninc = abs (gnulib::round (N_norm));
+    }
+  else
+    {
+      // Cannot integerize, and Ninc and toB derived.
+      lspc_fromA = rng_base;
+      lspc_scale = 1.0;
+      lspc_Ninc = gnulib::floor (abs (N_norm));
+      lspc_toB = lspc_fromA + lspc_Ninc;
+    }
+#undef SHOW_RANGE_CLASS
+}
+
+// When satisfiend with this result, move to header file as an inline function
+// for speed.  (Touching header file propogates through source tree.)
+double
+Range::compute_int_elem (octave_idx_type i) const
+{
+  // This algorithm is exact when base and increment are integers
+  // but skews the error residuals toward the end of the range when
+  // non integers.
+  return rng_base + i * rng_inc;
+}
+
+// When satisfiend with this result, move to header file as an inline function
+// for speed.  (Touching header file propogates through source tree.)
+double
+Range::compute_nonint_elem (octave_idx_type i) const
+{
+  // These are variations on the general OP's algorithm, which seems to
+  // be preferred by Internet dicussion boards.
+
+  // This algorithm has more flat distribution of residuals, but generally
+  // a higher residual norm which makes it fail a sine crossing test.
+//  return ( (lspc_fromA * (1.0 - (i/lspc_Ninc)) + lspc_toB * (i/lspc_Ninc)) * lspc_scale );
+
+  // This algorithm symmetrically distributes error residuals throughout
+  // the range and consistently creates one of the lower residual norms.
+  return ( ( (lspc_fromA * (lspc_Ninc - i) + lspc_toB * i) / lspc_Ninc ) * lspc_scale);
+
+  // There is a linear interpolation in the GNU library, but it is a bit
+  // of work to program, so we'll leave it as a future investigation.
+  // https://www.gnu.org/software/gsl/manual/html_node/Interpolation-Types.html
+// gsl_interp_linear somehow
+}
+
+// When satisfiend with this result, move to header file as an inline function
+// for speed.  (Touching header file propogates through source tree.)
+void
+Range::compute_elements (Matrix &outvec, octave_idx_type start, octave_idx_type end) const
+{
+#define SHOW_RANGE_ALGORITHM 0
+  if (all_elements_are_ints ())
+    {
+#if SHOW_RANGE_ALGORITHM
+      fprintf(stderr,"compute_int_elements, rng_base=%lg, rng_inc=%lg\n", rng_base, rng_inc);
+#endif
+      for (octave_idx_type i = start; i <= end; i++)
+        outvec(i) = compute_int_elem (i);
+    }
+  else
+    {
+#if SHOW_RANGE_ALGORITHM
+      fprintf(stderr,"compute_nonint_elements, lspc_Ninc=%lg, lspc_scale=%lg\n", lspc_Ninc, lspc_scale);
+#endif
+      for (octave_idx_type i = start; i <= end; i++)
+        outvec(i) = compute_nonint_elem (i);
+    }
+#undef SHOW_RANGE_ALGORITHM
+}
diff --git a/liboctave/array/Range.h b/liboctave/array/Range.h
--- a/liboctave/array/Range.h
+++ b/liboctave/array/Range.h
@@ -35,19 +35,24 @@
 public:
 
   Range (void)
-    : rng_base (0), rng_limit (0), rng_inc (0), rng_nelem (0), cache (1, 0) { }
+    : rng_base (0), rng_limit (0), rng_inc (0), rng_nelem (0), cache (1, 0)
+  {
+    construct_linear_spacing ();
+  }
 
   Range (const Range& r)
     : rng_base (r.rng_base), rng_limit (r.rng_limit), rng_inc (r.rng_inc),
-      rng_nelem (r.rng_nelem), cache (r.cache) { }
+      rng_nelem (r.rng_nelem), cache (r.cache)
+  {
+    construct_linear_spacing ();
+  }
 
-  Range (double b, double l)
-    : rng_base (b), rng_limit (l), rng_inc (1),
-      rng_nelem (nelem_internal ()), cache () { }
-
-  Range (double b, double l, double i)
+  Range (double b, double l, double i = 1)
     : rng_base (b), rng_limit (l), rng_inc (i),
-      rng_nelem (nelem_internal ()), cache () { }
+      rng_nelem (nelem_internal ()), cache ()
+  {
+    construct_linear_spacing ();
+  }
 
   // For operators' usage (to preserve element count).
   Range (double b, double i, octave_idx_type n)
@@ -56,6 +61,8 @@
   {
     if (! xfinite (b) || ! xfinite (i) | ! xfinite (rng_limit))
       rng_nelem = -2;
+
+    construct_linear_spacing ();
   }
 
   double base (void) const { return rng_base; }
@@ -90,33 +97,6 @@
 
   Array<double> index (const idx_vector& i) const;
 
-  void set_base (double b)
-  {
-    if (rng_base != b)
-      {
-        rng_base = b;
-        clear_cache ();
-      }
-  }
-
-  void set_limit (double l)
-  {
-    if (rng_limit != l)
-      {
-        rng_limit = l;
-        clear_cache ();
-      }
-  }
-
-  void set_inc (double i)
-  {
-    if (rng_inc != i)
-      {
-        rng_inc = i;
-        clear_cache ();
-      }
-  }
-
   friend OCTAVE_API std::ostream& operator << (std::ostream& os,
                                                const Range& r);
   friend OCTAVE_API std::istream& operator >> (std::istream& is, Range& r);
@@ -133,11 +113,17 @@
 
 private:
 
+  // User inputs
   double rng_base;
   double rng_limit;
   double rng_inc;
 
+  // Derived parameters
   octave_idx_type rng_nelem;
+  double lspc_fromA;
+  double lspc_toB;
+  double lspc_Ninc;
+  double lspc_scale;
 
   mutable Matrix cache;
 
@@ -145,6 +131,19 @@
 
   void clear_cache (void) const { cache.resize (0, 0); }
 
+  void construct_linear_spacing (void);
+
+  double compute_int_elem (octave_idx_type i) const;
+  double compute_nonint_elem (octave_idx_type i) const;
+  inline double compute_elem (octave_idx_type i) const
+    {
+      if (all_elements_are_ints ())
+        return compute_int_elem (i);
+      else
+        return compute_nonint_elem (i);
+    }
+  void compute_elements (Matrix &output, octave_idx_type start, octave_idx_type end) const;
+
 protected:
 
   // For operators' usage (to allow all values to be set directly).
@@ -154,6 +153,8 @@
   {
     if (! xfinite (b) || ! xfinite (i) || ! xfinite (l))
       rng_nelem = -2;
+
+    construct_linear_spacing ();
   }
 };
 
exporting patch:
<fdopen>
