diff -r 230c1e2a678d libinterp/corefcn/data.cc
--- a/libinterp/corefcn/data.cc	Sat Apr 04 14:12:10 2015 -0400
+++ b/libinterp/corefcn/data.cc	Sun May 03 16:47:23 2015 +0200
@@ -68,6 +68,7 @@
 #include "variables.h"
 #include "pager.h"
 #include "xnorm.h"
+#include "oct-rand.h"
 
 #if ! defined (CLOCKS_PER_SEC)
 #if defined (CLK_TCK)
@@ -574,6 +575,124 @@
 %!assert (all (log2 (pow2 (-1074:1023)) == -1074:1023))
 */
 
+
+DEFUN (__sprand_impl_basearg__, args, ,
+    "-*- texinfo -*-\n\
+@deftypefn  {Function File} {} __sprand_impl_basearg__ (@var{s}, @var{randfunname})\n\
+@deftypefnx {Function File} {} __sprand_impl_basearg__ (@var{m}, @var{n}, @var{d}, @var{randfunname})\n\
+Undocumented internal function.\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  if(args.length () != 2 && args.length () != 4)
+    {
+      error("__sprand_impl_basearg__: wrong number of arguments");
+      return retval;
+    }
+
+  std::string funname;
+  if(args.length () == 2)
+    {
+      funname = args(1).string_value();
+    }
+  else
+    {
+      funname = args(3).string_value();
+    }
+
+  std::string rand_dist = octave_rand::distribution ();//FIXME should add a stack unwind protection?
+
+  //FIXME: manage more than two distributions
+  if( funname == "rand" )
+    {
+      octave_rand::distribution ("uniform");
+    }
+  else if(funname == "randn")
+    {
+      octave_rand::distribution ("normal");
+    }
+  else
+    {
+      error ("__sprand_impl_basearg__: %s random function not recognized",
+          funname.c_str ());
+      octave_rand::distribution (rand_dist);
+      return retval;
+    }
+
+  SparseMatrix r;
+
+
+  if(args.length () == 2)
+    {
+      double * data;
+      octave_idx_type nelem;
+
+
+    //FIXME only working for dense matrices and sparse matrices, no bool or char
+      if(args(0).is_sparse_type ())
+        {
+          const SparseMatrix m = args(0).sparse_matrix_value ();
+          r = SparseMatrix (m);
+        }
+      else
+        {
+          const Matrix m = args(0).matrix_value ();
+          r = SparseMatrix (m);
+        }
+
+      r.maybe_compress();
+      data = r.data ();
+      nelem = r.nelem ();
+
+      octave_idx_type i;
+
+      Array<double> rd = octave_rand::vector(nelem);
+      for(i=0; i < nelem; i++)
+        {
+          data[i] = rd.xelem(i);
+        }
+    }
+  else if(args.length () == 4)
+    {
+      //  FIXME: declaring m and casting before checking for scalar and
+      //  integerness. Also is_idx_type_type does not exist
+      const octave_idx_type m = args(0).idx_type_value ();
+      const octave_idx_type n = args(1).idx_type_value ();
+      const double d = args(2).double_value ();
+
+      //original test: if (!(isscalar (m) && m == fix (m) && m > 0))
+      if (!(args(0).is_scalar_type() && double(m) == args(0).double_value() && m > 0))
+        {
+          error ("%s: M must be an integer greater than 0", funname.c_str ());
+          octave_rand::distribution (rand_dist);
+          return retval;
+        }
+      if (!(args(1).is_scalar_type() && double(n) == args(1).double_value() && n > 0))
+        {
+          error ("%s: N must be an integer greater than 0", funname.c_str ());
+          octave_rand::distribution (rand_dist);
+          return retval;
+        }
+
+      if (d < 0 || d > 1)
+        {
+          error ("%s: density D must be between 0 and 1", funname.c_str ());
+          octave_rand::distribution (rand_dist);
+          return retval;
+        }
+
+      r = SparseMatrix(
+            octave_rand::sparse_array ( dim_vector (m, n), d) );
+    }
+
+  octave_rand::distribution (rand_dist);
+
+  return retval = octave_value (r);
+}
+
+
+
 DEFUN (rem, args, ,
        "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} rem (@var{x}, @var{y})\n\
diff -r 230c1e2a678d liboctave/array/Sparse-f.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/array/Sparse-f.cc	Sun May 03 16:47:23 2015 +0200
@@ -0,0 +1,52 @@
+/*
+
+Copyright (C) 2004-2013 David Bateman
+Copyright (C) 1998-2004 Andy Adler
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+// Instantiate Sparse matrix of float values.
+
+#include "lo-mappers.h"
+#include "Sparse.h"
+#include "Sparse.cc"
+
+template <>
+bool
+sparse_ascending_compare<float> (float a, float b)
+{
+  return (xisnan (b) || (a < b));
+}
+
+template <>
+bool
+sparse_descending_compare<float> (float a, float b)
+{
+  return (xisnan (a) || (a > b));
+}
+
+INSTANTIATE_SPARSE (float, OCTAVE_API);
+
+#if 0
+template std::ostream& operator << (std::ostream&, const Sparse<float>&);
+#endif
diff -r 230c1e2a678d liboctave/array/module.mk
--- a/liboctave/array/module.mk	Sat Apr 04 14:12:10 2015 -0400
+++ b/liboctave/array/module.mk	Sun May 03 16:47:23 2015 +0200
@@ -115,6 +115,7 @@
   array/Sparse-b.cc \
   array/Sparse-C.cc \
   array/Sparse-d.cc \
+  array/Sparse-f.cc \
   array/uint16NDArray.cc \
   array/uint32NDArray.cc \
   array/uint64NDArray.cc \
diff -r 230c1e2a678d liboctave/numeric/oct-rand.cc
--- a/liboctave/numeric/oct-rand.cc	Sat Apr 04 14:12:10 2015 -0400
+++ b/liboctave/numeric/oct-rand.cc	Sun May 03 16:47:23 2015 +0200
@@ -568,6 +568,139 @@
   return retval;
 }
 
+template <typename T> Sparse<T>
+  octave_rand::do_sparse_array_templ (const dim_vector& dims, double d, double a /* = 1.0*/)
+{
+  Sparse<T> retval;
+
+  // Should use unwind_protect to restore distribuition?
+  //unwind_protect frame;
+  //frame.add_fcn (octave_rand::distribution, octave_rand::distribution ());
+
+  const octave_idx_type m = dims(0);
+  const octave_idx_type n = dims(1);
+
+  //FIXME: checks for overflows? but it should not happen because 2^128 is
+  //still inside double
+  const double mnd = static_cast<double>(m)*static_cast<double>(n)*d;
+
+  //FIXED: the previous implementation uses an approximate solution to control
+  //the density, it can fail tests
+  //when it does not generate k unique sample
+
+  if(gnulib::round ( mnd ) > dim_vector::dim_max())
+    {
+        (*current_liboctave_error_handler)
+            ("dimension too large for Octave's index type");
+        return retval;
+    }
+  const octave_idx_type k = gnulib::round ( mnd );
+
+  retval = Sparse<T> (m, n, k);
+  octave_idx_type i,j;
+
+  //FIXME: inefficient, iterate twice
+  //over the whole array, it is necessary to avoid switching
+  //the distribution at every step. Is it reasonable?
+
+  std::string rand_dist = octave_rand::distribution ();
+
+
+  octave_idx_type * ridx = retval.ridx ();
+  octave_idx_type * cidx = retval.cidx ();
+  T * data = retval.data ();
+
+  //build the uniform random column and row vectors
+
+  //FIXME: should call save state before switching distribution?
+  //save_state();
+
+  octave_rand::distribution ("uniform");
+
+  double dretval;
+
+  //filling the column buckets with samples
+  for(i = 0; i < k; i++)
+    {
+      do
+        {
+          if (use_old_generators)
+            {
+              F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, dretval);
+            }
+          else
+            {
+              dretval = oct_randu ();
+            }
+          j = static_cast<octave_idx_type>
+              ( gnulib::floor ( dretval * n) );//select a random bucket j
+        }
+      while(cidx[j+1] == m);//if the bucket is already full, reselect
+
+      cidx[j+1]++;//increase the bucket counter
+    }
+
+  //adapting to ccs notation
+  for(i = 0; i < n; i++)
+    {
+      cidx[i+1] += cidx[i];
+    }
+
+  for(i = 0; i < n; i++)
+    {
+
+      const octave_idx_type cn = cidx[i+1] - cidx[i];
+
+      std::set<octave_idx_type> rset;
+
+      //FIXME: build ordered set of sampled rows. in case of a dense, large row can be bingo sort slow
+      do
+        {
+          if (use_old_generators)
+            {
+              F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, dretval);
+            }
+          else
+            {
+              dretval = oct_randu ();
+            }
+          rset.insert(gnulib::floor ( dretval * m));
+        }
+      while( static_cast<octave_idx_type>(rset.size ()) < cn);
+
+      std::set<octave_idx_type>::const_iterator it = rset.begin();
+
+      for(j = cidx[i]; j < cidx[i+1]; j++)
+      {
+        ridx[j] = *it;
+        ++it;
+      }
+    }
+
+  //FIXME: should call save state before switching distribution?
+  //save_state();
+
+
+  //restore the distribution, and fill the data
+  octave_rand::distribution (rand_dist);
+
+  fill(k, data, a);
+
+  return retval;
+}
+
+Sparse<double>
+  octave_rand::do_sparse_array (const dim_vector& dims, double d, double a /* = 1.0*/)
+{
+  return do_sparse_array_templ<double>(dims,d,a);
+}
+
+Sparse<float>
+  octave_rand::do_float_sparse_array (const dim_vector& dims, double d, double a /*= 1.0*/)
+{
+  return do_sparse_array_templ<float>(dims,d,a);
+}
+
 // Make the random number generator give us a different sequence every
 // time we start octave unless we specifically set the seed.  The
 // technique used below will cycle monthly, but it it does seem to
diff -r 230c1e2a678d liboctave/numeric/oct-rand.h
--- a/liboctave/numeric/oct-rand.h	Sat Apr 04 14:12:10 2015 -0400
+++ b/liboctave/numeric/oct-rand.h	Sun May 03 16:47:23 2015 +0200
@@ -30,6 +30,7 @@
 #include "dNDArray.h"
 #include "fNDArray.h"
 #include "lo-ieee.h"
+#include "Sparse.h"
 
 class
 OCTAVE_API
@@ -170,6 +171,27 @@
                           : FloatNDArray ();
   }
 
+  // Return an sparse array of numbers from the sequence,
+  // uniformly selected from each row and column,
+  // with values according to the current distribution.
+  static Sparse<double>
+    sparse_array (const dim_vector& dims, double density, double a = 1.0)
+  {
+    return instance_ok () ? instance->do_sparse_array (dims, density, a)
+                            : Sparse<double> ();
+  }
+
+
+  // Return an sparse array of numbers from the sequence,
+  // uniformly selected from each row and column,
+  // with values according to the current distribution.
+  static Sparse<float>
+    float_sparse_array (const dim_vector& dims,  double density,float a = 1.0)
+  {
+    return instance_ok () ? instance->do_float_sparse_array (dims, density, a)
+                            : Sparse<float> ();
+  }
+
 private:
 
   static octave_rand *instance;
@@ -251,6 +273,22 @@
   // filled in column major order.
   FloatNDArray do_float_nd_array (const dim_vector& dims, float a = 1.);
 
+  template <typename T> Sparse<T>
+    do_sparse_array_templ (const dim_vector& dims, double d, double = 1.0);
+
+  // Return an sparse array of numbers from the sequence,
+  // uniformly selected from each row and column,
+  // with values according to the current distribution.
+  Sparse<double>
+    do_sparse_array (const dim_vector& dims, double d, double a = 1.0);
+
+  // Return an sparse array of numbers from the sequence,
+  // uniformly selected from each row and column,
+  // with values according to the current distribution.
+  Sparse<float>
+    do_float_sparse_array (const dim_vector& dims, double d, double a = 1.0);
+
+
   // Some helper functions.
 
   void initialize_ranlib_generators (void);
diff -r 230c1e2a678d scripts/sparse/module.mk
--- a/scripts/sparse/module.mk	Sat Apr 04 14:12:10 2015 -0400
+++ b/scripts/sparse/module.mk	Sun May 03 16:47:23 2015 +0200
@@ -1,7 +1,7 @@
 FCN_FILE_DIRS += sparse
 
 sparse_PRIVATE_FCN_FILES = \
-  sparse/private/__sprand_impl__.m
+  sparse/private/__sprand_impl_6arg__.m
 
 sparse_FCN_FILES = \
   sparse/bicg.m \
diff -r 230c1e2a678d scripts/sparse/private/__sprand_impl_6arg__.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/sparse/private/__sprand_impl_6arg__.m	Sun May 03 16:47:23 2015 +0200
@@ -0,0 +1,106 @@
+## Copyright (C) 2004-2013 Paul Kienzle
+## Copyright (C) 2012 Jordi Gutiérrez Hermoso
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+##
+## Original version by Paul Kienzle distributed as free software in the
+## public domain.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} __sprand_impl_6arg__ (@var{m}, @var{n}, @var{d}, @var{rc}, @var{fcnname}, @var{randfun})
+## Undocumented internal function.
+## @end deftypefn
+
+## Actual implementation of sprand and sprandn happens here.
+
+function S = __sprand_impl_6arg__ (varargin)
+
+
+	[m, n, d, rc, fcnname, randfun] = deal (varargin{:});
+
+	if (! (isscalar (m) && m == fix (m) && m > 0))
+		error ("%s: M must be an integer greater than 0", fcnname);
+	endif
+	if (! (isscalar (n) && n == fix (n) && n > 0))
+		error ("%s: N must be an integer greater than 0", fcnname);
+	endif
+	if (d < 0 || d > 1)
+		error ("%s: density D must be between 0 and 1", fcnname);
+	endif
+
+	## Create a matrix with specified reciprocal condition number.
+
+	if (! isscalar (rc) && ! isvector (rc))
+		error ("%s: RC must be a scalar or vector", fcnname);
+	endif
+
+	## We want to reverse singular valued decomposition A=U*S*V'.
+	## First, first S is constructed and then U = U1*U2*..Un and
+	## V' = V1*V2*..Vn are seen as Jacobi rotation matrices with angles and
+	## planes of rotation randomized.  Repeatedly apply rotations until the
+	## required density for A is achieved.
+
+	if (isscalar (rc))
+		if (rc < 0 || rc > 1)
+			error ("%s: reciprocal condition number RC must be between 0 and 1", fcnname);
+		endif
+		## Reciprocal condition number is ratio of smallest SV to largest SV
+		## Generate singular values randomly and sort them to build S
+		## Random singular values in range [rc, 1].
+		v = rand (1, min (m,n)) * (1 - rc) + rc;
+		v(1) = 1;
+		v(end) = rc;
+		v = sort (v, "descend");
+		S = sparse (diag (v, m, n));
+	else
+		## Only the min (m, n) greater singular values from rc vector are used.
+		if (length (rc) > min (m,n))
+			rc = rc(1:min(m, n));
+		endif
+		S = sparse (diag (sort (rc, "descend"), m, n));
+	endif
+
+	Uinit = speye (m);
+	Vinit = speye (n);
+	k = round (d*m*n);
+	while (nnz (S) < k)
+		if (m > 1)
+			## Construct U randomized rotation matrix
+			rot_angleu = 2 * pi * rand ();
+			cu = cos (rot_angleu); su = sin (rot_angleu);
+			rndtmp = randperm (m, 2);
+			i = rndtmp(1); j = rndtmp(2);
+			U = Uinit;
+			U(i, i) = cu; U(i, j) = -su;
+			U(j, i) = su; U(j, j) = cu;
+			S = U * S;
+		endif
+		if (n > 1)
+			## Construct V' randomized rotation matrix
+			rot_anglev = 2 * pi * rand ();
+			cv = cos (rot_anglev); sv = sin (rot_anglev);
+			rndtmp = randperm (n, 2);
+			i = rndtmp(1); j = rndtmp(2);
+			V = Vinit;
+			V(i, i) = cv;  V(i, j) = sv;
+			V(j, i) = -sv; V(j, j) = cv;
+			S = S * V;
+		endif
+	endwhile
+
+endfunction
+
diff -r 230c1e2a678d scripts/sparse/sprand.m
--- a/scripts/sparse/sprand.m	Sat Apr 04 14:12:10 2015 -0400
+++ b/scripts/sparse/sprand.m	Sun May 03 16:47:23 2015 +0200
@@ -53,11 +53,11 @@
 function s = sprand (m, n, d, rc)
 
   if (nargin == 1 )
-    s = __sprand_impl__ (m, @rand);
+    s = __sprand_impl_basearg__ (m, "rand");
   elseif ( nargin == 3)
-    s = __sprand_impl__ (m, n, d, "sprand", @rand);
+    s = __sprand_impl_basearg__ (m, n, d, "rand");
   elseif (nargin == 4)
-    s = __sprand_impl__ (m, n, d, rc, "sprand", @rand);
+    s = __sprand_impl_6arg__ (m, n, d, rc, "sprand", @rand);
   else
     print_usage ();
   endif
diff -r 230c1e2a678d scripts/sparse/sprandn.m
--- a/scripts/sparse/sprandn.m	Sat Apr 04 14:12:10 2015 -0400
+++ b/scripts/sparse/sprandn.m	Sun May 03 16:47:23 2015 +0200
@@ -45,11 +45,11 @@
 function s = sprandn (m, n, d, rc)
 
   if (nargin == 1 )
-    s = __sprand_impl__ (m, @randn);
+    s = __sprand_impl_basearg__ (m, "randn");
   elseif ( nargin == 3)
-    s = __sprand_impl__ (m, n, d, "sprandn", @randn);
+    s = __sprand_impl_basearg__ (m, n, d, "randn");
   elseif (nargin == 4)
-    s = __sprand_impl__ (m, n, d, rc, "sprandn", @randn);
+    s = __sprand_impl_6arg__ (m, n, d, rc, "sprand", @rand);
   else
     print_usage ();
   endif
