--- a/liboctave/oct-rand.cc
+++ b/liboctave/oct-rand.cc
@@ -419,19 +419,120 @@
   return retval;
 }
 
-Matrix
-octave_rand::do_matrix (octave_idx_type n, octave_idx_type m, double a)
+float
+octave_rand::do_float_scalar (float a)
 {
-  Matrix retval;
+  float retval = 0.0;
+
+  if (use_old_generators)
+    {
+      double da = a;
+      double dretval;
+      switch (current_distribution)
+        {
+        case uniform_dist:
+          F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, dretval);
+          break;
+
+        case normal_dist:
+          F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, dretval);
+          break;
+
+        case expon_dist:
+          F77_FUNC (dgenexp, DGENEXP) (1.0, dretval);
+          break;
 
-  if (n >= 0 && m >= 0)
-    {
-      retval.clear (n, m);
+        case poisson_dist:
+          if (da < 0.0 || xisnan(da) || xisinf(da))
+            dretval = octave_NaN;
+          else
+            {
+              // workaround bug in ignpoi, by calling with different Mu
+              F77_FUNC (dignpoi, DIGNPOI) (da + 1, dretval);
+              F77_FUNC (dignpoi, DIGNPOI) (da, dretval);
+            }
+          break;
 
-      if (n > 0 && m > 0)
-        fill (retval.capacity(), retval.fortran_vec(), a);
+        case gamma_dist:
+          if (da <= 0.0 || xisnan(da) || xisinf(da))
+            retval = octave_NaN;
+          else
+            F77_FUNC (dgengam, DGENGAM) (1.0, da, dretval);
+          break;
+
+        default:
+          (*current_liboctave_error_handler)
+            ("rand: invalid distribution ID = %d", current_distribution);
+          break;
+        }
+      retval = dretval;
     }
   else
+    {
+      switch (current_distribution)
+        {
+        case uniform_dist:
+          retval = oct_float_randu ();
+          break;
+
+        case normal_dist:
+          retval = oct_float_randn ();
+          break;
+
+        case expon_dist:
+          retval = oct_float_rande ();
+          break;
+
+        case poisson_dist:
+          // Keep poisson distribution in double precision for accuracy
+          retval = oct_randp (a);
+          break;
+
+        case gamma_dist:
+          retval = oct_float_randg (a);
+          break;
+
+        default:
+          (*current_liboctave_error_handler)
+            ("rand: invalid distribution ID = %d", current_distribution);
+          break;
+        }
+
+      save_state ();
+    }
+
+  return retval;
+}
+
+Array<double>
+octave_rand::do_vector (octave_idx_type n, double a)
+{
+  Array<double> retval;
+
+  if (n > 0)
+    {
+      retval.clear (n, 1);
+
+      fill (retval.capacity(), retval.fortran_vec(), a);
+    }
+  else if (n < 0)
+    (*current_liboctave_error_handler) ("rand: invalid negative argument");
+
+  return retval;
+}
+
+Array<float>
+octave_rand::do_float_vector (octave_idx_type n, float a)
+{
+  Array<float> retval;
+
+  if (n > 0)
+    {
+      retval.clear (n, 1);
+
+      fill (retval.capacity(), retval.fortran_vec(), a);
+    }
+  else if (n < 0)
     (*current_liboctave_error_handler) ("rand: invalid negative argument");
 
   return retval;
@@ -452,19 +553,17 @@
   return retval;
 }
 
-Array<double>
-octave_rand::do_vector (octave_idx_type n, double a)
+FloatNDArray
+octave_rand::do_float_nd_array (const dim_vector& dims, float a)
 {
-  Array<double> retval;
+  FloatNDArray retval;
 
-  if (n > 0)
+  if (! dims.all_zero ())
     {
-      retval.clear (n, 1);
+      retval.clear (dims);
 
-      fill (retval.capacity (), retval.fortran_vec (), a);
+      fill (retval.capacity(), retval.fortran_vec(), a);
     }
-  else if (n < 0)
-    (*current_liboctave_error_handler) ("rand: invalid negative argument");
 
   return retval;
 }
@@ -693,3 +792,94 @@
 
   return;
 }
+
+void
+octave_rand::fill (octave_idx_type len, float *v, float a)
+{
+  if (len < 1)
+    return;
+
+  switch (current_distribution)
+    {
+    case uniform_dist:
+      if (use_old_generators)
+        {
+#define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x)
+          MAKE_RAND (len);
+#undef RAND_FUNC
+        }
+      else
+        oct_fill_float_randu (len, v);
+      break;
+
+    case normal_dist:
+      if (use_old_generators)
+        {
+#define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x)
+          MAKE_RAND (len);
+#undef RAND_FUNC
+        }
+      else
+        oct_fill_float_randn (len, v);
+      break;
+
+    case expon_dist:
+      if (use_old_generators)
+        {
+#define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x)
+          MAKE_RAND (len);
+#undef RAND_FUNC
+        }
+      else
+        oct_fill_float_rande (len, v);
+      break;
+
+    case poisson_dist:
+      if (use_old_generators)
+        {
+          double da = a;
+          if (da < 0.0 || xisnan(da) || xisinf(da))
+#define RAND_FUNC(x) x = octave_NaN;
+            MAKE_RAND (len);
+#undef RAND_FUNC
+          else
+            {
+              // workaround bug in ignpoi, by calling with different Mu
+              double tmp;
+              F77_FUNC (dignpoi, DIGNPOI) (da + 1, tmp);
+#define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (da, x)
+                MAKE_RAND (len);
+#undef RAND_FUNC
+            }
+        }
+      else
+        oct_fill_float_randp (a, len, v);
+      break;
+
+    case gamma_dist:
+      if (use_old_generators)
+        {
+          double da = a;
+          if (da <= 0.0 || xisnan(da) || xisinf(da))
+#define RAND_FUNC(x) x = octave_NaN;
+            MAKE_RAND (len);
+#undef RAND_FUNC
+          else
+#define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, da, x)
+            MAKE_RAND (len);
+#undef RAND_FUNC
+        }
+      else
+        oct_fill_float_randg (a, len, v);
+      break;
+
+    default:
+      (*current_liboctave_error_handler)
+        ("rand: invalid distribution ID = %d", current_distribution);
+      break;
+    }
+
+  save_state ();
+
+  return;
+}