diff --git a/liboctave/array/Range.cc b/liboctave/array/Range.cc
--- a/liboctave/array/Range.cc
+++ b/liboctave/array/Range.cc
@@ -219,14 +219,38 @@ xinit (T base, T limit, T inc, bool reve
   // The following case also catches Inf values for increment when
   // there will be only one element.
 
-  if ((limit <= base && base + inc < limit)
-      || (limit >= base && base + inc > limit))
+  if (((limit <= base && base + inc < limit)
+       || (limit >= base && base + inc > limit))
+      && ! xteq (base + inc, limit))
     {
       final_val = base;
       nel = 1;
       return;
     }
 
+  // Catch integer ranges.
+  if (math::nint_big (base) == base && math::nint_big (inc) == inc
+      && ! math::isinf (limit))
+    {
+      nel = (limit - base + inc) / inc;
+      final_val = base + (nel - 1) * inc;
+
+      // The computed final value should never be beyond the specified
+      // limit of the range.
+      if ((inc > 0 && final_val > limit)
+          || (inc < 0 && final_val < limit))
+        {
+          nel--;
+          final_val = base + (nel - 1) * inc;
+        }
+
+      // Preserve -0.
+      if (final_val == limit)
+        final_val = limit;
+
+      return;
+    }
+
   // Any other calculations with Inf will give us either a NaN range
   // or an infinite nember of elements.
 
@@ -253,7 +277,12 @@ xinit (T base, T limit, T inc, bool reve
   // the number of elements and the final value in a way that attempts
   // to avoid rounding errors as much as possible.
 
+  (*current_liboctave_warning_with_id_handler)
+    ("Octave:floating-point-range",
+     "floating point ranges may produce unexpected results");
+
   nel = xnumel_internal (base, limit, inc);
+
   final_val = xfinal_value (base, limit, inc, nel);
 }
 
diff --git a/scripts/general/interpn.m b/scripts/general/interpn.m
--- a/scripts/general/interpn.m
+++ b/scripts/general/interpn.m
@@ -131,7 +131,8 @@ function vi = interpn (varargin)
     y = cell (1, nd);
     for i = 1 : nd
       x{i} = 1 : sz(i);
-      y{i} = 1 : (1 / (2 ^ m)) : sz(i);
+      pow2m = 2 ^ m;
+      y{i} = (pow2m : (sz(i) * pow2m)) / pow2m;
     endfor
     y{1} = y{1}.';
     [y{:}] = ndgrid (y{:});
diff --git a/test/range.tst b/test/range.tst
--- a/test/range.tst
+++ b/test/range.tst
@@ -153,18 +153,25 @@
 ## Test final value within eps of an integer (bug #46859)
 %!test <*46859>
 %! rng = 1 : (1001/250)/(1/250);
-%! assert (rng(end), 1001);
+%! assert (numel (rng), 1000);
+%! assert (rng(end), 1000);
 
 %!test <*46859>
 %! rng = 2000: -1 : (1001/250)/(1/250);
+%! assert (numel (rng), 1000);
 %! assert (rng(end), 1001);
 
-## This is not Matlab compatible (stops at 1000 with 999 elements)
-## Octave prefers the more intuitive "pure math" approach where
-## (1001/250) / (1/250) => (1001/250)*(250/1) => 1001.
 %!test <*46859>
 %! rng = 1 : (1001/250)/(1/250);
-%! assert (numel (1000));
+%! assert (numel (rng), (1000));
+%! assert (rng(end), 1000);
+
+%!assert <64692> ([-2:-eps], [-2, -1])
+%!assert <64692> ([2:-1:eps], [2, 1])
+%!assert <64692> ([-2:eps], [-2, -1, 0])
+%!assert <64692> ([2:-1:-eps], [2, 1, 0])
+%!assert <64692> ([1.8:0.05:1.9], [1.8, 1.85, 1.9])
+%!assert <64692> ([1.85:0.05:1.9], [1.85, 1.9])
 
 ## Combinations of exceptional values and a few normal ones.