--- a/liboctave/array/Range.cc	Tue Aug 04 10:56:14 2020 -0400
+++ b/liboctave/array/Range.cc	Thu Aug 06 13:53:07 2020 -0400
@@ -39,6 +39,225 @@
 #include "lo-mappers.h"
 #include "lo-utils.h"
 
+namespace octave
+{
+  template <typename T>
+  T xtfloor (T x, T 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...
+
+    T q = 1;
+
+    if (x < 0)
+      q = 1 - ct;
+
+    T rmax = q / (2 - ct);
+
+    T t1 = 1 + std::floor (x);
+    t1 = (ct / q) * (t1 < 0 ? -t1 : t1);
+    t1 = (rmax < t1 ? rmax : t1);
+    t1 = (ct > t1 ? ct : t1);
+    t1 = std::floor (x + t1);
+
+    if (x <= 0 || (t1 - x) < rmax)
+      return t1;
+    else
+      return t1 - 1;
+  }
+
+  template <typename T>
+  bool xteq (T u, T v, T ct = 3 * std::numeric_limits<T>::epsilon ())
+  {
+    T tu = std::abs (u);
+    T tv = std::abs (v);
+
+    return std::abs (u - v) < ((tu > tv ? tu : tv) * ct);
+  }
+
+  template <typename T>
+  octave_idx_type xnumel_internal (T base, T limit, T inc)
+  {
+    octave_idx_type retval = -1;
+    if (! octave::math::isfinite (base) || ! octave::math::isfinite (inc)
+        || octave::math::isnan (limit))
+      retval = -2;
+    else if (octave::math::isinf (limit)
+             && ((inc > 0 && limit > 0)
+                 || (inc < 0 && limit < 0)))
+      retval = std::numeric_limits<octave_idx_type>::max () - 1;
+    else if (inc == 0
+             || (limit > base && inc < 0)
+             || (limit < base && inc > 0))
+      {
+        retval = 0;
+      }
+    else
+      {
+        T ct = 3 * std::numeric_limits<T>::epsilon ();
+
+        T tmp = xtfloor ((limit - base + inc) / inc, ct);
+
+        octave_idx_type n_elt
+          = (tmp > 0 ? static_cast<octave_idx_type> (tmp) : 0);
+
+        // 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
+        // fewer element.  If that fails, try again with one more
+        // element.
+        //
+        // I'm not sure this is very good, but it seems to work better
+        // than just using tfloor as above.  For example, without it,
+        // the expression 1.8:0.05:1.9 fails to produce the expected
+        // result of [1.8, 1.85, 1.9].
+
+        if (! xteq (base + (n_elt - 1) * inc, limit))
+          {
+            if (xteq (base + (n_elt - 2) * inc, limit))
+              n_elt--;
+            else if (xteq (base + n_elt * inc, limit))
+              n_elt++;
+          }
+
+        retval = (n_elt < std::numeric_limits<octave_idx_type>::max () - 1
+                  ? n_elt : -1);
+      }
+
+    return retval;
+  }
+
+  template <typename T>
+  bool xall_elements_are_ints (T base, T inc, octave_idx_type nel)
+  {
+    // If the base and increment are ints, the final value in the range
+    // will also be an integer, even if the limit is not.  If the range
+    // has only one or zero elements, then the base needs to be an integer.
+
+    return (! (octave::math::isnan (base) || octave::math::isnan (inc))
+            && (octave::math::nint_big (base) == base || nel < 1)
+            && (octave::math::nint_big (inc) == inc || nel <= 1));
+  }
+
+  template <typename T>
+  T
+  xfinal_value (T base, T limit, T inc, octave_idx_type nel)
+  {
+    T retval = T (0);
+
+    if (nel <= 1)
+      return base;
+
+    // If increment is 0, then numel should also be zero.
+
+    retval = base + (nel - 1) * inc;
+
+    // On some machines (x86 with extended precision floating point
+    // arithmetic, for example) it is possible that we can overshoot
+    // the limit by approximately the machine precision even though
+    // we were very careful in our calculation of the number of
+    // elements.  Therefore, we clip the result to the limit if it
+    // overshoots.
+
+    // NOTE: The test also includes equality (>= limit) to have
+    // expressions such as -5:1:-0 result in a -0 endpoint.
+
+    if ((inc > T (0) && retval >= limit) || (inc < T (0) && retval <= limit))
+      retval = limit;
+
+    // If all elements are integers, then ensure the final value is.
+
+    if (xall_elements_are_ints (base, inc, nel))
+      retval = std::round (retval);
+
+    return retval;
+  }
+
+  template <>
+  bool
+  range<double>::all_elements_are_ints (void) const
+  {
+    return xall_elements_are_ints (m_base, m_increment, m_numel);
+  }
+
+  template <>
+  bool
+  range<float>::all_elements_are_ints (void) const
+  {
+    return xall_elements_are_ints (m_base, m_increment, m_numel);
+  }
+
+  template <>
+  octave_idx_type
+  range<double>::get_numel (void) const
+  {
+    return xnumel_internal (m_base, m_limit, m_increment);
+  }
+
+  template <>
+  octave_idx_type
+  range<float>::get_numel (void) const
+  {
+    return xnumel_internal (m_base, m_limit, m_increment);
+  }
+
+  template <>
+  double
+  range<double>::get_final_value (void) const
+  {
+    return xfinal_value (m_base, m_limit, m_increment, m_numel);
+  }
+
+  template <>
+  float
+  range<float>::get_final_value (void) const
+  {
+    return xfinal_value (m_base, m_limit, m_increment, m_numel);
+  }
+
+  template <>
+  octave_idx_type
+  range<double>::nnz (void) const
+  {
+    octave_idx_type retval = 0;
+
+    if (! isempty ())
+      {
+        if ((m_base > 0 && m_limit > 0)
+            || (m_base < 0 && m_limit < 0))
+          {
+            // All elements have the same sign, hence there are no zeros.
+            retval = m_numel;
+          }
+        else if (m_increment != 0)
+          {
+            if (m_base == 0 || m_limit == 0)
+              // Exactly one zero at beginning or end of range.
+              retval = m_numel - 1;
+            else if (octave::math::mod (-m_base, m_increment) != 0)
+              // Range crosses negative/positive without hitting zero.
+              retval = m_numel;
+            else
+              // Range crosses negative/positive and hits zero.
+              retval = m_numel - 1;
+          }
+        else
+          {
+            // All elements are equal (m_increment = 0) but not
+            // positive or negative, therefore all elements are zero.
+            retval = 0;
+          }
+      }
+
+    return retval;
+  }
+}
+
 bool
 Range::all_elements_are_ints (void) const
 {

--- a/liboctave/array/Range.h	Tue Aug 04 10:56:14 2020 -0400
+++ b/liboctave/array/Range.h	Thu Aug 06 13:53:07 2020 -0400
@@ -32,8 +32,323 @@
 
 #include "dMatrix.h"
 #include "dim-vector.h"
+#include "lo-error.h"
 #include "oct-sort.h"
 
+template <typename T> class Array;
+
+namespace octave
+{
+  // Helper class used solely for idx_vector.loop () function call
+
+  template <typename T>
+  class rangeidx_helper
+  {
+  public:
+
+    rangeidx_helper (T *a, T b, T i, T l, octave_idx_type n)
+      : array (a), base (b), inc (i), limit (l), nmax (n-1) { }
+
+    void operator () (octave_idx_type i)
+    {
+      if (i == 0)
+        *array++ = base;
+      else if (i < nmax)
+        *array++ = base + T (i) * inc;
+      else
+        *array++ = limit;
+    }
+
+  private:
+
+    T *array, base, inc, limit;
+    octave_idx_type nmax;
+
+  };
+
+  template <typename T>
+  class range
+  {
+  public:
+
+    range (void)
+      : m_base (0), m_increment (0), m_limit (0), m_final (0), m_numel (0)
+    { }
+
+    // LIMIT is an upper limit and may be outside the range of actual
+    // values.  For floating point ranges, we perform a tolerant check
+    // to attempt to capture limit in the set of values if it is "close"
+    // to the value of base + a multiple of the increment.
+
+    range (const T& base, const T& increment, const T& limit)
+      : m_base (base), m_increment (increment), m_limit (limit),
+        m_final (), m_numel ()
+    {
+      init ();
+    }
+
+    range (const T& base, const T& limit)
+      : m_base (base), m_increment (1), m_limit (limit), m_final (), m_numel ()
+    {
+      init ();
+    }
+
+    // Allow conversion from (presumably) properly constructed Range
+    // objects and to create constant ranges (see the static
+    // make_constant method).  The values of base, limit, increment,
+    // and numel must be consistent.
+
+    // FIXME: Actually check that base, limit, increment, and numel are
+    // consistent.
+
+    // FIXME: Is there a way to limit this to T == double?
+
+    range (const T& base, const T& increment, const T& limit,
+           octave_idx_type numel)
+      : m_base (base), m_increment (increment), m_limit (limit),
+        m_final (limit), m_numel (numel)
+    { }
+
+    // We don't use a constructor for this because it will conflict with
+    // range<T> (base, limit) when T is octave_idx_type.
+
+    static range<T> make_constant (const T& base, octave_idx_type numel)
+    {
+      // We could just make this constructor public, but it allows
+      // inconsistent ranges to be constructed.  And it is probably much
+      // clearer to see "make_constant" instead of puzzling over the
+      // purpose of this strange constructor form.
+
+      return range<T> (base, T (), base, numel);
+    }
+
+    // We don't use a constructor for this because it will conflict with
+    // range<T> (base, limit, increment) when T is octave_idx_type.
+
+    static range<T> make_n_element_range (const T& base, const T& increment,
+                                          octave_idx_type numel)
+    {
+      // We could just make this constructor public, but it allows
+      // inconsistent ranges to be constructed.  And it is probably much
+      // clearer to see "make_constant" instead of puzzling over the
+      // purpose of this strange constructor form.
+
+      T final_val = base + double (numel - 1) * increment;
+
+      return range<T> (base, increment, final_val, numel);
+    }
+
+    range (const range<T>&) = default;
+
+    range<T>& operator = (const range<T>&) = default;
+
+    ~range (void) = default;
+
+    T base (void) const { return m_base; }
+    T increment (void) const { return m_increment; }
+    T limit (void) const { return m_limit; }
+
+    T final_value (void) const { return m_final; }
+
+    T min (void) const
+    {
+      return (m_numel > 0
+              ? m_increment > T (0) ? base () : final_value ()
+              : T (0));
+    }
+
+    T max (void) const
+    {
+      return (m_numel > 0
+              ? m_increment > T (0) ? final_value () : base ()
+              : T (0));
+    }
+
+    octave_idx_type numel (void) const { return m_numel; }
+
+    dim_vector dims (void) const { return dim_vector (1, m_numel); }
+
+    octave_idx_type rows (void) const { return 1; }
+
+    octave_idx_type cols (void) const { return numel (); }
+    octave_idx_type columns (void) const { return numel (); }
+
+    bool isempty (void) const { return numel () == 0; }
+
+    bool all_elements_are_ints (void) const { return true; }
+
+    sortmode issorted (sortmode mode = ASCENDING) const
+    {
+      if (m_numel > 1 && m_increment > T (0))
+        mode = (mode == DESCENDING) ? UNSORTED : ASCENDING;
+      else if (m_numel > 1 && m_increment < T (0))
+        mode = (mode == ASCENDING) ? UNSORTED : DESCENDING;
+      else
+        mode = (mode == UNSORTED) ? ASCENDING : mode;
+
+      return mode;
+    }
+
+    octave_idx_type nnz (void) const;
+
+    // Support for single-index subscripting, without generating matrix cache.
+
+    T checkelem (octave_idx_type i) const
+    {
+      if (i < 0 || i >= m_numel)
+        octave::err_index_out_of_range (2, 2, i+1, m_numel, dims ());
+
+      if (i == 0)
+        return m_base;
+      else if (i < m_numel - 1)
+        return m_base + T (i) * m_increment;
+      else
+        return final_value ();
+    }
+
+    T checkelem (octave_idx_type i, octave_idx_type j) const
+    {
+      // Ranges are *always* row vectors.
+      if (i != 0)
+        octave::err_index_out_of_range (1, 1, i+1, m_numel, dims ());
+
+      return checkelem (j);
+    }
+
+    T elem (octave_idx_type i) const
+    {
+      if (i == 0)
+        return m_base;
+      else if (i < m_numel - 1)
+        return m_base + T (i) * m_increment;
+      else
+        return final_value ();
+    }
+
+    T elem (octave_idx_type /* i */, octave_idx_type j) const
+    {
+      return elem (j);
+    }
+
+    T operator () (octave_idx_type i) const
+    {
+      return elem (i);
+    }
+
+    T operator () (octave_idx_type i, octave_idx_type j) const
+    {
+      return elem (i, j);
+    }
+
+    Array<T> index (const idx_vector& i) const
+    {
+      Array<T> retval;
+
+      octave_idx_type n = m_numel;
+
+      if (i.is_colon ())
+        {
+          retval = array_value ().reshape (dim_vector (m_numel, 1));
+        }
+      else
+        {
+          if (i.extent (n) != n)
+            octave::err_index_out_of_range (1, 1, i.extent (n), n, dims ());
+
+          dim_vector rd = i.orig_dimensions ();
+          octave_idx_type il = i.length (n);
+
+          // taken from Array.cc.
+          if (n != 1 && rd.isvector ())
+            rd = dim_vector (1, il);
+
+          retval.clear (rd);
+
+          // idx_vector loop across all values in i,
+          // executing __rangeidx_helper (i) for each i
+          i.loop (n, rangeidx_helper<T> (retval.fortran_vec (),
+                                         m_base, m_increment, final_value (),
+                                         m_numel));
+        }
+
+      return retval;
+    }
+
+    Array<T> diag (octave_idx_type k) const
+    {
+      return array_value ().diag (k);
+    }
+
+    Array<T> array_value (void) const
+    {
+      octave_idx_type nel = numel ();
+
+      Array<T> retval (dim_vector (1, nel));
+
+      if (nel > 0)
+        {
+          // The first element must always be *exactly* the base.
+          // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment).
+          retval(0) = m_base;
+
+          for (octave_idx_type i = 1; i < nel - 1; i++)
+            retval.xelem (i) = m_base + double (i) * m_increment;
+
+          retval.xelem (nel - 1) = final_value ();
+        }
+
+      return retval;
+    }
+
+  private:
+
+    T m_base;
+    T m_increment;
+    T m_limit;
+    T m_final;
+    octave_idx_type m_numel;
+
+    void init (void)
+    {
+      m_numel = get_numel ();
+      m_final = get_final_value ();
+    }
+
+    // Setting the number of elements to zero when the increment is zero
+    // is intentional and matches the behavior of Matlab's colon
+    // operator.
+
+    octave_idx_type get_numel (void) const
+    {
+      return ((m_increment == T (0)
+               || (m_limit > m_base && m_increment < T (0))
+               || (m_limit < m_base && m_increment > T (0)))
+              ? T (0)
+              : (m_limit - m_base + m_increment) / m_increment);
+    }
+
+    // This calculation is appropriate for integer ranges.
+
+    T get_final_value (void) const
+    {
+      return m_base + double (m_numel - 1) * m_increment;
+    }
+  };
+
+  // Specializations defined externally.
+
+  template <> bool range<double>::all_elements_are_ints (void) const;
+  template <> bool range<float>::all_elements_are_ints (void) const;
+
+  template <> octave_idx_type range<double>::get_numel (void) const;
+  template <> octave_idx_type range<float>::get_numel (void) const;
+
+  template <> double range<double>::get_final_value (void) const;
+  template <> float range<float>::get_final_value (void) const;
+
+  template <> octave_idx_type range<double>::nnz (void) const;
+}
+
 class
 OCTAVE_API
 Range
@@ -44,6 +359,16 @@
     : m_base (0), m_limit (0), m_inc (0), m_numel (0)
   { }
 
+  // Assume range is already properly constructed, so just copy internal
+  // values.  However, we set LIMIT to the computed final value because
+  // that mimics the behavior of the other Range class constructors that
+  // reset limit to the computed final value.
+
+  Range (const octave::range<double>& r)
+    : m_base (r.base ()), m_limit (r.final_value ()), m_inc (r.increment ()),
+      m_numel (r.numel ())
+  { }
+
   Range (const Range& r) = default;
 
   Range& operator = (const Range& r) = default;
@@ -67,7 +392,8 @@
   // NOTE: The following constructor may be deprecated and removed after
   // the arithmetic operators are removed.
 
-  // For operators' usage (to preserve element count).
+  // For operators' usage (to preserve element count) and to create
+  // constant row vectors (obsolete usage).
 
   Range (double b, double i, octave_idx_type n)
     : m_base (b), m_limit (b + (n-1) * i), m_inc (i), m_numel (n)