Mercurial > octave-nkf

changeset 14655:43db83eff9db

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293) * oct-rand.cc (octave_rand:do_matrix): Remove method. (octave_rand::do_float_scalar, octave_rand::do_float_nd_array, octave_rand::do_float_vector): New methods. * oct-rand.h (octave_rand:do_matrix, octave_rand::matrix): Remove methods. (octave_rand::do_float_scalar, octave_rand::do_float_nd_array, octave_rand::do_float_vector, octave_rand::float_scalar, octave_rand::do_nd_array, octave_rand::float_vector): Declare new methods. * randgamma.c (oct_fill_float_randg, oct_float_randg): New functions. * randgamma.h (oct_fill_float_randg, oct_float_randg): Declare new functions. * randpoisson.c (oct_fill_float_randp, oct_float_randp): New functions. (poisson_cdf_lookup_float, poisson_rejection_float): New static functions. * randpoisson.h (oct_fill_float_randp, oct_float_randp): Declare new functions. * randmtzig.c (randu64) : Remove static function. (ALLBITS): Remove compile flag logic. (randu32): Change return type to float. (oct_float_randu, oct_float_randn, oct_float_rande, oct_fill_float_randu, oct_fill_float_randn, oct_fill_float_rande): New functions. (create_ziggurat_float_tables): New static function. * randmtzig.h (oct_float_randu, oct_float_randn, oct_float_rande, oct_fill_float_randu, oct_fill_float_randn, oct_fill_float_rande): Declare new functions. * rand.cc (do_rand): Add logic to parse "double" or "single" as final argument.
author David Bateman <dbateman@free.fr>
date Sat, 19 May 2012 00:14:59 +0200
parents 389f49a88656
children 8d2ce821e38a
files liboctave/oct-rand.cc liboctave/oct-rand.h liboctave/randgamma.c liboctave/randgamma.h liboctave/randmtzig.c liboctave/randmtzig.h liboctave/randpoisson.c liboctave/randpoisson.h src/DLD-FUNCTIONS/rand.cc
diffstat 9 files changed, 747 insertions(+), 104 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/oct-rand.cc	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/oct-rand.cc	Sat May 19 00:14:59 2012 +0200
@@ -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;
+}
--- a/liboctave/oct-rand.h	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/oct-rand.h	Sat May 19 00:14:59 2012 +0200
@@ -27,8 +27,8 @@
 #include <string>
 
 #include "dColVector.h"
-#include "dMatrix.h"
 #include "dNDArray.h"
+#include "fNDArray.h"
 #include "lo-ieee.h"
 
 class
@@ -136,13 +136,24 @@
     return instance_ok () ? instance->do_scalar (a) : octave_NaN;
   }
 
-  // Return a matrix of numbers from the sequence, filled in column
-  // major order.
-  static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.0)
+  // Return the next number from the sequence.
+  static float float_scalar (float a = 1.0)
   {
-    return instance_ok () ? instance->do_matrix (r, c, a) : Matrix ();
+    return instance_ok () ? instance->do_float_scalar (a) : octave_Float_NaN;
   }
 
+  // Return an array of numbers from the sequence.
+  static Array<double> vector (octave_idx_type n, double a = 1.0)
+   {
+     return instance_ok () ? instance->do_vector (n, a) : Array<double> ();
+   }
+
+  // Return an array of numbers from the sequence.
+  static Array<float> float_vector (octave_idx_type n, float a = 1.0)
+   {
+     return instance_ok () ? instance->do_float_vector (n, a) : Array<float> ();
+   }
+
   // Return an N-dimensional array of numbers from the sequence,
   // filled in column major order.
   static NDArray nd_array (const dim_vector& dims, double a = 1.0)
@@ -150,10 +161,12 @@
     return instance_ok () ? instance->do_nd_array (dims, a) : NDArray ();
   }
 
-  // Return an array of numbers from the sequence.
-  static Array<double> vector (octave_idx_type n, double a = 1.0)
+
+  // Return an N-dimensional array of numbers from the sequence,
+  // filled in column major order.
+  static FloatNDArray float_nd_array (const dim_vector& dims, float a = 1.0)
   {
-    return instance_ok () ? instance->do_vector (n, a) : Array<double> ();
+    return instance_ok () ? instance->do_float_nd_array (dims, a) : FloatNDArray ();
   }
 
 private:
@@ -220,16 +233,22 @@
   // Return the next number from the sequence.
   double do_scalar (double a = 1.);
 
-  // Return a matrix of numbers from the sequence, filled in column
-  // major order.
-  Matrix do_matrix (octave_idx_type r, octave_idx_type c, double a = 1.);
+  // Return the next number from the sequence.
+  float do_float_scalar (float a = 1.);
+
+  // Return an array of numbers from the sequence.
+  Array<double> do_vector (octave_idx_type n, double a = 1.);
+
+  // Return an array of numbers from the sequence.
+  Array<float> do_float_vector (octave_idx_type n, float a = 1.);
 
   // Return an N-dimensional array of numbers from the sequence,
   // filled in column major order.
   NDArray do_nd_array (const dim_vector& dims, double a = 1.);
 
-  // Return an array of numbers from the sequence.
-  Array<double> do_vector (octave_idx_type n, double a = 1.);
+  // Return an N-dimensional array of numbers from the sequence,
+  // filled in column major order.
+  FloatNDArray do_float_nd_array (const dim_vector& dims, float a = 1.);
 
   // Some helper functions.
 
@@ -248,6 +267,8 @@
   void switch_to_generator (int dist);
 
   void fill (octave_idx_type len, double *v, double a);
+
+  void fill (octave_idx_type len, float *v, float a);
 };
 
 #endif
--- a/liboctave/randgamma.c	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/randgamma.c	Sat May 19 00:14:59 2012 +0200
@@ -141,3 +141,59 @@
   oct_fill_randg(a,1,&ret);
   return ret;
 }
+
+#undef NAN
+#undef RUNI
+#undef RNOR
+#undef REXP
+#define NAN octave_Float_NaN
+#define RUNI oct_float_randu()
+#define RNOR oct_float_randn()
+#define REXP oct_float_rande()
+
+void
+oct_fill_float_randg (float a, octave_idx_type n, float *r)
+{
+  octave_idx_type i;
+  /* If a < 1, start by generating gamma(1+a) */
+  const float d =  (a < 1. ? 1.+a : a) - 1./3.;
+  const float c = 1./sqrt(9.*d);
+
+  /* Handle invalid cases */
+  if (a <= 0 || INFINITE(a))
+    {
+      for (i=0; i < n; i++)
+        r[i] = NAN;
+      return;
+    }
+
+  for (i=0; i < n; i++)
+    {
+      float x, xsq, v, u;
+    frestart:
+      x = RNOR;
+      v = (1+c*x);
+      v *= v*v;
+      if (v <= 0)
+        goto frestart; /* rare, so don't bother moving up */
+      u = RUNI;
+      xsq = x*x;
+      if (u >= 1.-0.0331*xsq*xsq && log(u) >= 0.5*xsq + d*(1-v+log(v)))
+        goto frestart;
+      r[i] = d*v;
+    }
+  if (a < 1)
+    { /* Use gamma(a) = gamma(1+a)*U^(1/a) */
+      /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
+      for (i = 0; i < n; i++)
+        r[i] *= exp(-REXP/a);
+    }
+}
+
+float
+oct_float_randg (float a)
+{
+  float ret;
+  oct_fill_float_randg(a,1,&ret);
+  return ret;
+}
--- a/liboctave/randgamma.h	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/randgamma.h	Sat May 19 00:14:59 2012 +0200
@@ -32,6 +32,9 @@
 extern OCTAVE_API double oct_randg (double a);
 extern OCTAVE_API void oct_fill_randg (double a, octave_idx_type n, double *p);
 
+extern OCTAVE_API float oct_float_randg (float a);
+extern OCTAVE_API void oct_fill_float_randg (float a, octave_idx_type n, float *p);
+
 #ifdef  __cplusplus
 }
 #endif
--- a/liboctave/randmtzig.c	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/randmtzig.c	Sat May 19 00:14:59 2012 +0200
@@ -64,10 +64,6 @@
    http://www.math.keio.ac.jp/matumoto/emt.html
    email: matumoto@math.keio.ac.jp
 
-   * 2006-04-01 David Bateman
-   * * convert for use in octave, declaring static functions only used
-   *   here and adding oct_ to functions visible externally
-   * * inverse sense of ALLBITS
    * 2004-01-19 Paul Kienzle
    * * comment out main
    * add init_by_entropy, get_state, set_state
@@ -87,6 +83,15 @@
    *
    * 2005-02-23 Paul Kienzle
    * * fix -DHAVE_X86_32 flag and add -DUSE_X86_32=0|1 for explicit control
+   *
+   * 2006-04-01 David Bateman
+   * * convert for use in octave, declaring static functions only used
+   *   here and adding oct_ to functions visible externally
+   * * inverse sense of ALLBITS
+   *
+   * 2012-05-18 David Bateman
+   * * Remove randu64 and ALLBIT option
+   * * Add the single precision generators
    */
 
 /*
@@ -96,9 +101,6 @@
    available.  This is not necessary if your architecture has
    /dev/urandom defined.
 
-   Compile with -DALLBITS to disable 53-bit random numbers. This is about
-   50% slower than using 32-bit random numbers.
-
    Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86.
    You can force X86 behaviour with -DUSE_X86_32=1, or suppress it with
    -DUSE_X86_32=0. You should also consider -march=i686 or similar for
@@ -126,21 +128,28 @@
    static uint32_t randi32(void)   returns 32-bit unsigned int
    static uint64_t randi53(void)   returns 53-bit unsigned int
    static uint64_t randi54(void)   returns 54-bit unsigned int
-   static uint64_t randi64(void)   returns 64-bit unsigned int
-   static double randu32(void)     returns 32-bit uniform in (0,1)
+   static float randu32(void)     returns 32-bit uniform in (0,1)
    static double randu53(void)     returns 53-bit uniform in (0,1)
 
    double oct_randu(void)       returns M-bit uniform in (0,1)
    double oct_randn(void)       returns M-bit standard normal
    double oct_rande(void)       returns N-bit standard exponential
 
+   float oct_float_randu(void)       returns M-bit uniform in (0,1)
+   float oct_float_randn(void)       returns M-bit standard normal
+   float oct_float_rande(void)       returns N-bit standard exponential
+
    === Array generators ===
    void oct_fill_randi32(octave_idx_type, uint32_t [])
    void oct_fill_randi64(octave_idx_type, uint64_t [])
+
    void oct_fill_randu(octave_idx_type, double [])
    void oct_fill_randn(octave_idx_type, double [])
    void oct_fill_rande(octave_idx_type, double [])
 
+   void oct_fill_float_randu(octave_idx_type, float [])
+   void oct_fill_float_randn(octave_idx_type, float [])
+   void oct_fill_float_rande(octave_idx_type, float [])
 */
 
 #if defined (HAVE_CONFIG_H)
@@ -180,6 +189,7 @@
 static int left = 1;
 static int initf = 0;
 static int initt = 1;
+static int inittf = 1;
 
 /* initializes state[MT_N] with a seed */
 void
@@ -378,34 +388,14 @@
 #endif
 }
 
-#if 0
-// FIXME -- this doesn't seem to be used anywhere; should it be removed?
-static uint64_t
-randi64 (void)
-{
-  const uint32_t lo = randi32();
-  const uint32_t hi = randi32();
-#if HAVE_X86_32
-  uint64_t u;
-  uint32_t *p = (uint32_t *)&u;
-  p[0] = lo;
-  p[1] = hi;
-  return u;
-#else
-  return (((uint64_t)hi<<32)|lo);
-#endif
-}
-#endif
-
-#ifdef ALLBITS
 /* generates a random number on (0,1)-real-interval */
-static double
+static float
 randu32 (void)
 {
-  return ((double)randi32() + 0.5) * (1.0/4294967296.0);
+  return ((float)randi32() + 0.5) * (1.0/4294967296.0);
   /* divided by 2^32 */
 }
-#else
+
 /* generates a random number on (0,1) with 53-bit resolution */
 static double
 randu53 (void)
@@ -414,35 +404,22 @@
   const uint32_t b=randi32()>>6;
   return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
 }
-#endif
 
 /* Determine mantissa for uniform doubles */
 double
 oct_randu (void)
 {
-#ifdef ALLBITS
+  return randu53 ();
+}
+
+/* Determine mantissa for uniform floats */
+float
+oct_float_randu (void)
+{
   return randu32 ();
-#else
-  return randu53 ();
-#endif
 }
 
 /* ===== Ziggurat normal and exponential generators ===== */
-#ifdef ALLBITS
-# define ZIGINT uint32_t
-# define EMANTISSA 4294967296.0 /* 32 bit mantissa */
-# define ERANDI randi32() /* 32 bits for mantissa */
-# define NMANTISSA 2147483648.0 /* 31 bit mantissa */
-# define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
-# define RANDU randu32()
-#else
-# define ZIGINT uint64_t
-# define EMANTISSA 9007199254740992.0  /* 53 bit mantissa */
-# define ERANDI randi53() /* 53 bits for mantissa */
-# define NMANTISSA EMANTISSA
-# define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */
-# define RANDU randu53()
-#endif
 
 #define ZIGGURAT_TABLE_SIZE 256
 
@@ -454,6 +431,13 @@
 #define ZIGGURAT_EXP_INV_R 0.129918765548341586
 #define EXP_SECTION_AREA 0.0039496598225815571993
 
+#define ZIGINT uint64_t
+#define EMANTISSA 9007199254740992.0  /* 53 bit mantissa */
+#define ERANDI randi53() /* 53 bits for mantissa */
+#define NMANTISSA EMANTISSA
+#define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */
+#define RANDU randu53()
+
 static ZIGINT ki[ZIGGURAT_TABLE_SIZE];
 static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE];
 static ZIGINT ke[ZIGGURAT_TABLE_SIZE];
@@ -592,7 +576,6 @@
        * Of course, different compilers and operating systems may
        * have something to do with this.
        */
-#if !defined(ALLBITS)
 # if HAVE_X86_32
       /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */
       double x;
@@ -615,14 +598,6 @@
       const double x = ( r&1 ? -rabs : rabs) * wi[idx];
 # endif /* !HAVE_X86_32 */
       if (rabs < (int64_t)ki[idx])
-#else /* ALLBITS */
-      /* 32-bit mantissa */
-      const uint32_t r = randi32();
-      const uint32_t rabs = r&LMASK;
-      const int idx = (int)(r&0xFF);
-      const double x = ((int32_t)r) * wi[idx];
-      if (rabs < ki[idx])
-#endif /* ALLBITS */
         return x;        /* 99.3% of the time we return here 1st try */
       else if (idx == 0)
         {
@@ -677,6 +652,173 @@
     }
 }
 
+#undef ZIGINT
+#undef EMANTISSA
+#undef ERANDI
+#undef NMANTISSA
+#undef NRANDI
+#undef RANDU
+
+#define ZIGINT uint32_t
+#define EMANTISSA 4294967296.0 /* 32 bit mantissa */
+#define ERANDI randi32() /* 32 bits for mantissa */
+#define NMANTISSA 2147483648.0 /* 31 bit mantissa */
+#define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
+#define RANDU randu32()
+
+static ZIGINT fki[ZIGGURAT_TABLE_SIZE];
+static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE];
+static ZIGINT fke[ZIGGURAT_TABLE_SIZE];
+static float fwe[ZIGGURAT_TABLE_SIZE], ffe[ZIGGURAT_TABLE_SIZE];
+
+
+static void
+create_ziggurat_float_tables (void)
+{
+  int i;
+  float x, x1;
+
+  /* Ziggurat tables for the normal distribution */
+  x1 = ZIGGURAT_NOR_R;
+  fwi[255] = x1 / NMANTISSA;
+  ffi[255] = exp (-0.5 * x1 * x1);
+
+  /* Index zero is special for tail strip, where Marsaglia and Tsang
+   * defines this as
+   * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
+   * where v is the area of each strip of the ziggurat.
+   */
+  fki[0] = (ZIGINT) (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA);
+  fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA;
+  ffi[0] = 1.;
+
+  for (i = 254; i > 0; i--)
+    {
+      /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
+       * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
+       */
+      x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + ffi[i+1]));
+      fki[i+1] = (ZIGINT)(x / x1 * NMANTISSA);
+      fwi[i] = x / NMANTISSA;
+      ffi[i] = exp (-0.5 * x * x);
+      x1 = x;
+    }
+
+  fki[1] = 0;
+
+  /* Zigurrat tables for the exponential distribution */
+  x1 = ZIGGURAT_EXP_R;
+  fwe[255] = x1 / EMANTISSA;
+  ffe[255] = exp (-x1);
+
+  /* Index zero is special for tail strip, where Marsaglia and Tsang
+   * defines this as
+   * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
+   * where v is the area of each strip of the ziggurat.
+   */
+  fke[0] = (ZIGINT) (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA);
+  fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA;
+  ffe[0] = 1.;
+
+  for (i = 254; i > 0; i--)
+    {
+      /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
+       * need inverse operator of y = exp(-x) -> x = -ln(y)
+       */
+      x = - log(EXP_SECTION_AREA / x1 + ffe[i+1]);
+      fke[i+1] = (ZIGINT)(x / x1 * EMANTISSA);
+      fwe[i] = x / EMANTISSA;
+      ffe[i] = exp (-x);
+      x1 = x;
+    }
+  fke[1] = 0;
+
+  inittf = 0;
+}
+
+/*
+ * Here is the guts of the algorithm. As Marsaglia and Tsang state the
+ * algorithm in their paper
+ *
+ * 1) Calculate a random signed integer j and let i be the index
+ *     provided by the rightmost 8-bits of j
+ * 2) Set x = j * w_i. If j < k_i return x
+ * 3) If i = 0, then return x from the tail
+ * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
+ * 5) goto step 1
+ *
+ * Where f is the functional form of the distribution, which for a normal
+ * distribution is exp(-0.5*x*x)
+ */
+
+float
+oct_float_randn (void)
+{
+  if (inittf)
+    create_ziggurat_float_tables();
+
+  while (1)
+    {
+      /* 32-bit mantissa */
+      const uint32_t r = randi32();
+      const uint32_t rabs = r&LMASK;
+      const int idx = (int)(r&0xFF);
+      const float x = ((int32_t)r) * fwi[idx];
+      if (rabs < fki[idx])
+        return x;        /* 99.3% of the time we return here 1st try */
+      else if (idx == 0)
+        {
+          /* As stated in Marsaglia and Tsang
+           *
+           * For the normal tail, the method of Marsaglia[5] provides:
+           * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
+           * then return r+x. Except that r+x is always in the positive
+           * tail!!!! Any thing random might be used to determine the
+           * sign, but as we already have r we might as well use it
+           *
+           * [PAK] but not the bottom 8 bits, since they are all 0 here!
+           */
+          float xx, yy;
+          do
+            {
+              xx = - ZIGGURAT_NOR_INV_R * log (RANDU);
+              yy = - log (RANDU);
+            }
+          while ( yy+yy <= xx*xx);
+          return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
+        }
+      else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp(-0.5*x*x))
+        return x;
+    }
+}
+
+float
+oct_float_rande (void)
+{
+  if (inittf)
+    create_ziggurat_float_tables();
+
+  while (1)
+    {
+      ZIGINT ri = ERANDI;
+      const int idx = (int)(ri & 0xFF);
+      const float x = ri * fwe[idx];
+      if (ri < fke[idx])
+        return x;               // 98.9% of the time we return here 1st try
+      else if (idx == 0)
+        {
+          /* As stated in Marsaglia and Tsang
+           *
+           * For the exponential tail, the method of Marsaglia[5] provides:
+           * x = r - ln(U);
+           */
+          return ZIGGURAT_EXP_R - log(RANDU);
+        }
+      else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp(-x))
+        return x;
+    }
+}
+
 /* Array generators */
 void
 oct_fill_randu (octave_idx_type n, double *p)
@@ -701,3 +843,27 @@
   for (i = 0; i < n; i++)
     p[i] = oct_rande();
 }
+
+void
+oct_fill_float_randu (octave_idx_type n, float *p)
+{
+  octave_idx_type i;
+  for (i = 0; i < n; i++)
+    p[i] = oct_float_randu();
+}
+
+void
+oct_fill_float_randn (octave_idx_type n, float *p)
+{
+  octave_idx_type i;
+  for (i = 0; i < n; i++)
+    p[i] = oct_float_randn();
+}
+
+void
+oct_fill_float_rande (octave_idx_type n, float *p)
+{
+  octave_idx_type i;
+  for (i = 0; i < n; i++)
+    p[i] = oct_float_rande();
+}
--- a/liboctave/randmtzig.h	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/randmtzig.h	Sat May 19 00:14:59 2012 +0200
@@ -82,11 +82,19 @@
 extern OCTAVE_API double oct_randn (void);
 extern OCTAVE_API double oct_rande (void);
 
+extern OCTAVE_API float oct_float_randu (void);
+extern OCTAVE_API float oct_float_randn (void);
+extern OCTAVE_API float oct_float_rande (void);
+
 /* === Array generators === */
 extern OCTAVE_API void oct_fill_randu (octave_idx_type n, double *p);
 extern OCTAVE_API void oct_fill_randn (octave_idx_type n, double *p);
 extern OCTAVE_API void oct_fill_rande (octave_idx_type n, double *p);
 
+extern OCTAVE_API void oct_fill_float_randu (octave_idx_type n, float *p);
+extern OCTAVE_API void oct_fill_float_randn (octave_idx_type n, float *p);
+extern OCTAVE_API void oct_fill_float_rande (octave_idx_type n, float *p);
+
 #ifdef  __cplusplus
 }
 #endif
--- a/liboctave/randpoisson.c	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/randpoisson.c	Sat May 19 00:14:59 2012 +0200
@@ -368,6 +368,46 @@
   }
 }
 
+static void
+poisson_cdf_lookup_float(double lambda, float *p, size_t n)
+{
+  double t[TABLESIZE];
+
+  /* Precompute the table for the u up to and including 0.458.
+   * We will almost certainly need it. */
+  int intlambda = (int)floor(lambda);
+  double P;
+  int tableidx;
+  size_t i = n;
+
+  t[0] = P = exp(-lambda);
+  for (tableidx = 1; tableidx <= intlambda; tableidx++) {
+    P = P*lambda/(double)tableidx;
+    t[tableidx] = t[tableidx-1] + P;
+  }
+
+  while (i-- > 0) {
+    double u = RUNI;
+    int k = (u > 0.458 ? intlambda : 0);
+  nextk:
+    if ( u <= t[k] ) {
+      p[i] = (float) k;
+      continue;
+    }
+    if (++k < tableidx) goto nextk;
+
+    while (tableidx < TABLESIZE) {
+      P = P*lambda/(double)tableidx;
+      t[tableidx] = t[tableidx-1] + P;
+      if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
+      tableidx++;
+      if (u <= t[tableidx-1]) break;
+    }
+
+    p[i] = (float)(tableidx-1);
+  }
+}
+
 /* From Press, et al., Numerical Recipes */
 static void
 poisson_rejection (double lambda, double *p, size_t n)
@@ -392,6 +432,30 @@
     }
 }
 
+/* From Press, et al., Numerical Recipes */
+static void
+poisson_rejection_float (double lambda, float *p, size_t n)
+{
+  double sq = sqrt(2.0*lambda);
+  double alxm = log(lambda);
+  double g = lambda*alxm - LGAMMA(lambda+1.0);
+  size_t i;
+
+  for (i = 0; i < n; i++)
+    {
+      double y, em, t;
+      do {
+        do {
+          y = tan(M_PI*RUNI);
+          em = sq * y + lambda;
+        } while (em < 0.0);
+        em = floor(em);
+        t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g);
+      } while (RUNI > t);
+      p[i] = em;
+    }
+}
+
 /* The cutoff of L <= 1e8 in the following two functions before using
  * the normal approximation is based on:
  *   > L=1e8; x=floor(linspace(0,2*L,1000));
@@ -463,3 +527,68 @@
   }
   return ret;
 }
+
+/* Generate a set of poisson numbers with the same distribution */
+void
+oct_fill_float_randp (float FL, octave_idx_type n, float *p)
+{
+  double L = FL;
+  octave_idx_type i;
+  if (L < 0.0 || INFINITE(L))
+    {
+      for (i=0; i<n; i++)
+        p[i] = NAN;
+    }
+  else if (L <= 10.0)
+    {
+      poisson_cdf_lookup_float(L, p, n);
+    }
+  else if (L <= 1e8)
+    {
+      for (i=0; i<n; i++)
+        p[i] = pprsc(L);
+    }
+  else
+    {
+      /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
+      const double sqrtL = sqrt(L);
+      for (i = 0; i < n; i++)
+        {
+          p[i] = floor(RNOR*sqrtL + L + 0.5);
+          if (p[i] < 0.0)
+            p[i] = 0.0; /* will probably never happen */
+        }
+    }
+}
+
+/* Generate one poisson variate */
+float
+oct_float_randp (float FL)
+{
+  double L = FL;
+  float ret;
+  if (L < 0.0) ret = NAN;
+  else if (L <= 12.0) {
+    /* From Press, et al. Numerical recipes */
+    double g = exp(-L);
+    int em = -1;
+    double t = 1.0;
+    do {
+      ++em;
+      t *= RUNI;
+    } while (t > g);
+    ret = em;
+  } else if (L <= 1e8) {
+    /* numerical recipes */
+    poisson_rejection_float(L, &ret, 1);
+  } else if (INFINITE(L)) {
+    /* FIXME R uses NaN, but the normal approx. suggests that as
+     * limit should be inf. Which is correct? */
+    ret = NAN;
+  } else {
+    /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
+    ret = floor(RNOR*sqrt(L) + L + 0.5);
+    if (ret < 0.0) ret = 0.0; /* will probably never happen */
+  }
+  return ret;
+}
--- a/liboctave/randpoisson.h	Fri May 18 13:16:04 2012 -0700
+++ b/liboctave/randpoisson.h	Sat May 19 00:14:59 2012 +0200
@@ -32,6 +32,9 @@
 extern OCTAVE_API double oct_randp (double L);
 extern OCTAVE_API void oct_fill_randp (double L, octave_idx_type n, double *p);
 
+extern OCTAVE_API float oct_float_randp (float L);
+extern OCTAVE_API void oct_fill_float_randp (float L, octave_idx_type n, float *p);
+
 #ifdef  __cplusplus
 }
 #endif
--- a/src/DLD-FUNCTIONS/rand.cc	Fri May 18 13:16:04 2012 -0700
+++ b/src/DLD-FUNCTIONS/rand.cc	Sat May 19 00:14:59 2012 +0200
@@ -1,3 +1,4 @@
+
 /*
 
 Copyright (C) 1996-2012 John W. Eaton
@@ -60,6 +61,7 @@
   NDArray a;
   int idx = 0;
   dim_vector dims;
+  bool is_single = false;
 
   unwind_protect frame;
   // Restore current distribution on any exit.
@@ -68,6 +70,19 @@
 
   octave_rand::distribution (distribution);
 
+  if (nargin > 0 && args(nargin-1).is_string())
+    {
+      std::string s_arg = args(nargin-1).string_value ();
+
+      if (s_arg == "single")
+        {
+          is_single = true;
+          nargin--;
+        }
+      else if (s_arg == "double")
+        nargin--;
+    }
+
   if (additional_arg)
     {
       if (nargin == 0)
@@ -298,27 +313,54 @@
 
   dims.chop_trailing_singletons ();
 
-  if (additional_arg)
+  if (is_single)
     {
-      if (a.length() == 1)
-        return octave_rand::nd_array (dims, a(0));
-      else
+      if (additional_arg)
         {
-          if (a.dims() != dims)
+          if (a.length() == 1)
+            return octave_rand::float_nd_array (dims, a(0));
+          else
             {
-              error ("%s: mismatch in argument size", fcn);
-              return retval;
+              if (a.dims() != dims)
+                {
+                  error ("%s: mismatch in argument size", fcn);
+                  return retval;
+                }
+              octave_idx_type len = a.length ();
+              FloatNDArray m (dims);
+              float *v = m.fortran_vec ();
+              for (octave_idx_type i = 0; i < len; i++)
+                v[i] = octave_rand::float_scalar (a(i));
+              return m;
             }
-          octave_idx_type len = a.length ();
-          NDArray m (dims);
-          double *v = m.fortran_vec ();
-          for (octave_idx_type i = 0; i < len; i++)
-            v[i] = octave_rand::scalar (a(i));
-          return m;
         }
+      else
+        return octave_rand::float_nd_array (dims);
     }
   else
-    return octave_rand::nd_array (dims);
+    {
+      if (additional_arg)
+        {
+          if (a.length() == 1)
+            return octave_rand::nd_array (dims, a(0));
+          else
+            {
+              if (a.dims() != dims)
+                {
+                  error ("%s: mismatch in argument size", fcn);
+                  return retval;
+                }
+              octave_idx_type len = a.length ();
+              NDArray m (dims);
+              double *v = m.fortran_vec ();
+              for (octave_idx_type i = 0; i < len; i++)
+                v[i] = octave_rand::scalar (a(i));
+              return m;
+            }
+        }
+      else
+        return octave_rand::nd_array (dims);
+    }
 }
 
 DEFUN_DLD (rand, args, ,
@@ -332,6 +374,8 @@
 @deftypefnx {Loadable Function} {@var{v} =} rand (\"seed\")\n\
 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{v})\n\
 @deftypefnx {Loadable Function} {} rand (\"seed\", \"reset\")\n\
+@deftypefnx {Loadable Function} {} rand (@dots{}, \"single\")\n\
+@deftypefnx {Loadable Function} {} rand (@dots{}, \"double\")\n\
 Return a matrix with random elements uniformly distributed on the\n\
 interval (0, 1).  The arguments are handled the same as the arguments\n\
 for @code{eye}.\n\
@@ -402,6 +446,9 @@
 \n\
 The state or seed of the generator can be reset to a new random value\n\
 using the \"reset\" keyword.\n\
+\n\
+The class of the value returned can be controlled by a trailing \"double\"\n\
+or \"single\" argument. These are the only valid classes.\n\
 @seealso{randn, rande, randg, randp}\n\
 @end deftypefn")
 {
@@ -499,6 +546,8 @@
 @deftypefnx {Loadable Function} {@var{v} =} randn (\"seed\")\n\
 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{v})\n\
 @deftypefnx {Loadable Function} {} randn (\"seed\", \"reset\")\n\
+@deftypefnx {Loadable Function} {} randn (@dots{}, \"single\")\n\
+@deftypefnx {Loadable Function} {} randn (@dots{}, \"double\")\n\
 Return a matrix with normally distributed random\n\
 elements having zero mean and variance one.  The arguments are\n\
 handled the same as the arguments for @code{rand}.\n\
@@ -506,6 +555,9 @@
 By default, @code{randn} uses the Marsaglia and Tsang ``Ziggurat technique''\n\
 to transform from a uniform to a normal distribution.\n\
 \n\
+The class of the value returned can be controlled by a trailing \"double\"\n\
+or \"single\" argument. These are the only valid classes.\n\
+\n\
 Reference: G. Marsaglia and W.W. Tsang,\n\
 @cite{Ziggurat Method for Generating Random Variables},\n\
 J. Statistical Software, vol 5, 2000,\n\
@@ -565,12 +617,17 @@
 @deftypefnx {Loadable Function} {@var{v} =} rande (\"seed\")\n\
 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{v})\n\
 @deftypefnx {Loadable Function} {} rande (\"seed\", \"reset\")\n\
+@deftypefnx {Loadable Function} {} rande (@dots{}, \"single\")\n\
+@deftypefnx {Loadable Function} {} rande (@dots{}, \"double\")\n\
 Return a matrix with exponentially distributed random elements.  The\n\
 arguments are handled the same as the arguments for @code{rand}.\n\
 \n\
 By default, @code{randn} uses the Marsaglia and Tsang ``Ziggurat technique''\n\
 to transform from a uniform to an exponential distribution.\n\
 \n\
+The class of the value returned can be controlled by a trailing \"double\"\n\
+or \"single\" argument. These are the only valid classes.\n\
+\n\
 Reference: G. Marsaglia and W.W. Tsang,\n\
 @cite{Ziggurat Method for Generating Random Variables},\n\
 J. Statistical Software, vol 5, 2000,\n\
@@ -632,6 +689,8 @@
 @deftypefnx {Loadable Function} {@var{v} =} randg (\"seed\")\n\
 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{v})\n\
 @deftypefnx {Loadable Function} {} randg (\"seed\", \"reset\")\n\
+@deftypefnx {Loadable Function} {} randg (@dots{}, \"single\")\n\
+@deftypefnx {Loadable Function} {} randg (@dots{}, \"double\")\n\
 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\
 The arguments are handled the same as the arguments for @code{rand},\n\
 except for the argument @var{a}.\n\
@@ -711,6 +770,9 @@
 @end example\n\
 \n\
 @end table\n\
+\n\
+The class of the value returned can be controlled by a trailing \"double\"\n\
+or \"single\" argument. These are the only valid classes.\n\
 @seealso{rand, randn, rande, randp}\n\
 @end deftypefn")
 {
@@ -898,6 +960,8 @@
 @deftypefnx {Loadable Function} {@var{v} =} randp (\"seed\")\n\
 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{v})\n\
 @deftypefnx {Loadable Function} {} randp (\"seed\", \"reset\")\n\
+@deftypefnx {Loadable Function} {} randp (@dots{}, \"single\")\n\
+@deftypefnx {Loadable Function} {} randp (@dots{}, \"double\")\n\
 Return a matrix with Poisson distributed random elements with mean value\n\
 parameter given by the first argument, @var{l}.  The arguments\n\
 are handled the same as the arguments for @code{rand}, except for the\n\
@@ -928,6 +992,9 @@
 L. Montanet, et al., @cite{Review of Particle Properties}, Physical Review\n\
 D 50 p1284, 1994.\n\
 @end table\n\
+\n\
+The class of the value returned can be controlled by a trailing \"double\"\n\
+or \"single\" argument. These are the only valid classes.\n\
 @seealso{rand, randn, rande, randg}\n\
 @end deftypefn")
 {