--- a/libinterp/interpfcn/data.cc
+++ b/libinterp/interpfcn/data.cc
@@ -3342,6 +3342,11 @@
   return retval;
 }
 
+/*
+%% Debian bug #706376
+%!assert (isempty (speye(2^16)), false)
+*/
+
 DEFUN (isnumeric, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} isnumeric (@var{x})\n\

--- a/libinterp/octave-value/ov-base-sparse.cc
+++ b/libinterp/octave-value/ov-base-sparse.cc
@@ -278,6 +278,15 @@
 
 template <class T>
 bool
+octave_base_sparse<T>::is_empty (void) const
+{
+  dim_vector dv = dims ();
+
+  return (dv.any_zero ());
+}
+
+template <class T>
+bool
 octave_base_sparse<T>::print_as_scalar (void) const
 {
   dim_vector dv = dims ();

--- a/libinterp/octave-value/ov-base-sparse.h
+++ b/libinterp/octave-value/ov-base-sparse.h
@@ -72,6 +72,8 @@
 
   ~octave_base_sparse (void) { }
 
+  octave_idx_type numel (void) const { return dims ().safe_numel (); }
+
   octave_idx_type nnz (void) const { return matrix.nnz (); }
 
   octave_idx_type nzmax (void) const { return matrix.nzmax (); }
@@ -137,6 +139,8 @@
 
   bool is_defined (void) const { return true; }
 
+  bool is_empty (void) const;
+
   bool is_constant (void) const { return true; }
 
   bool is_true (void) const;

--- a/libinterp/octave-value/ov-base.h
+++ b/libinterp/octave-value/ov-base.h
@@ -331,7 +331,7 @@
 
   virtual bool is_defined (void) const { return false; }
 
-  bool is_empty (void) const { return numel () == 0; }
+  virtual bool is_empty (void) const { return numel () == 0; }
 
   virtual bool is_cell (void) const { return false; }
 

--- a/libinterp/octave-value/ov-bool-sparse.h
+++ b/libinterp/octave-value/ov-bool-sparse.h
@@ -83,9 +83,9 @@
 
   octave_base_value *try_narrowing_conversion (void);
 
-  // FIXME Adapt idx_vector to allow sparse logical indexing!!
+  // FIXME Adapt idx_vector to allow sparse logical indexing without overflow!!
   idx_vector index_vector (void) const
-    { return idx_vector (bool_array_value ()); }
+    { return idx_vector (matrix); }
 
   builtin_type_t builtin_type (void) const { return btyp_bool; }
 

--- a/liboctave/array/idx-vector.cc
+++ b/liboctave/array/idx-vector.cc
@@ -423,30 +423,27 @@
 }
 
 idx_vector::idx_vector_rep::idx_vector_rep (const Sparse<bool>& bnda)
-  : data (0), len (0), ext (0), aowner (0), orig_dims ()
+  : data (0), len (bnda.nnz ()), ext (0), aowner (0), orig_dims ()
 {
-  for (octave_idx_type i = 0, l = bnda.nnz (); i < l; i++)
-    if (bnda.data (i)) len++;
+  const dim_vector dv = bnda.dims ();
 
-  dim_vector dv = bnda.dims ();
-
-  orig_dims = ((dv.length () == 2 && dv(0) == 1)
-               ? dim_vector (1, len) : orig_dims = dim_vector (len, 1));
+  if (! dv.all_zero ())
+    orig_dims = ((dv.length () == 2 && dv(0) == 1)
+                 ? dim_vector (1, len) : dim_vector (len, 1));
 
   if (len != 0)
     {
       octave_idx_type *d = new octave_idx_type [len];
 
-      octave_idx_type nnz = bnda.nnz ();
+      octave_idx_type k = 0;
+      octave_idx_type nc = bnda.cols ();
+      octave_idx_type nr = bnda.rows ();
 
-      octave_idx_type k = 0;
-      // FIXME: I hope this is OK, i.e. the element iterated this way are correctly ordered.
-      for (octave_idx_type i = 0; i < nnz; i++)
-        {
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = bnda.cidx(j); i < bnda.cidx(j+1); i++)
           if (bnda.data (i))
-            d[k++] = bnda.cidx (i) + bnda.rows () * bnda.ridx (i);
-        }
-
+            d[k++] = j * nr + bnda.ridx (i);
+ 
       data = d;
 
       ext = d[k-1] + 1;

--- a/liboctave/util/lo-array-gripes.cc
+++ b/liboctave/util/lo-array-gripes.cc
@@ -130,7 +130,11 @@
   const char *err_id = error_id_invalid_index;
 
   (*current_liboctave_error_with_id_handler)
-    (err_id, "subscript indices must be either positive integers or logicals");
+#ifdef USE_64_BIT_IDX_T
+    (err_id, "subscript indices must be either positive integers less than 2^63 or logicals");
+#else
+    (err_id, "subscript indices must be either positive integers less than 2^31 or logicals");
+#endif
 }
 
 // FIXME -- the following is a common error message to resize,

--- a/scripts/sparse/private/__sprand_impl__.m
+++ b/scripts/sparse/private/__sprand_impl__.m
@@ -55,9 +55,27 @@
 
   mn = m*n;
   k = round (d*mn);
-  idx = randperm (mn, k);
+  if (mn > sizemax ())
+    ## randperm will overflow, so use alternative methods
+
+    idx = unique (fix (rand (min (k*1.01, k+10), 1) * mn)) + 1;
 
-  [i, j] = ind2sub ([m, n], idx);
+    ## idx contains random numbers in [1,mn]
+    ## generate 1% or 10 more random values than necessary in order to
+    ## reduce the probability that there are less than k distinct
+    ## values; maybe a better strategy could be used but I don't think
+    ## it's worth the price
+    
+    ## actual number of entries in S
+    k = min (length (idx), k);
+    j = floor ((idx(1:k) - 1) / m);
+    i = idx(1:k) - j * m;
+  else
+    idx = randperm (mn, k);
+    [i, j] = ind2sub ([m, n], idx);
+  endif
+
+
   S = sparse (i, j, randfun (k, 1), m, n);
 
 endfunction

--- a/scripts/sparse/sprand.m
+++ b/scripts/sparse/sprand.m
@@ -81,3 +81,7 @@
 %!error sprand (3, 3, -1)
 %!error sprand (3, 3, 2)
 
+%% Test very large, very low density matrix doesn't fail 
+%!test
+%! s = sprand(1e6,1e6,1e-7);
+

--- a/scripts/sparse/sprandn.m
+++ b/scripts/sparse/sprandn.m
@@ -72,3 +72,6 @@
 %!error sprandn (3, 3, -1)
 %!error sprandn (3, 3, 2)
 
+%% Test very large, very low density matrix doesn't fail 
+%!test
+%! s = sprandn(1e6,1e6,1e-7);