--- a/doc/ChangeLog	Mon Apr 03 19:03:31 2006 +0000
+++ b/doc/ChangeLog	Mon Apr 03 19:18:26 2006 +0000
@@ -1,3 +1,8 @@
+2006-04-03  David Bateman  <dbateman@free.fr>
+
+	* interpreter/matrix.txi: Add rande, randp, randg and update
+          for different random generator behavior.
+
 2006-03-28  John W. Eaton  <jwe@octave.org>
 
 	* texinfo.tex: Update FSF address.

--- a/doc/interpreter/matrix.txi	Mon Apr 03 19:03:31 2006 +0000
+++ b/doc/interpreter/matrix.txi	Mon Apr 03 19:18:26 2006 +0000
@@ -143,8 +143,15 @@
 
 @DOCSTRING(randn)
 
-The @code{rand} and @code{randn} functions use separate generators.
-This ensures that
+@DOCSTRING(rande)
+
+@DOCSTRING(randp)
+
+@DOCSTRING(randg)
+
+The new random generators all use a common Mersenne Twister generator,
+and so the state of only one of the generators needs to be reset.
+The old generator function use separate generators. This ensures that
 
 @example
 @group

--- a/libcruft/ChangeLog	Mon Apr 03 19:03:31 2006 +0000
+++ b/libcruft/ChangeLog	Mon Apr 03 19:18:26 2006 +0000
@@ -1,3 +1,7 @@
+2006-04-03  David Bateman  <dbateman@free.fr>
+
+	* ranlib/wrap.f (dgenexp, dgengam, dignpoi): New functions.
+
 2006-03-21  John W. Eaton  <jwe@octave.org>
 
 	* misc/f77-fcn.h (F77_XFCN): Save octave_interrupt_immediately and

--- a/libcruft/ranlib/wrap.f	Mon Apr 03 19:03:31 2006 +0000
+++ b/libcruft/ranlib/wrap.f	Mon Apr 03 19:18:26 2006 +0000
@@ -8,3 +8,18 @@
       result = genunf (real (low), real (high))
       return
       end
+      subroutine dgenexp (av, result)
+      double precision av, result
+      result = genexp (real (av))
+      return
+      end
+      subroutine dgengam (a, r, result)
+      double precision a, r, result
+      result = gengam (real (a), real (r))
+      return
+      end
+      subroutine dignpoi (mu, result)
+      double precision mu, result
+      result = ignpoi (real (mu))
+      return
+      end

--- a/liboctave/ChangeLog	Mon Apr 03 19:03:31 2006 +0000
+++ b/liboctave/ChangeLog	Mon Apr 03 19:18:26 2006 +0000
@@ -1,3 +1,32 @@
+2006-04-03  David Bateman  <dbateman@free.fr>
+
+	* Makefile.in (INCLUDES): Add randgamma.h, randpoisson.h and
+        randmtzig.h to the list.
+	(LIBOCTAVE_C_SOURCES): Add randgamma.c, randpoisson.c and
+        randmtzig.c to the list.
+	* oct-rand.cc (do_old_initialization): Rename from do_initialization.
+	(use_old_generators): New variable.
+	(old_initialized): Rename from initialized.
+	(new_initialized): New variable.
+	(oct_init_by_entropy): New function.
+	(maybe_initialize): Initialize new or old generator depending on
+	value of use_old_generators.
+	(octave_rand::state): New functions.
+	(octave_rand::distribution): Add gamma, exponential and poisson
+	distributions.
+	(octave_rand::exponential_distribution,
+	octave_rand::poisson_distribution,
+	octave_rand::gamma_distribution): New methods to select
+        exponential, poisson or gamma distribution.
+	(octave_rand::scalar, octave_rand::matrix, octave_rand::nd_array,
+	octave_rand::vector): Add new distributions.
+	* oct-rand.h: Provide decls for new functions.
+	(octave_rand::matrix, octave_rand::scalar, octave_rand::
+	(octave_rand::scalar, octave_rand::matrix, octave_rand::nd_array,
+	octave_rand::vector): New arg A, for gamma and poisson distributions.
+	* randpoisson.c, randpoisson.h, randgamma.c, randmtzig.c,
+        randmtzig.h: New files.
+
 2006-03-24  John W. Eaton  <jwe@octave.org>
 
 	* dSparse.cc (SparseMatrix::bsolve): Integer work vector is

--- a/liboctave/Makefile.in	Mon Apr 03 19:03:31 2006 +0000
+++ b/liboctave/Makefile.in	Mon Apr 03 19:18:26 2006 +0000
@@ -69,8 +69,9 @@
 	oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h \
 	oct-shlib.h oct-sort.h oct-spparms.h oct-syscalls.h \
 	oct-sparse.h oct-time.h oct-types.h oct-uname.h \
-	pathlen.h pathsearch.h \
-	prog-args.h so-array.h sparse-sort.h statdefs.h str-vec.h \
+	pathlen.h pathsearch.h prog-args.h \
+	randgamma.h randmtzig.h randpoisson.h \
+	so-array.h sparse-sort.h statdefs.h str-vec.h \
 	sparse-util.h sun-utils.h sysdir.h systime.h syswait.h \
 	$(OPTS_INC) \
 	$(MATRIX_INC) \
@@ -128,7 +129,8 @@
 	$(SPARSE_MX_OP_SRC)
 
 LIBOCTAVE_C_SOURCES := f2c-main.c filemode.c getopt.c getopt1.c \
-	lo-cieee.c lo-cutils.c mkdir.c oct-getopt.c rename.c \
+	lo-cieee.c lo-cutils.c mkdir.c oct-getopt.c \
+	randgamma.c randmtzig.c randpoisson.c rename.c \
 	rmdir.c strftime.c strptime.c tempname.c tempnam.c
 
 LIBOCTAVE_SOURCES := $(LIBOCTAVE_CXX_SOURCES) $(LIBOCTAVE_C_SOURCES)

--- a/liboctave/oct-rand.cc	Mon Apr 03 19:03:31 2006 +0000
+++ b/liboctave/oct-rand.cc	Mon Apr 03 19:18:26 2006 +0000
@@ -24,22 +24,34 @@
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
+#include <vector>
 
 #include "f77-fcn.h"
+#include "lo-ieee.h"
 #include "lo-error.h"
+#include "lo-mappers.h"
 #include "oct-rand.h"
 #include "oct-time.h"
+#include "data-conv.h"
+#include "randmtzig.h"
+#include "randpoisson.h"
+#include "randgamma.h"
 
 // Possible distributions of random numbers.  This was handled with an
 // enum, but unwind_protecting that doesn't work so well.
 #define uniform_dist 1
 #define normal_dist 2
+#define expon_dist 3
+#define poisson_dist 4
+#define gamma_dist 5
 
 // Current distribution of random numbers.
 static int current_distribution = uniform_dist;
 
-// Has the seed been set yet?
-static bool initialized = false;
+// Has the seed/state been set yet?
+static bool old_initialized = false;
+static bool new_initialized = false;
+static bool use_old_generators = false;
 
 extern "C"
 {
@@ -50,6 +62,15 @@
   F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&);
 
   F77_RET_T
+  F77_FUNC (dgenexp, DGENEXP) (const double&, double&);
+
+  F77_RET_T
+  F77_FUNC (dignpoi, DIGNPOI) (const double&, double&);
+
+  F77_RET_T
+  F77_FUNC (dgengam, DGENGAM) (const double&, const double&, double&);
+
+  F77_RET_T
   F77_FUNC (setall, SETALL) (const octave_idx_type&, const octave_idx_type&);
 
   F77_RET_T
@@ -83,7 +104,7 @@
 // work ok to give fairly different seeds each time Octave starts.
 
 static void
-do_initialization (void)
+do_old_initialization (void)
 {
   octave_localtime tm;
  
@@ -99,20 +120,32 @@
 
   F77_FUNC (setall, SETALL) (s0, s1);
 
-  initialized = true;
+  old_initialized = true;
 }
 
 static inline void
 maybe_initialize (void)
 {
-  if (! initialized)
-    do_initialization ();
+  if (use_old_generators)
+    {
+      if (! old_initialized)
+	do_old_initialization ();
+    }
+  else
+    {
+      if (! new_initialized)
+	{
+	  oct_init_by_entropy ();
+	  new_initialized = true;
+	}
+    }
 }
 
 double
 octave_rand::seed (void)
 {
-  maybe_initialize ();
+  if (! old_initialized)
+    do_old_initialization ();
 
   union d2i { double d; octave_idx_type i[2]; };
   union d2i u;
@@ -123,6 +156,7 @@
 void
 octave_rand::seed (double s)
 {
+  use_old_generators = true;
   maybe_initialize ();
 
   union d2i { double d; octave_idx_type i[2]; };
@@ -133,6 +167,41 @@
   F77_FUNC (setsd, SETSD) (i0, i1);
 }
 
+ColumnVector
+octave_rand::state (void)
+{
+  ColumnVector s (MT_N + 1);
+  if (! new_initialized)
+    {
+      oct_init_by_entropy ();
+      new_initialized = true;
+    }
+
+  OCTAVE_LOCAL_BUFFER (unsigned FOUR_BYTE_INT, 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::state (const ColumnVector &s)
+{
+  use_old_generators = false;
+  maybe_initialize ();
+
+  octave_idx_type len = s.length();
+  octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
+  OCTAVE_LOCAL_BUFFER (unsigned FOUR_BYTE_INT, tmp, MT_N + 1);
+  for (octave_idx_type i = 0; i < n; i++)
+    tmp[i] = static_cast<unsigned FOUR_BYTE_INT> (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);
+}
+
 std::string
 octave_rand::distribution (void)
 {
@@ -142,6 +211,12 @@
     return "uniform";
   else if (current_distribution == normal_dist)
     return "normal";
+  else if (current_distribution == expon_dist)
+    return "exponential";
+  else if (current_distribution == poisson_dist)
+    return "poisson";
+  else if (current_distribution == gamma_dist)
+    return "gamma";
   else
     {
       abort ();
@@ -156,6 +231,12 @@
     octave_rand::uniform_distribution ();
   else if (d == "normal")
     octave_rand::normal_distribution ();
+  else if (d == "exponential")
+    octave_rand::exponential_distribution ();
+  else if (d == "poisson")
+    octave_rand::poisson_distribution ();
+  else if (d == "gamma")
+    octave_rand::gamma_distribution ();
   else
     (*current_liboctave_error_handler) ("rand: invalid distribution");
 }
@@ -180,114 +261,105 @@
   F77_FUNC (setcgn, SETCGN) (normal_dist);
 }
 
+void
+octave_rand::exponential_distribution (void)
+{
+  maybe_initialize ();
+
+  current_distribution = expon_dist;
+
+  F77_FUNC (setcgn, SETCGN) (expon_dist);
+}
+
+void
+octave_rand::poisson_distribution (void)
+{
+  maybe_initialize ();
+
+  current_distribution = poisson_dist;
+
+  F77_FUNC (setcgn, SETCGN) (poisson_dist);
+}
+
+void
+octave_rand::gamma_distribution (void)
+{
+  maybe_initialize ();
+
+  current_distribution = gamma_dist;
+
+  F77_FUNC (setcgn, SETCGN) (gamma_dist);
+}
+
+
 double
-octave_rand::scalar (void)
+octave_rand::scalar (double a)
 {
   maybe_initialize ();
 
   double retval = 0.0;
 
-  switch (current_distribution)
+  if (use_old_generators)
     {
-    case uniform_dist:
-      F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
-      break;
+      switch (current_distribution)
+	{
+	case uniform_dist:
+	  F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
+	  break;
 
-    case normal_dist:
-      F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
-      break;
+	case normal_dist:
+	  F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
+	  break;
 
-    default:
-      abort ();
-      break;
-    }
-
-  return retval;
-}
+	case expon_dist:
+	  F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
+	  break;
 
-#define MAKE_RAND_MAT(mat, nr, nc, f, F) \
-  do \
-    { \
-      double val; \
-      for (volatile octave_idx_type j = 0; j < nc; j++) \
-	for (volatile octave_idx_type i = 0; i < nr; i++) \
-	  { \
-	    OCTAVE_QUIT; \
-	    F77_FUNC (f, F) (0.0, 1.0, val); \
-	    mat(i,j) = val; \
-	  } \
-    } \
-  while (0)
-
-Matrix
-octave_rand::matrix (octave_idx_type n, octave_idx_type m)
-{
-  maybe_initialize ();
-
-  Matrix retval;
+	case poisson_dist:
+	  if (a < 0.0 || xisnan(a) || xisinf(a))
+	    retval = octave_NaN;
+	  else
+	    {
+	      // workaround bug in ignpoi, by calling with different Mu
+	      F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
+	      F77_FUNC (dignpoi, DIGNPOI) (a, retval);
+	    }
+	  break;
 
-  if (n >= 0 && m >= 0)
-    {
-      retval.resize (n, m);
+	case gamma_dist:
+	  if (a <= 0.0 || xisnan(a) || xisinf(a))
+	    retval = octave_NaN;
+	  else
+	    F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
+	  break;
 
-      if (n > 0 && m > 0)
-	{
-	  switch (current_distribution)
-	    {
-	    case uniform_dist:
-	      MAKE_RAND_MAT (retval, n, m, dgenunf, DGENUNF);
-	      break;
-
-	    case normal_dist:
-	      MAKE_RAND_MAT (retval, n, m, dgennor, DGENNOR);
-	      break;
-
-	    default:
-	      abort ();
-	      break;
-	    }
+	default:
+	  abort ();
+	  break;
 	}
     }
   else
-    (*current_liboctave_error_handler) ("rand: invalid negative argument");
-
-  return retval;
-}
-
-#define MAKE_RAND_ND_ARRAY(mat, len, f, F) \
-  do \
-    { \
-      double val; \
-      for (volatile octave_idx_type i = 0; i < len; i++) \
-	{ \
-	  OCTAVE_QUIT; \
-	  F77_FUNC (f, F) (0.0, 1.0, val); \
-	  mat(i) = val; \
-	} \
-    } \
-  while (0)
-
-NDArray
-octave_rand::nd_array (const dim_vector& dims)
-{
-  maybe_initialize ();
-
-  NDArray retval;
-
-  if (! dims.all_zero ())
     {
-      retval.resize (dims);
-
-      octave_idx_type len = retval.length ();
-
       switch (current_distribution)
 	{
 	case uniform_dist:
-	  MAKE_RAND_ND_ARRAY (retval, len, dgenunf, DGENUNF);
+	  retval = oct_randu();
 	  break;
 
 	case normal_dist:
-	  MAKE_RAND_ND_ARRAY (retval, len, dgennor, DGENNOR);
+	  retval = oct_randn();
+	  break;
+
+	case expon_dist:
+	  retval = oct_rande();
+	  break;
+
+	case poisson_dist:
+	  retval = oct_randp(a);
+	  break;
+
+	case gamma_dist:
+	  retval = oct_randg(a);
 	  break;
 
 	default:
@@ -299,21 +371,142 @@
   return retval;
 }
 
-#define MAKE_RAND_ARRAY(array, n, f, F) \
+#define MAKE_RAND(len) \
   do \
     { \
       double val; \
-      for (volatile octave_idx_type i = 0; i < n; i++) \
+      for (volatile octave_idx_type i = 0; i < len; i++) \
 	{ \
 	  OCTAVE_QUIT; \
-	  F77_FUNC (f, F) (0.0, 1.0, val); \
-	  array(i) = val; \
+	  RAND_FUNC (val); \
+	  v[i] = val; \
 	} \
     } \
   while (0)
 
+static void
+fill_rand (octave_idx_type len, double *v, double a)
+{
+  maybe_initialize ();
+
+  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_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_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_rande (len, v);
+      break;
+
+    case poisson_dist:
+      if (use_old_generators)
+	{
+	  if (a < 0.0 || xisnan(a) || xisinf(a))
+#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) (a + 1, tmp);
+#define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (a, x)
+		MAKE_RAND (len);
+#undef RAND_FUNC
+	    }
+	}
+      else
+	oct_fill_randp (a, len, v);
+      break;
+
+    case gamma_dist:
+      if (use_old_generators)
+	{
+	  if (a <= 0.0 || xisnan(a) || xisinf(a))
+#define RAND_FUNC(x) x = octave_NaN;
+	    MAKE_RAND (len);
+#undef RAND_FUNC
+	  else
+#define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, a, x)
+	    MAKE_RAND (len);
+#undef RAND_FUNC
+	}
+      else
+	oct_fill_randg (a, len, v);
+      break;
+
+    default:
+      abort ();
+      break;
+    }
+
+  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)
+octave_rand::vector (octave_idx_type n, double a)
 {
   maybe_initialize ();
 
@@ -323,20 +516,7 @@
     {
       retval.resize (n);
 
-      switch (current_distribution)
-	{
-	case uniform_dist:
-	  MAKE_RAND_ARRAY (retval, n, dgenunf, DGENUNF);
-	  break;
-
-	case normal_dist:
-	  MAKE_RAND_ARRAY (retval, n, dgennor, DGENNOR);
-	  break;
-
-	default:
-	  abort ();
-	  break;
-	}
+      fill_rand (retval.capacity(), retval.fortran_vec(), a);
     }
   else if (n < 0)
     (*current_liboctave_error_handler) ("rand: invalid negative argument");

--- a/liboctave/oct-rand.h	Mon Apr 03 19:03:31 2006 +0000
+++ b/liboctave/oct-rand.h	Mon Apr 03 19:18:26 2006 +0000
@@ -26,6 +26,7 @@
 
 #include <string>
 
+#include "dColVector.h"
 #include "dMatrix.h"
 #include "dNDArray.h"
 
@@ -38,6 +39,12 @@
   // Set the seed.
   static void seed (double s);
 
+  // Return the current state.
+  static ColumnVector state (void);
+
+  // Set the current state/
+  static void state (const ColumnVector &s);
+  
   // Return the current distribution.
   static std::string distribution (void);
 
@@ -49,19 +56,25 @@
 
   static void normal_distribution (void);
 
+  static void exponential_distribution (void);
+
+  static void poisson_distribution (void);
+
+  static void gamma_distribution (void);
+
   // Return the next number from the sequence.
-  static double scalar (void);
+  static double scalar (double a = 1.);
 
   // Return a matrix of numbers from the sequence, filled in column
   // major order.
-  static Matrix matrix (octave_idx_type r, octave_idx_type c);
+  static Matrix 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.
-  static NDArray nd_array (const dim_vector& dims);
+  static NDArray nd_array (const dim_vector& dims, double a = 1.);
 
   // Return an array of numbers from the sequence.
-  static Array<double> vector (octave_idx_type n);
+  static Array<double> vector (octave_idx_type n, double a = 1.);
 };
 
 #endif

--- a/src/ChangeLog	Mon Apr 03 19:03:31 2006 +0000
+++ b/src/ChangeLog	Mon Apr 03 19:18:26 2006 +0000
@@ -1,5 +1,12 @@
 2006-04-03  David Bateman  <dbateman@free.fr>
 
+	* DLD-FUNCTIONS/rand.cc (do_rand): Additional argument for
+        gamma and poisson distributions.  Change "state" and "seed"
+        arguments so that they choose between generators.
+	Add, poisson, gamma and exponential generators.
+	(Frand, Frandn): Update docs for new generators, add tests.
+	(Frande, Frandp, Frandg): New generators, with test code.
+
 	* DLD-FUNCTIONS/daspk.cc (Fdaspk): Allow functions to be passed
 	using function handles, inline functions, and cell arrays of
 	strings, inline and function handles.

--- a/src/DLD-FUNCTIONS/rand.cc	Mon Apr 03 19:03:31 2006 +0000
+++ b/src/DLD-FUNCTIONS/rand.cc	Mon Apr 03 19:18:26 2006 +0000
@@ -42,28 +42,56 @@
 #include "utils.h"
 
 static octave_value
-do_rand (const octave_value_list& args, int nargin, const char *fcn)
+do_rand (const octave_value_list& args, int nargin, const char *fcn,
+	 bool additional_arg = false)
 {
   octave_value retval;
+  NDArray a;
+  int idx = 0;
+  dim_vector dims;
 
-  dim_vector dims;
+  if (additional_arg)
+    {
+      if (nargin == 0)
+	{
+	  error ("%s: expecting at least one argument", fcn);
+	  goto done;
+	}
+      else if (args(0).is_string())
+	additional_arg = false;
+      else
+	{
+	  a = args(0).array_value ();
+	  if (error_state)
+	    {
+	      error ("%s: expecting scalar or matrix arguments", fcn);
+	      goto done;
+	    }
+	  idx++;
+	  nargin--;
+	}
+    }
 
   switch (nargin)
     {
     case 0:
       {
-	dims.resize (2);
+	if (additional_arg)
+	  dims = a.dims ();
+	else
+	  {
+	    dims.resize (2);
 
-	dims(0) = 1;
-	dims(1) = 1;
-
+	    dims(0) = 1;
+	    dims(1) = 1;
+	  }
 	goto gen_matrix;
       }
       break;
 
     case 1:
       {
-	octave_value tmp = args(0);
+	octave_value tmp = args(idx);
 
 	if (tmp.is_string ())
 	  {
@@ -73,10 +101,14 @@
 	      {
 		retval = octave_rand::distribution ();
 	      }
-	    else if (s_arg == "seed" || s_arg == "state")
+	    else if (s_arg == "seed")
 	      {
 		retval = octave_rand::seed ();
 	      }
+	    else if (s_arg == "state" || s_arg == "twister")
+	      {
+		retval = octave_rand::state ();
+	      }
 	    else if (s_arg == "uniform")
 	      {
 		octave_rand::uniform_distribution ();
@@ -85,6 +117,18 @@
 	      {
 		octave_rand::normal_distribution ();
 	      }
+	    else if (s_arg == "exponential")
+	      {
+		octave_rand::exponential_distribution ();
+	      }
+	    else if (s_arg == "poisson")
+	      {
+		octave_rand::poisson_distribution ();
+	      }
+	    else if (s_arg == "gamma")
+	      {
+		octave_rand::gamma_distribution ();
+	      }
 	    else
 	      error ("%s: unrecognized string argument", fcn);
 	  }
@@ -176,19 +220,27 @@
 
     default:
       {
-	octave_value tmp = args(0);
+	octave_value tmp = args(idx);
 
 	if (nargin == 2 && tmp.is_string ())
 	  {
 	    std::string ts = tmp.string_value ();
 
-	    if (ts == "seed" || ts == "state")
+	    if (ts == "seed")
 	      {
-		double d = args(1).double_value ();
+		double d = args(idx+1).double_value ();
 
 		if (! error_state)
 		  octave_rand::seed (d);
 	      }
+	    else if (ts == "state" || ts == "twister")
+	      {
+		ColumnVector s = 
+		  ColumnVector (args(idx+1).vector_value(false, true));
+
+		if (! error_state)
+		  octave_rand::state (s);
+	      }
 	    else
 	      error ("%s: unrecognized string argument", fcn);
 	  }
@@ -198,7 +250,7 @@
 
 	    for (int i = 0; i < nargin; i++)
 	      {
-		dims(i) = (octave_idx_type)args(i).int_value ();
+		dims(i) = (octave_idx_type)args(idx+i).int_value ();
 
 		if (error_state)
 		  {
@@ -221,33 +273,97 @@
 
   dims.chop_trailing_singletons ();
 
-  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, ,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {} rand (@var{x})\n\
 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\
-@deftypefnx {Loadable Function} {} rand (@code{\"seed\"}, @var{x})\n\
+@deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\
+@deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\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}.  In\n\
-addition, you can set the seed for the random number generator using the\n\
+for @code{eye}.\n\
+\n\
+You can query the state of the random number generator using the\n\
 form\n\
 \n\
 @example\n\
-rand (\"seed\", @var{x})\n\
+v = rand (\"state\")\n\
+@end example\n\
+\n\
+This returns a column vector @var{v} of length 625. Later, you can\n\
+restore the random number generator to the state @var{v}\n\
+using the form\n\
+\n\
+@example\n\
+rand (\"state\", v)\n\
 @end example\n\
 \n\
 @noindent\n\
-where @var{x} is a scalar value.  If called as\n\
+You may also initialize the state vector from an arbitrary vector of\n\
+length <= 625 for @var{v}.  This new state will be a hash based on the\n\
+the value of @var{v}, not @var{v} itself.\n\
+\n\
+By default, the generator is initialized from @code{/dev/urandom} if it is\n\
+available, otherwise from cpu time, wall clock time and the current\n\
+fraction of a second.\n\
+\n\
+@code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\
+(See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\
+equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\
+Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\
+@url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\
+Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\
+values together, otherwise the generator state can be learned after\n\
+reading 624 consecutive values.\n\
+\n\
+@code{rand} includes a second random number generator, that was the\n\
+previous generator used in octave. The new generator is used by default\n\
+as it is significantly faster than the old generator, and produces\n\
+random numebrs with a significantly longer cycle time. However, in\n\
+some circumstances it might be desireable to obtain the same random\n\
+sequences as used by the old generators. To do this the keyword\n\
+\"seed\" is used to specify that the old generators should be use,\n\
+as in\n\
 \n\
 @example\n\
-rand (\"seed\")\n\
+rand (\"seed\", val)\n\
 @end example\n\
 \n\
-@noindent\n\
-@code{rand} returns the current value of the seed.\n\
+which sets the seed of the generator to @var{val}. The seed of the\n\
+generator can be queried with\n\
+\n\
+@example\n\
+s = rand (\"seed\")\n\
+@end example\n\
+\n\
+However, it should be noted that querying the seed will not cause\n\
+@code{rand} to use the old generators, only setting the seed will.\n\
+To cause @code{rand} to once again use the new generators, the\n\
+keyword \"state\" should be used to reset the state of the @code{rand}.\n\
+@seealso{randn,rande,randg,randp}\n\
 @end deftypefn")
 {
   octave_value retval;
@@ -259,6 +375,64 @@
   return retval;
 }
 
+/*
+%!test # 'state' can be a scalar
+%! rand('state',12); x = rand(1,4);
+%! rand('state',12); y = rand(1,4);
+%! assert(x,y);
+%!test # 'state' can be a vector
+%! rand('state',[12,13]); x=rand(1,4);
+%! rand('state',[12;13]); y=rand(1,4);
+%! assert(x,y);
+%!test # querying 'state' doesn't disturb sequence
+%! rand('state',12); rand(1,2); x=rand(1,2);
+%! rand('state',12); rand(1,2);
+%! s=rand('state'); y=rand(1,2);
+%! assert(x,y);
+%! rand('state',s); z=rand(1,2);
+%! assert(x,z);
+%!test # 'seed' must be a scalar
+%! rand('seed',12); x = rand(1,4);
+%! rand('seed',12); y = rand(1,4);
+%! assert(x,y);
+%!error(rand('seed',[12,13]))
+%!test # querying 'seed' returns a value which can be used later
+%! s=rand('seed'); x=rand(1,2);
+%! rand('seed',s); y=rand(1,2);
+%! assert(x,y);
+%!test # querying 'seed' doesn't disturb sequence
+%! rand('seed',12); rand(1,2); x=rand(1,2);
+%! rand('seed',12); rand(1,2);
+%! s=rand('seed'); y=rand(1,2);
+%! assert(x,y);
+%! rand('seed',s); z=rand(1,2);
+%! assert(x,z);
+*/
+
+/*
+%!test
+%! % statistical tests may fail occasionally.
+%! rand("state",12);
+%! x = rand(100000,1);
+%! assert(max(x)<1.); %*** Please report this!!! ***
+%! assert(min(x)>0.); %*** Please report this!!! ***
+%! assert(mean(x),0.5,0.0024);
+%! assert(var(x),1/48,0.0632);
+%! assert(skewness(x),0,0.012); 
+%! assert(kurtosis(x),-6/5,0.0094);
+%!test
+%! % statistical tests may fail occasionally.
+%! rand("seed",12);
+%! x = rand(100000,1);
+%! assert(max(x)<1.); %*** Please report this!!! ***
+%! assert(min(x)>0.); %*** Please report this!!! ***
+%! assert(mean(x),0.5,0.0024);
+%! assert(var(x),1/48,0.0632);
+%! assert(skewness(x),0,0.012); 
+%! assert(kurtosis(x),-6/5,0.0094);
+*/
+
+
 static std::string current_distribution = octave_rand::distribution ();
 
 static void
@@ -271,25 +445,18 @@
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {} randn (@var{x})\n\
 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\
-@deftypefnx {Loadable Function} {} randn (@code{\"seed\"}, @var{x})\n\
-Return a matrix with normally distributed random elements.  The\n\
-arguments are handled the same as the arguments for @code{eye}.  In\n\
-addition, you can set the seed for the random number generator using the\n\
-form\n\
-\n\
-@example\n\
-randn (\"seed\", @var{x})\n\
-@end example\n\
+@deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\
+@deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\
+Return a matrix with normally distributed random elements. The\n\
+arguments are handled the same as the arguments for @code{rand}.\n\
 \n\
-@noindent\n\
-where @var{x} is a scalar value.  If called as\n\
+By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\
+transform from a uniform to a normal distribution. (G. Marsaglia and\n\
+W.W. Tsang, 'Ziggurat method for generating random variables',\n\
+J. Statistical Software, vol 5, 2000,\n\
+@url{http://www.jstatsoft.org/v05/i08/})\n\
 \n\
-@example\n\
-randn (\"seed\")\n\
-@end example\n\
-\n\
-@noindent\n\
-@code{randn} returns the current value of the seed.\n\
+@seealso{rand,rande,randg,randp}\n\
 @end deftypefn")
 {
   octave_value retval;
@@ -318,6 +485,369 @@
 }
 
 /*
+%!test
+%! % statistical tests may fail occasionally.
+%! rand("state",12);
+%! x = randn(100000,1);
+%! assert(mean(x),0,0.01);
+%! assert(var(x),1,0.02);
+%! assert(skewness(x),0,0.02);
+%! assert(kurtosis(x),0,0.04);
+%!test
+%! % statistical tests may fail occasionally.
+%! rand("seed",12);
+%! x = randn(100000,1);
+%! assert(mean(x),0,0.01);
+%! assert(var(x),1,0.02);
+%! assert(skewness(x),0,0.02);
+%! assert(kurtosis(x),0,0.04);
+*/
+
+DEFUN_DLD (rande, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} rande (@var{x})\n\
+@deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\
+@deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\
+@deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\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 a Marsaglia and Tsang Ziggurat technique to\n\
+transform from a uniform to a exponential distribution. (G. Marsaglia and\n\
+W.W. Tsang, 'Ziggurat method for generating random variables',\n\
+J. Statistical Software, vol 5, 2000,\n\
+@url{http://www.jstatsoft.org/v05/i08/})\n\
+@seealso{rand,randn,randg,randp}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  unwind_protect::begin_frame ("rande");
+
+  // This relies on the fact that elements are popped from the unwind
+  // stack in the reverse of the order they are pushed
+  // (i.e. current_distribution will be reset before calling
+  // reset_rand_generator()).
+
+  unwind_protect::add (reset_rand_generator, 0);
+  unwind_protect_str (current_distribution);
+
+  current_distribution = "exponential";
+
+  octave_rand::distribution (current_distribution);
+
+  retval = do_rand (args, nargin, "rande");
+
+  unwind_protect::run_frame ("rande");
+
+  return retval;
+}
+
+/*
+%!test
+%! % statistical tests may fail occasionally
+%! rand("state",12);
+%! x = rande(100000,1);
+%! assert(min(x)>0); % *** Please report this!!! ***
+%! assert(mean(x),1,0.01);
+%! assert(var(x),1,0.03);
+%! assert(skewness(x),2,0.06);
+%! assert(kurtosis(x),6,0.7);
+%!test
+%! % statistical tests may fail occasionally
+%! rand("seed",12);
+%! x = rande(100000,1);
+%! assert(min(x)>0); % *** Please report this!!! ***
+%! assert(mean(x),1,0.01);
+%! assert(var(x),1,0.03);
+%! assert(skewness(x),2,0.06);
+%! assert(kurtosis(x),6,0.7);
+*/
+
+DEFUN_DLD (randg, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\
+@deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\
+@deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\
+@deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\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\
+\n\
+This can be used to generate many distributions:\n\
+\n\
+@table @asis\n\
+@item @code{gamma (a,b)} for @code{a > -1}, @code{b > 0}\n\
+@example\n\
+r = b*randg(a)\n\
+@end example\n\
+@item @code{beta(a,b)} for @code{a > -1}, @code{b > -1}\n\
+@example\n\
+r1 = randg(a,1)\n\
+r = r1 / (r1 + randg(b,1))\n\
+@end example\n\
+@item @code{Erlang(a, n)}\n\
+@example\n\
+r = a*randg(n)\n\
+@end example\n\
+@item @code{chisq(df)} for @code{df > 0}\n\
+@example\n\
+r = 2*randg(df/2)\n\
+@end example\n\
+@item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\
+@example\n\
+r = randn() / sqrt(2*randg(df/2)/df)\n\
+@end example\n\
+@item @code{F(n1,n2)} for @code{0 < n1}, @code{0 < n2}\n\
+@example\n\
+r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite\n\
+r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite\n\
+r = r1 / r2\n\n\
+@end example\n\
+@item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\
+@example\n\
+r = randp((1-p)/p * randg(n))\n\
+@end example\n\
+@item non-central @code{chisq(df,L)}, for @code{df >= 0} and @code{L > 0}\n\
+(use chisq if @code{L = 0})\n\
+@example\n\
+r = randp(L/2)\n\
+r(r > 0) = 2*randg(r(r > 0))\n\
+r(df > 0) += 2*randg(df(df > 0)/2)\n\
+@end example\n\
+@item @code{Dirichlet(a1,...,ak)}\n\
+@example\n\
+r = (randg(a1),...,randg(ak))\n\
+r = r / sum(r)\n\
+@end example\n\
+@end table\n\
+@seealso{rand,randn,rande,randp}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin < 1)
+    error ("randg: insufficient arguments");
+  else
+    {
+      unwind_protect::begin_frame ("randg");
+
+      // This relies on the fact that elements are popped from the unwind
+      // stack in the reverse of the order they are pushed
+      // (i.e. current_distribution will be reset before calling
+      // reset_rand_generator()).
+
+      unwind_protect::add (reset_rand_generator, 0);
+      unwind_protect_str (current_distribution);
+
+      current_distribution = "gamma";
+
+      octave_rand::distribution (current_distribution);
+
+      retval = do_rand (args, nargin, "randg", true);
+
+      unwind_protect::run_frame ("randg");
+    }
+
+  return retval;
+}
+
+/*
+%!test
+%! rand("state",12)
+%!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report
+%!test
+%! % statistical tests may fail occasionally.
+%! a=0.1; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.01);
+%! assert(var(x),     a,         0.01);
+%! assert(skewness(x),2/sqrt(a), 1.);
+%! assert(kurtosis(x),6/a,       50.);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=0.95; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.01);
+%! assert(var(x),     a,         0.04);
+%! assert(skewness(x),2/sqrt(a), 0.2);
+%! assert(kurtosis(x),6/a,       2.);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=1; x = randg(a,100000,1);
+%! assert(mean(x),a,             0.01);
+%! assert(var(x),a,              0.04);
+%! assert(skewness(x),2/sqrt(a), 0.2);
+%! assert(kurtosis(x),6/a,       2.);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=10; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.1);
+%! assert(var(x),     a,         0.5);
+%! assert(skewness(x),2/sqrt(a), 0.1);
+%! assert(kurtosis(x),6/a,       0.5);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=100; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.2);
+%! assert(var(x),     a,         2.);
+%! assert(skewness(x),2/sqrt(a), 0.05);
+%! assert(kurtosis(x),6/a,       0.2);
+%!test
+%! rand("seed",12)
+%!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report
+%!test
+%! % statistical tests may fail occasionally.
+%! a=0.1; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.01);
+%! assert(var(x),     a,         0.01);
+%! assert(skewness(x),2/sqrt(a), 1.);
+%! assert(kurtosis(x),6/a,       50.);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=0.95; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.01);
+%! assert(var(x),     a,         0.04);
+%! assert(skewness(x),2/sqrt(a), 0.2);
+%! assert(kurtosis(x),6/a,       2.);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=1; x = randg(a,100000,1);
+%! assert(mean(x),a,             0.01);
+%! assert(var(x),a,              0.04);
+%! assert(skewness(x),2/sqrt(a), 0.2);
+%! assert(kurtosis(x),6/a,       2.);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=10; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.1);
+%! assert(var(x),     a,         0.5);
+%! assert(skewness(x),2/sqrt(a), 0.1);
+%! assert(kurtosis(x),6/a,       0.5);
+%!test
+%! % statistical tests may fail occasionally.
+%! a=100; x = randg(a,100000,1);
+%! assert(mean(x),    a,         0.2);
+%! assert(var(x),     a,         2.);
+%! assert(skewness(x),2/sqrt(a), 0.05);
+%! assert(kurtosis(x),6/a,       0.2);
+*/
+
+
+DEFUN_DLD (randp, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\
+@deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\
+@deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\
+@deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\
+Return a matrix with Poisson distributed random elements. The arguments\n\
+are handled the same as the arguments for @code{rand}, except for the\n\
+argument @var{l}.\n\
+\n\
+Five different algorithms are used depending on the range of @var{l}\n\
+and whether or not @var{l} is a scalar or a matrix.\n\
+\n\
+@table @asis\n\
+@item For scalar @var{l} <= 12, use direct method.\n\
+Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\
+@item For scalar @var{l} > 12, use rejection method.[1]\n\
+Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\
+@item For matrix @var{l} <= 10, use inversion method.[2]\n\
+Stadlober E., et al., WinRand source code, available via FTP.\n\
+@item For matrix @var{l} > 10, use patchwork rejection method.\n\
+Stadlober E., et al., WinRand source code, available via FTP, or\n\
+H. Zechner, 'Efficient sampling from continuous and discrete\n\
+unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\
+University Graz, Austria, 1994.\n\
+@item For @var{l} > 1e8, use normal approximation.\n\
+L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\
+D 50 p1284, 1994\n\
+@end table\n\
+@seealso{rand,randn,rande,randg}\n\
+@end deftypefn")
+{
+  octave_value retval;
+
+  int nargin = args.length ();
+
+  if (nargin < 1)
+    error ("randp: insufficient arguments");
+  else
+    {
+      unwind_protect::begin_frame ("randp");
+
+      // This relies on the fact that elements are popped from the unwind
+      // stack in the reverse of the order they are pushed
+      // (i.e. current_distribution will be reset before calling
+      // reset_rand_generator()).
+
+      unwind_protect::add (reset_rand_generator, 0);
+      unwind_protect_str (current_distribution);
+
+      current_distribution = "poisson";
+
+      octave_rand::distribution (current_distribution);
+
+      retval = do_rand (args, nargin, "randp", true);
+
+      unwind_protect::run_frame ("randp");
+    }
+
+  return retval;
+}
+
+/*
+%!test
+%! rand("state",12)
+%!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report
+%!test
+%! % statistical tests may fail occasionally.
+%! for a=[5 15]
+%!   x = randp(a,100000,1);
+%!   assert(min(x)>=0); % *** Please report this!!! ***
+%!   assert(mean(x),a,0.03);
+%!   assert(var(x),a,0.2);
+%!   assert(skewness(x),1/sqrt(a),0.03);
+%!   assert(kurtosis(x),1/a,0.08);
+%! end
+%!test
+%! % statistical tests may fail occasionally.
+%! for a=[5 15]
+%!   x = randp(a*ones(100000,1),100000,1);
+%!   assert(min(x)>=0); % *** Please report this!!! ***
+%!   assert(mean(x),a,0.03);
+%!   assert(var(x),a,0.2);
+%!   assert(skewness(x),1/sqrt(a),0.03);
+%!   assert(kurtosis(x),1/a,0.08);
+%! end
+%!test
+%! rand("seed",12)
+%!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report
+%!test
+%! % statistical tests may fail occasionally.
+%! for a=[5 15]
+%!   x = randp(a,100000,1);
+%!   assert(min(x)>=0); % *** Please report this!!! ***
+%!   assert(mean(x),a,0.03);
+%!   assert(var(x),a,0.2);
+%!   assert(skewness(x),1/sqrt(a),0.03);
+%!   assert(kurtosis(x),1/a,0.08);
+%! end
+%!test
+%! % statistical tests may fail occasionally.
+%! for a=[5 15]
+%!   x = randp(a*ones(100000,1),100000,1);
+%!   assert(min(x)>=0); % *** Please report this!!! ***
+%!   assert(mean(x),a,0.03);
+%!   assert(var(x),a,0.2);
+%!   assert(skewness(x),1/sqrt(a),0.03);
+%!   assert(kurtosis(x),1/a,0.08);
+%! end
+*/
+
+/*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
 ;;; End: ***