--- a/libinterp/interpfcn/pr-output.cc
+++ b/libinterp/interpfcn/pr-output.cc
@@ -2724,14 +2724,17 @@
                 {
                   octave_quit ();
 
-                  double val = base + i * increment;
+                  double val;
+                  if (i == 0)
+                    val = base;
+                  else
+                    val = base + i * increment;
 
                   if (i == num_elem - 1)
                     {
                       // See the comments in Range::matrix_value.
-
-                      if ((increment > 0 && val > limit)
-                          || (increment < 0 && val < limit))
+                      if ((increment > 0 && val >= limit)
+                          || (increment < 0 && val <= limit))
                         val = limit;
                     }
 

--- a/liboctave/array/Range.cc
+++ b/liboctave/array/Range.cc
@@ -36,14 +36,6 @@
 #include "lo-utils.h"
 #include "Array-util.h"
 
-Range::Range (double b, double i, octave_idx_type n)
-  : rng_base (b), rng_limit (b + (n-1) * i), rng_inc (i),
-    rng_nelem (n), cache ()
-{
-  if (! xfinite (b) || ! xfinite (i))
-    rng_nelem = -2;
-}
-
 bool
 Range::all_elements_are_ints (void) const
 {
@@ -64,17 +56,24 @@
       cache.resize (1, rng_nelem);
       double b = rng_base;
       double increment = rng_inc;
-      for (octave_idx_type i = 0; i < rng_nelem; i++)
-        cache(i) = b + i * increment;
+      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;
+        }
 
       // 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.
+      // elements.  The tests need equality (>= rng_limit or <= rng_limit)
+      // to have expressions like -5:1:-0 result in a -0 endpoint.
 
-      if ((rng_inc > 0 && cache(rng_nelem-1) > rng_limit)
-          || (rng_inc < 0 && cache(rng_nelem-1) < rng_limit))
+      if ((rng_inc > 0 && cache(rng_nelem-1) >= rng_limit)
+          || (rng_inc < 0 && cache(rng_nelem-1) <= rng_limit))
         cache(rng_nelem-1) = rng_limit;
     }
 
@@ -87,16 +86,68 @@
   if (i < 0 || i >= rng_nelem)
     gripe_index_out_of_range (1, 1, i+1, rng_nelem);
 
-  return rng_base + rng_inc * i;
+  if (i == 0)
+    return rng_base;
+  else if (i < rng_nelem - 1) 
+    return rng_base + i * rng_inc;
+  else
+    {
+      double end = rng_base + i * rng_inc;
+      if ((rng_inc > 0 && end >= rng_limit)
+          || (rng_inc < 0 && end <= rng_limit))
+        return rng_limit;
+      else
+        return end;
+    }
 }
 
+double
+Range::elem (octave_idx_type i) const
+{
+#if defined (BOUNDS_CHECKING)
+  return checkelem (i);
+#else
+  if (i == 0)
+    return rng_base;
+  else if (i < rng_nelem - 1) 
+    return rng_base + i * rng_inc;
+  else
+    {
+      double end = rng_base + i * rng_inc;
+      if ((rng_inc > 0 && end >= rng_limit)
+          || (rng_inc < 0 && end <= rng_limit))
+        return rng_limit;
+      else
+        return end;
+    }
+#endif
+}
+
+// Pseudo-class used for idx_vector.loop () function call
 struct _rangeidx_helper
 {
-  double *array, base, inc;
-  _rangeidx_helper (double *a, double b, double i)
-    : array (a), base (b), inc (i) { }
+  double *array, base, inc, limit;
+  octave_idx_type nmax;
+
+  _rangeidx_helper (double *a, double b, double i, double l, octave_idx_type n)
+    : array (a), base (b), inc (i), limit (l), nmax (n-1) { }
+
   void operator () (octave_idx_type i)
-    { *array++ = base + i * inc; }
+    {
+      if (i == 0)
+        *array++ = base;
+      else if (i < nmax) 
+        *array++ = base + i * inc;
+      else
+        {
+          double end = base + i * inc;
+          if ((inc > 0 && end >= limit)
+              || (inc < 0 && end <= limit))
+            *array++ = limit;
+          else
+            *array++ = end;
+        }
+    }
 };
 
 Array<double>
@@ -119,13 +170,14 @@
       octave_idx_type il = i.length (n);
 
       // taken from Array.cc.
-
       if (n != 1 && rd.is_vector ())
         rd = dim_vector (1, il);
 
       retval.clear (rd);
 
-      i.loop (n, _rangeidx_helper (retval.fortran_vec (), rng_base, rng_inc));
+      // idx_vector loop across all values in i, executing _rangeidx_helper (i) foreach i
+      i.loop (n, _rangeidx_helper (retval.fortran_vec (),
+                                   rng_base, rng_inc, rng_limit, rng_nelem));
     }
 
   return retval;
@@ -146,8 +198,7 @@
           retval = rng_base + (rng_nelem - 1) * rng_inc;
 
           // See the note in the matrix_value method above.
-
-          if (retval < rng_limit)
+          if (retval <= rng_limit)
             retval = rng_limit;
         }
 
@@ -166,8 +217,7 @@
           retval = rng_base + (rng_nelem - 1) * rng_inc;
 
           // See the note in the matrix_value method above.
-
-          if (retval > rng_limit)
+          if (retval >= rng_limit)
             retval = rng_limit;
         }
       else
@@ -268,7 +318,7 @@
   if (dim == 1)
     {
       if (mode == ASCENDING)
-          retval.sort_internal (sidx, true);
+        retval.sort_internal (sidx, true);
       else if (mode == DESCENDING)
         retval.sort_internal (sidx, false);
     }
@@ -281,10 +331,10 @@
 sortmode
 Range::is_sorted (sortmode mode) const
 {
-  if (rng_nelem > 1 && rng_inc < 0)
+  if (rng_nelem > 1 && rng_inc > 0)
+    mode = (mode == DESCENDING) ? UNSORTED : ASCENDING;
+  else if (rng_nelem > 1 && rng_inc < 0)
     mode = (mode == ASCENDING) ? UNSORTED : DESCENDING;
-  else if (rng_nelem > 1 && rng_inc > 0)
-    mode = (mode == DESCENDING) ? UNSORTED : ASCENDING;
   else
     mode = mode ? mode : ASCENDING;
 
@@ -298,12 +348,15 @@
   double increment = a.inc ();
   octave_idx_type num_elem = a.nelem ();
 
-  for (octave_idx_type i = 0; i < num_elem-1; i++)
-    os << b + i * increment << " ";
+  if (num_elem > 1)
+    {
+      // First element must be the base *exactly* (-0).
+      os << b << " ";
+      for (octave_idx_type i = 1; i < num_elem-1; i++)
+        os << b + i * increment << " ";
+    }
 
-  // Prevent overshoot.  See comment in the matrix_value method
-  // above.
-
+  // Prevent overshoot.  See comment in the matrix_value method above.
   os << (increment > 0 ? a.max () : a.min ()) << "\n";
 
   return os;
@@ -502,7 +555,7 @@
             n_elt++;
         }
 
-      retval = (n_elt >= std::numeric_limits<octave_idx_type>::max () - 1) ? -1 : n_elt;
+      retval = (n_elt < std::numeric_limits<octave_idx_type>::max () - 1) ? n_elt : -1;
     }
 
   return retval;

--- a/liboctave/array/Range.h
+++ b/liboctave/array/Range.h
@@ -50,7 +50,13 @@
       rng_nelem (nelem_internal ()), cache () { }
 
   // For operators' usage (to preserve element count).
-  Range (double b, double i, octave_idx_type n);
+  Range (double b, double i, octave_idx_type n)
+    : rng_base (b), rng_limit (b + (n-1) * i), rng_inc (i),
+      rng_nelem (n), cache ()
+    {
+      if (! xfinite (b) || ! xfinite (i))
+        rng_nelem = -2;
+    }
 
   double base (void) const { return rng_base; }
   double limit (void) const { return rng_limit; }
@@ -80,14 +86,7 @@
 
   double checkelem (octave_idx_type i) const;
 
-  double elem (octave_idx_type i) const
-    {
-#if defined (BOUNDS_CHECKING)
-      return checkelem (i);
-#else
-      return rng_base + rng_inc * i;
-#endif
-    }
+  double elem (octave_idx_type i) const;
 
   Array<double> index (const idx_vector& i) const;
 

--- a/test/range.tst
+++ b/test/range.tst
@@ -72,3 +72,38 @@
 %!assert ([ r ; uint32(z)          ], uint32(expect))
 %!assert ([ r ; uint64(z)          ], uint64(expect))
 
+## Test corner cases of ranges (base and limit)
+
+%!shared r, rrev, rneg
+%! r = -0:3;
+%! rrev = 3:-1:-0;
+%! rneg = -3:-0;
+
+%!assert (full (r), [-0 1 2 3])
+%!assert (signbit (full (r)), logical ([1 0 0 0]))
+%!assert (r(1), -0)
+%!assert (signbit (r(1)), true)
+%!assert (signbit (r(1:2)), logical ([1 0]))
+%!assert (signbit (r(2:-1:1)), logical ([0 1]))
+%!assert (signbit (r([2 1 1 3])), logical ([0 1 1 0]))
+
+%!assert (full (rrev), [3 2 1 -0])
+%!assert (signbit (full (rrev)), logical ([0 0 0 1]))
+%!assert (rrev(4), -0)
+%!assert (signbit (rrev(4)), true)
+%!assert (signbit (rrev(3:4)), logical ([0 1]))
+%!assert (signbit (rrev(4:-1:3)), logical ([1 0]))
+%!assert (signbit (rrev([1 4 4 2])), logical ([0 1 1 0]))
+
+%!assert (min (r), -0)
+%!assert (signbit (min (r)), true)
+%!assert (min (rrev), -0)
+%!assert (signbit (min (rrev)), true)
+
+%!assert (max (rneg), -0)
+%!assert (signbit (min (rneg)), true)
+
+%!assert (sort (r, "descend"), [3 2 1 -0])
+%!assert (signbit (sort (r, "descend")), logical ([0 0 0 1]))
+%!assert (signbit (sort (rrev, "ascend")), logical ([1 0 0 0]))
+