* src/misc/random.hh (nrand, bmrand, prand): New functions.

(barand): New class.
* src/misc/random.cc (nrand, bmrand, prand): New functions.
* wrap/python/spot.i: Process src/misc/random.hh.

ChangeLog

@@ -1,5 +1,10 @@
 2004-12-07  Alexandre Duret-Lutz  <adl@src.lip6.fr>
 
+	* src/misc/random.hh (nrand, bmrand, prand): New functions.
+	(barand): New class.
+	* src/misc/random.cc (nrand, bmrand, prand): New functions.
+	* wrap/python/spot.i: Process src/misc/random.hh.
+
 	* src/misc/timer.cc: Do not include cassert, then.
 
 2004-12-07  Denis Poitrenaud  <Denis.Poitrenaud@lip6.fr>

src/misc/random.cc

@@ -57,4 +57,82 @@ namespace spot
     return min + static_cast<int>((max - min + 1) * drand());
   }
 
+  double
+  nrand()
+  {
+    const double r = drand();
+
+    const double lim = 1.e-20;
+    if (r < lim)
+      return -1./lim;
+    if (r > 1.0 - lim)
+      return 1./lim;
+
+    double t;
+    if (r < 0.5)
+      t = sqrt(-2.0 * log(r));
+    else
+      t = sqrt(-2.0 * log(1.0 - r));
+
+    const double p0 = 0.322232431088;
+    const double p1 = 1.0;
+    const double p2 = 0.342242088547;
+    const double p3 = 0.204231210245e-1;
+    const double p4 = 0.453642210148e-4;
+    const double q0 = 0.099348462606;
+    const double q1 = 0.588581570495;
+    const double q2 = 0.531103462366;
+    const double q3 = 0.103537752850;
+    const double q4 = 0.385607006340e-2;
+    const double p = p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)));
+    const double q = q0 + t * (q1 + t * (q2 + t * (q3 + t * q4)));
+
+    if (r < 0.5)
+      return (p / q) - t;
+    else
+      return t - (p / q);
+  }
+
+  double
+  bmrand()
+  {
+    static double next;
+    static bool has_next = false;
+
+    if (has_next)
+      {
+	has_next = false;
+	return next;
+      }
+
+    double x;
+    double y;
+    double r;
+    do
+      {
+	x = 2.0 * drand() - 1.0;
+	y = 2.0 * drand() - 1.0;
+	r = x * x + y * y;
+      }
+    while (r >= 1.0 || r == 0.0);
+    r = sqrt(-2 * log(r) / r);
+    next = y * r;
+    has_next = true;
+    return x * r;
+  }
+
+  int
+  prand(double p)
+  {
+    double s = 0.0;
+    long x = 0;
+
+    while (s < p)
+      {
+	s -= log(1.0 - drand());
+	++x;
+      }
+    return x - 1;
+  }
+
 }

src/misc/random.hh

@@ -19,6 +19,8 @@
 // Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 // 02111-1307, USA.
 
+#include <cmath>
+
 namespace spot
 {
   /// \addtogroup random Random functions
@@ -36,7 +38,7 @@ namespace spot
   /// \see drand, mrand, srand
   int rrand(int min, int max);
 
-  /// \brief Compute a pseudo-random integer value between  0 and
+  /// \brief Compute a pseudo-random integer value between 0 and
   /// \a max-1 included.
   ///
   /// \see drand, rrand, srand
@@ -48,5 +50,67 @@ namespace spot
   /// \see mrand, rrand, srand
   double drand();
 
+  /// \brief Compute a pseudo-random double value
+  /// following a standard normal distribution.  (Odeh & Evans)
+  ///
+  /// This uses a polynomial approximation of the inverse cumulated
+  /// density function from Odeh & Evans, Journal of Applied
+  /// Statistics, 1974, vol 23, pp 96-97.
+  double nrand();
+
+  /// \brief Compute a pseudo-random double value
+  /// following a standard normal distribution.  (Box-Muller)
+  ///
+  /// This uses the polar form of the Box-Muller transform
+  /// to generate random values.
+  double bmrand();
+
+  /// \brief Compute pseudo-random integer value between 0
+  /// and \a n included, following a binomial distribution
+  /// for probability \a p.
+  ///
+  /// \a gen must be a random function computing a pseudo-random
+  /// double value following a standard normal distribution.
+  /// Use nrand() or bmrand().
+  ///
+  /// Usually approximating a binomial distribution using a normal
+  /// distribution and is accurate only if <code>n*p</code> and
+  /// <code>n*(1-p)</code> are greater than 5.
+  template<double (*gen)()>
+  class barand
+  {
+  public:
+    barand(int n, double p)
+      : n_(n), m_(n * p), s_(sqrt(n * p * (1 - p)))
+    {
+    }
+
+    int
+    rand() const
+    {
+      int res;
+
+      for (;;)
+	{
+	  double x = gen() * s_ + m_;
+	  if (x < 0.0)
+	    continue;
+	  res = static_cast<int> (x);
+          if (res <= n_)
+	    break;
+        }
+      return res;
+    }
+  protected:
+    const int n_;
+    const double m_;
+    const double s_;
+  };
+
+  /// \brief Return a pseudo-random positive integer value
+  /// following a Poisson distribution with parameter \a p.
+  ///
+  /// \pre <code>p > 0</code>
+  int prand(double p);
   /// @}
 }

wrap/python/spot.i

@@ -36,6 +36,7 @@
 #include "misc/bddalloc.hh"
 #include "misc/minato.hh"
 #include "misc/modgray.hh"
+#include "misc/random.hh"
 
 #include "ltlast/formula.hh"
 #include "ltlast/refformula.hh"
@@ -95,6 +96,7 @@ using namespace spot;
 %include "misc/version.hh"
 %include "misc/bddalloc.hh"
 %include "misc/minato.hh"
+%include "misc/random.hh"
 
 %feature("director") spot::loopless_modular_mixed_radix_gray_code;
 %include "misc/modgray.hh"