From 541705a36ade2aef4859392c49fb43ed7856ca58 Mon Sep 17 00:00:00 2001 From: Alexandre Duret-Lutz Date: Tue, 7 Dec 2004 18:52:10 +0000 Subject: [PATCH] * 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 | 5 +++ src/misc/random.cc | 78 ++++++++++++++++++++++++++++++++++++++++++++++ src/misc/random.hh | 66 ++++++++++++++++++++++++++++++++++++++- wrap/python/spot.i | 2 ++ 4 files changed, 150 insertions(+), 1 deletion(-) diff --git a/ChangeLog b/ChangeLog index 80c7102ec..e72fd2249 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,10 @@ 2004-12-07 Alexandre Duret-Lutz + * 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 diff --git a/src/misc/random.cc b/src/misc/random.cc index 5682f3167..7e3635a72 100644 --- a/src/misc/random.cc +++ b/src/misc/random.cc @@ -57,4 +57,82 @@ namespace spot return min + static_cast((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; + } + } diff --git a/src/misc/random.hh b/src/misc/random.hh index 684239961..9427e3d1f 100644 --- a/src/misc/random.hh +++ b/src/misc/random.hh @@ -19,6 +19,8 @@ // Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA // 02111-1307, USA. +#include + 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 n*p and + /// n*(1-p) are greater than 5. + template + 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 (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 p > 0 + int prand(double p); /// @} } diff --git a/wrap/python/spot.i b/wrap/python/spot.i index 9418592e4..f7317cd53 100644 --- a/wrap/python/spot.i +++ b/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"