--- a/liboctave/ChangeLog	Tue Feb 26 20:29:58 2008 -0500
+++ b/liboctave/ChangeLog	Wed Feb 27 04:33:39 2008 -0500
@@ -1,3 +1,7 @@
+2008-02-27  John W. Eaton  <jwe@octave.org>
+
+	* oct-rand.cc (class octave_rand): Make it a proper singleton class.
+
 2008-02-26  John W. Eaton  <jwe@octave.org>
 
 	* oct-rand.cc (get_dist_id): Fix typo.

--- a/liboctave/oct-rand.cc	Tue Feb 26 20:29:58 2008 -0500
+++ b/liboctave/oct-rand.cc	Wed Feb 27 04:33:39 2008 -0500
@@ -39,23 +39,6 @@
 #include "randgamma.h"
 #include "mach-info.h"
 
-static const int unknown_dist = 0;
-static const int uniform_dist = 1;
-static const int normal_dist = 2;
-static const int expon_dist = 3;
-static const int poisson_dist = 4;
-static const int gamma_dist = 5;
-
-// Current distribution of random numbers.
-static int current_distribution = uniform_dist;
-
-// Has the seed/state been set yet?
-static bool old_initialized = false;
-static bool new_initialized = false;
-static bool use_old_generators = false;
-
-std::map<int, ColumnVector> rand_states;
-
 extern "C"
 {
   F77_RET_T
@@ -86,160 +69,39 @@
   F77_FUNC (setcgn, SETCGN) (const octave_idx_type&);
 }
 
-static octave_idx_type
-force_to_fit_range (octave_idx_type i, octave_idx_type lo, octave_idx_type hi)
-{
-  assert (hi > lo && lo >= 0 && hi > lo);
-
-  i = i > 0 ? i : -i;
-
-  if (i < lo)
-    i = lo;
-  else if (i > hi)
-    i = i % hi;
-
-  return i;
-}
-
-// 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
-// work ok to give fairly different seeds each time Octave starts.
-
-static void
-do_old_initialization (void)
-{
-  octave_localtime tm;
-  int stored_distribution = current_distribution;
-  F77_FUNC (setcgn, SETCGN) (uniform_dist);
+octave_rand *octave_rand::instance = 0;
 
-  int hour = tm.hour() + 1;
-  int minute = tm.min() + 1;
-  int second = tm.sec() + 1;
-
-  octave_idx_type s0 = tm.mday() * hour * minute * second;
-  octave_idx_type s1 = hour * minute * second;
-
-  s0 = force_to_fit_range (s0, 1, 2147483563);
-  s1 = force_to_fit_range (s1, 1, 2147483399);
-
-  F77_FUNC (setall, SETALL) (s0, s1);
-  F77_FUNC (setcgn, SETCGN) (stored_distribution);
+octave_rand::octave_rand (void)
+  : current_distribution (uniform_dist), use_old_generators (false),
+    rand_states ()
+{
+  initialize_ranlib_generators ();
 
-  old_initialized = true;
-}
-
-static ColumnVector
-get_internal_state (void)
-{
-  ColumnVector s (MT_N + 1);
-
-  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
-
-  oct_get_state (tmp);
-
-  for (octave_idx_type i = 0; i <= MT_N; i++)
-    s.elem (i) = static_cast<double> (tmp [i]);
-
-  return s;
+  initialize_mersenne_twister ();
 }
 
-static inline void
-save_state (void)
-{
-  rand_states[current_distribution] = get_internal_state ();;
-}
-
-static void
-initialize_rand_states (void)
-{
-  if (! new_initialized)
-    {
-      oct_init_by_entropy ();
-
-      ColumnVector s = get_internal_state ();
-
-      rand_states[uniform_dist] = s;
-      rand_states[normal_dist] = s;
-      rand_states[expon_dist] = s;
-      rand_states[poisson_dist] = s;
-      rand_states[gamma_dist] = s;
-
-      new_initialized = true;
-    }
-}
-
-static inline void
-maybe_initialize (void)
+bool
+octave_rand::instance_ok (void)
 {
-  if (use_old_generators)
-    {
-      if (! old_initialized)
-	do_old_initialization ();
-    }
-  else
-    {
-      if (! new_initialized)
-	initialize_rand_states ();
-    }
-}
+  bool retval = true;
+
+  if (! instance)
+    instance = new octave_rand ();
 
-static int
-get_dist_id (const std::string& d)
-{
-  int retval = unknown_dist;
+  if (! instance)
+    {
+      (*current_liboctave_error_handler)
+	("unable to create octave_rand object!");
 
-  if (d == "uniform" || d == "rand")
-    retval = uniform_dist;
-  else if (d == "normal" || d == "randn")
-    retval = normal_dist;
-  else if (d == "exponential" || d == "rande")
-    retval = expon_dist;
-  else if (d == "poisson" || d == "randp")
-    retval = poisson_dist;
-  else if (d == "gamma" || d == "randg")
-    retval = gamma_dist;
-  else
-    (*current_liboctave_error_handler)
-      ("rand: invalid distribution `%s'", d.c_str ());
+      retval = false;
+    }
 
   return retval;
 }
 
-static void
-set_internal_state (const ColumnVector& s)
+double
+octave_rand::do_seed (void)
 {
-  octave_idx_type len = s.length ();
-  octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
-
-  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
-
-  for (octave_idx_type i = 0; i < n; i++)
-    tmp[i] = static_cast<uint32_t> (s.elem(i));
-
-  if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0)
-    oct_set_state (tmp);
-  else
-    oct_init_by_array (tmp, len);
-}
-
-static inline void
-switch_to_generator (int dist)
-{
-  if (dist != current_distribution)
-    {
-      current_distribution = dist;
-
-      set_internal_state (rand_states[dist]);
-    }
-}
-
-double
-octave_rand::seed (void)
-{
-  if (! old_initialized)
-    do_old_initialization ();
-
   union d2i { double d; octave_idx_type i[2]; };
   union d2i u;
     
@@ -258,13 +120,26 @@
   return u.d;
 }
 
+static octave_idx_type
+force_to_fit_range (octave_idx_type i, octave_idx_type lo, octave_idx_type hi)
+{
+  assert (hi > lo && lo >= 0 && hi > lo);
+
+  i = i > 0 ? i : -i;
+
+  if (i < lo)
+    i = lo;
+  else if (i > hi)
+    i = i % hi;
+
+  return i;
+}
+
 void
-octave_rand::seed (double s)
+octave_rand::do_seed (double s)
 {
   use_old_generators = true;
 
-  maybe_initialize ();
-
   int i0, i1;
   union d2i { double d; octave_idx_type i[2]; };
   union d2i u;
@@ -288,21 +163,16 @@
 }
 
 ColumnVector
-octave_rand::state (const std::string& d)
+octave_rand::do_state (const std::string& d)
 {
-  if (! new_initialized)
-    initialize_rand_states ();
-
   return rand_states[d.empty () ? current_distribution : get_dist_id (d)];
 }
 
 void
-octave_rand::state (const ColumnVector& s, const std::string& d)
+octave_rand::do_state (const ColumnVector& s, const std::string& d)
 {
   use_old_generators = false;
 
-  maybe_initialize ();
-
   int old_dist = current_distribution;
 
   int new_dist = d.empty () ? current_distribution : get_dist_id (d);
@@ -321,12 +191,10 @@
 }
 
 std::string
-octave_rand::distribution (void)
+octave_rand::do_distribution (void)
 {
   std::string retval;
 
-  maybe_initialize ();
-
   switch (current_distribution)
     {
     case uniform_dist:
@@ -359,7 +227,7 @@
 }
 
 void
-octave_rand::distribution (const std::string& d)
+octave_rand::do_distribution (const std::string& d)
 {
   int id = get_dist_id (d);
 
@@ -393,50 +261,40 @@
 }
 
 void
-octave_rand::uniform_distribution (void)
+octave_rand::do_uniform_distribution (void)
 {
-  maybe_initialize ();
-
   switch_to_generator (uniform_dist);
 
   F77_FUNC (setcgn, SETCGN) (uniform_dist);
 }
 
 void
-octave_rand::normal_distribution (void)
+octave_rand::do_normal_distribution (void)
 {
-  maybe_initialize ();
-
   switch_to_generator (normal_dist);
 
   F77_FUNC (setcgn, SETCGN) (normal_dist);
 }
 
 void
-octave_rand::exponential_distribution (void)
+octave_rand::do_exponential_distribution (void)
 {
-  maybe_initialize ();
-
   switch_to_generator (expon_dist);
 
   F77_FUNC (setcgn, SETCGN) (expon_dist);
 }
 
 void
-octave_rand::poisson_distribution (void)
+octave_rand::do_poisson_distribution (void)
 {
-  maybe_initialize ();
-
   switch_to_generator (poisson_dist);
 
   F77_FUNC (setcgn, SETCGN) (poisson_dist);
 }
 
 void
-octave_rand::gamma_distribution (void)
+octave_rand::do_gamma_distribution (void)
 {
-  maybe_initialize ();
-
   switch_to_generator (gamma_dist);
 
   F77_FUNC (setcgn, SETCGN) (gamma_dist);
@@ -444,10 +302,8 @@
 
 
 double
-octave_rand::scalar (double a)
+octave_rand::do_scalar (double a)
 {
-  maybe_initialize ();
-
   double retval = 0.0;
 
   if (use_old_generators)
@@ -526,6 +382,167 @@
   return retval;
 }
 
+Matrix
+octave_rand::do_matrix (octave_idx_type n, octave_idx_type m, double a)
+{
+  Matrix retval;
+
+  if (n >= 0 && m >= 0)
+    {
+      retval.resize (n, m);
+
+      if (n > 0 && m > 0)
+	fill (retval.capacity(), retval.fortran_vec(), a);
+    }
+  else
+    (*current_liboctave_error_handler) ("rand: invalid negative argument");
+
+  return retval;
+}
+
+NDArray
+octave_rand::do_nd_array (const dim_vector& dims, double a)
+{
+  NDArray retval;
+
+  if (! dims.all_zero ())
+    {
+      retval.resize (dims);
+
+      fill (retval.capacity(), retval.fortran_vec(), a);
+    }
+
+  return retval;
+}
+
+Array<double>
+octave_rand::do_vector (octave_idx_type n, double a)
+{
+  Array<double> retval;
+
+  if (n > 0)
+    {
+      retval.resize (n);
+
+      fill (retval.capacity (), retval.fortran_vec (), a);
+    }
+  else if (n < 0)
+    (*current_liboctave_error_handler) ("rand: invalid negative argument");
+
+  return retval;
+}
+
+// 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
+// work ok to give fairly different seeds each time Octave starts.
+
+void
+octave_rand::initialize_ranlib_generators (void)
+{
+  octave_localtime tm;
+  int stored_distribution = current_distribution;
+  F77_FUNC (setcgn, SETCGN) (uniform_dist);
+
+  int hour = tm.hour() + 1;
+  int minute = tm.min() + 1;
+  int second = tm.sec() + 1;
+
+  octave_idx_type s0 = tm.mday() * hour * minute * second;
+  octave_idx_type s1 = hour * minute * second;
+
+  s0 = force_to_fit_range (s0, 1, 2147483563);
+  s1 = force_to_fit_range (s1, 1, 2147483399);
+
+  F77_FUNC (setall, SETALL) (s0, s1);
+  F77_FUNC (setcgn, SETCGN) (stored_distribution);
+}
+
+void
+octave_rand::initialize_mersenne_twister (void)
+{
+  oct_init_by_entropy ();
+
+  ColumnVector s = get_internal_state ();
+
+  rand_states[uniform_dist] = s;
+  rand_states[normal_dist] = s;
+  rand_states[expon_dist] = s;
+  rand_states[poisson_dist] = s;
+  rand_states[gamma_dist] = s;
+}
+
+ColumnVector
+octave_rand::get_internal_state (void)
+{
+  ColumnVector s (MT_N + 1);
+
+  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
+
+  oct_get_state (tmp);
+
+  for (octave_idx_type i = 0; i <= MT_N; i++)
+    s.elem (i) = static_cast<double> (tmp [i]);
+
+  return s;
+}
+
+void
+octave_rand::save_state (void)
+{
+  rand_states[current_distribution] = get_internal_state ();;
+}
+
+int
+octave_rand::get_dist_id (const std::string& d)
+{
+  int retval = unknown_dist;
+
+  if (d == "uniform" || d == "rand")
+    retval = uniform_dist;
+  else if (d == "normal" || d == "randn")
+    retval = normal_dist;
+  else if (d == "exponential" || d == "rande")
+    retval = expon_dist;
+  else if (d == "poisson" || d == "randp")
+    retval = poisson_dist;
+  else if (d == "gamma" || d == "randg")
+    retval = gamma_dist;
+  else
+    (*current_liboctave_error_handler)
+      ("rand: invalid distribution `%s'", d.c_str ());
+
+  return retval;
+}
+
+void
+octave_rand::set_internal_state (const ColumnVector& s)
+{
+  octave_idx_type len = s.length ();
+  octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
+
+  OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
+
+  for (octave_idx_type i = 0; i < n; i++)
+    tmp[i] = static_cast<uint32_t> (s.elem(i));
+
+  if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0)
+    oct_set_state (tmp);
+  else
+    oct_init_by_array (tmp, len);
+}
+
+void
+octave_rand::switch_to_generator (int dist)
+{
+  if (dist != current_distribution)
+    {
+      current_distribution = dist;
+
+      set_internal_state (rand_states[dist]);
+    }
+}
+
 #define MAKE_RAND(len) \
   do \
     { \
@@ -539,11 +556,9 @@
     } \
   while (0)
 
-static void
-fill_rand (octave_idx_type len, double *v, double a)
+void
+octave_rand::fill (octave_idx_type len, double *v, double a)
 {
-  maybe_initialize ();
-
   if (len < 1)
     return;
 
@@ -630,58 +645,6 @@
   return;
 }
 
-Matrix
-octave_rand::matrix (octave_idx_type n, octave_idx_type m, double a)
-{
-  Matrix retval;
-
-  if (n >= 0 && m >= 0)
-    {
-      retval.resize (n, m);
-
-      if (n > 0 && m > 0)
-	fill_rand (retval.capacity(), retval.fortran_vec(), a);
-    }
-  else
-    (*current_liboctave_error_handler) ("rand: invalid negative argument");
-
-  return retval;
-}
-
-NDArray
-octave_rand::nd_array (const dim_vector& dims, double a)
-{
-  NDArray retval;
-
-  if (! dims.all_zero ())
-    {
-      retval.resize (dims);
-
-      fill_rand (retval.capacity(), retval.fortran_vec(), a);
-    }
-
-  return retval;
-}
-
-Array<double>
-octave_rand::vector (octave_idx_type n, double a)
-{
-  maybe_initialize ();
-
-  Array<double> retval;
-
-  if (n > 0)
-    {
-      retval.resize (n);
-
-      fill_rand (retval.capacity(), retval.fortran_vec(), a);
-    }
-  else if (n < 0)
-    (*current_liboctave_error_handler) ("rand: invalid negative argument");
-
-  return retval;
-}
-
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***

--- a/liboctave/oct-rand.h	Tue Feb 26 20:29:58 2008 -0500
+++ b/liboctave/oct-rand.h	Wed Feb 27 04:33:39 2008 -0500
@@ -23,59 +23,209 @@
 #if !defined (octave_rand_h)
 #define octave_rand_h 1
 
+#include <map>
 #include <string>
 
 #include "dColVector.h"
 #include "dMatrix.h"
 #include "dNDArray.h"
+#include "lo-ieee.h"
 
-struct
+class
 OCTAVE_API
 octave_rand
 {
+protected:
+
+  octave_rand (void);
+
+public:
+
+  ~octave_rand (void) { }
+
+  static bool instance_ok (void);
+
   // Return the current seed.
-  static double seed (void);
+  static double seed (void)
+  {
+    return instance_ok () ? instance->do_seed () : octave_NaN;
+  }
 
   // Set the seed.
-  static void seed (double s);
+  static void seed (double s)
+  {
+    if (instance_ok ())
+      instance->do_seed (s);
+  }
 
   // Return the current state.
-  static ColumnVector state (const std::string& d = std::string ());
+  static ColumnVector state (const std::string& d = std::string ())
+  {
+    return instance_ok () ? instance->do_state (d) : ColumnVector ();
+  }
 
   // Set the current state/
   static void state (const ColumnVector &s,
-		     const std::string& d = std::string ());
+		     const std::string& d = std::string ())
+  {
+    if (instance_ok ())
+      instance->do_state (s, d);
+  }
   
   // Return the current distribution.
-  static std::string distribution (void);
+  static std::string distribution (void)
+  {
+    return instance_ok () ? instance->do_distribution () : std::string ();
+  }
 
   // Set the current distribution.  May be either "uniform" (the
   // default), "normal", "exponential", "poisson", or "gamma".
-  static void distribution (const std::string& d);
+  static void distribution (const std::string& d)
+  {
+    if (instance_ok ())
+      instance->do_distribution (d);
+  }
 
-  static void uniform_distribution (void);
+  static void uniform_distribution (void)
+  {
+    if (instance_ok ())
+      instance->do_uniform_distribution ();
+  }
 
-  static void normal_distribution (void);
+  static void normal_distribution (void)
+  {
+    if (instance_ok ())
+      instance->do_normal_distribution ();
+  }
 
-  static void exponential_distribution (void);
+  static void exponential_distribution (void)
+  {
+    if (instance_ok ())
+      instance->do_exponential_distribution ();
+  }
 
-  static void poisson_distribution (void);
+  static void poisson_distribution (void)
+  {
+    if (instance_ok ())
+      instance->do_poisson_distribution ();
+  }
 
-  static void gamma_distribution (void);
+  static void gamma_distribution (void)
+  {
+    if (instance_ok ())
+      instance->do_gamma_distribution ();
+  }
 
   // Return the next number from the sequence.
-  static double scalar (double a = 1.);
+  static double scalar (double a = 1.0)
+  {
+    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.);
+  static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.0)
+  {
+    return instance_ok () ? instance->do_matrix (r, c, a) : Matrix ();
+  }
 
   // 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.);
+  static NDArray nd_array (const dim_vector& dims, double a = 1.0)
+  {
+    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.);
+  static Array<double> vector (octave_idx_type n, double a = 1.0)
+  {
+    return instance_ok () ? instance->do_vector (n, a) : Array<double> ();
+  }
+
+private:
+
+  static octave_rand *instance;
+
+  enum
+  {
+    unknown_dist,
+    uniform_dist,
+    normal_dist,
+    expon_dist,
+    poisson_dist,
+    gamma_dist
+  };
+
+  // Current distribution of random numbers.
+  int current_distribution;
+
+  // If TRUE, use old RANLIB generators.  Otherwise, use Mersenne
+  // Twister generator.
+  bool use_old_generators;
+
+  // Saved MT states.
+  std::map<int, ColumnVector> rand_states;
+
+  // Return the current seed.
+  double do_seed (void);
+
+  // Set the seed.
+  void do_seed (double s);
+
+  // Return the current state.
+  ColumnVector do_state (const std::string& d);
+
+  // Set the current state/
+  void do_state (const ColumnVector &s, const std::string& d);
+  
+  // Return the current distribution.
+  std::string do_distribution (void);
+
+  // Set the current distribution.  May be either "uniform" (the
+  // default), "normal", "exponential", "poisson", or "gamma".
+  void do_distribution (const std::string& d);
+
+  void do_uniform_distribution (void);
+
+  void do_normal_distribution (void);
+
+  void do_exponential_distribution (void);
+
+  void do_poisson_distribution (void);
+
+  void do_gamma_distribution (void);
+
+  // 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 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.);
+
+  // Some helper functions.
+
+  void initialize_ranlib_generators (void);
+
+  void initialize_mersenne_twister (void);
+
+  ColumnVector get_internal_state (void);
+
+  void save_state (void);
+
+  int get_dist_id (const std::string& d);
+
+  void set_internal_state (const ColumnVector& s);
+
+  void switch_to_generator (int dist);
+
+  void fill (octave_idx_type len, double *v, double a);
 };
 
 #endif