diff --git a/src/graph/graph.hh b/src/graph/graph.hh
index 73eacb6bf..fec5e0ac7 100644
--- a/src/graph/graph.hh
+++ b/src/graph/graph.hh
@@ -831,7 +831,7 @@ namespace spot
     {
       //std::cerr << "\nbefore\n";
       //dump_storage(std::cerr);
-      std::sort(transitions_.begin() + 1, transitions_.end(), p);
+      std::stable_sort(transitions_.begin() + 1, transitions_.end(), p);
     }
 
     // Should be called only when it is known that all transitions
diff --git a/src/ltlvisit/simplify.cc b/src/ltlvisit/simplify.cc
index f5e55b79f..16f228c10 100644
--- a/src/ltlvisit/simplify.cc
+++ b/src/ltlvisit/simplify.cc
@@ -533,26 +533,36 @@ namespace spot
 				    const formula* f1,
 				    const formula* f2)
 	{
-	  // Rewrite a<=>b as (a&b)|(!a&!b)
 	  if (equiv)
-	    return
-	      multop::instance(multop::Or,
-			       multop::instance(multop::And,
-						recurse_(f1, false),
-						recurse_(f2, false)),
-			       multop::instance(multop::And,
-						recurse_(f1, true),
-						recurse_(f2, true)));
+	    {
+	      // Rewrite a<=>b as (a&b)|(!a&!b)
+	      auto recurse_f1_true = recurse_(f1, true);
+	      auto recurse_f1_false = recurse_(f1, false);
+	      auto recurse_f2_true =  recurse_(f2, true);
+	      auto recurse_f2_false =  recurse_(f2, false);
+	      auto left =  multop::instance(multop::And,
+					    recurse_f1_false,
+					    recurse_f2_false);
+	      auto right =  multop::instance(multop::And,
+					    recurse_f1_true,
+					    recurse_f2_true);
+	      return multop::instance(multop::Or, left, right);
+	    }
 	  else
-	    // Rewrite a^b as (a&!b)|(!a&b)
-	    return
-	      multop::instance(multop::Or,
-			       multop::instance(multop::And,
-						recurse_(f1, false),
-						recurse_(f2, true)),
-			       multop::instance(multop::And,
-						recurse_(f1, true),
-						recurse_(f2, false)));
+	    {
+	      // Rewrite a^b as (a&!b)|(!a&b)
+	      auto recurse_f1_true = recurse_(f1, true);
+	      auto recurse_f1_false = recurse_(f1, false);
+	      auto recurse_f2_true =  recurse_(f2, true);
+	      auto recurse_f2_false =  recurse_(f2, false);
+	      auto left  =  multop::instance(multop::And,
+					     recurse_f1_false,
+					     recurse_f2_true);
+	      auto right = multop::instance(multop::And,
+					    recurse_f1_true,
+					    recurse_f2_false);
+	      return multop::instance(multop::Or, left, right);
+	    }
 	}
 
 	void
diff --git a/src/misc/random.cc b/src/misc/random.cc
index 6eac054f4..1f963a514 100644
--- a/src/misc/random.cc
+++ b/src/misc/random.cc
@@ -1,5 +1,5 @@
 // -*- coding: utf-8 -*-
-// Copyright (C) 2011, 2012, 2013, 2014 Laboratoire de Recherche et
+// Copyright (C) 2011, 2012, 2013, 2014, 2015 Laboratoire de Recherche et
 // Développement de l'Epita (LRDE).
 // Copyright (C) 2004 Laboratoire d'Informatique de Paris 6 (LIP6),
 // département Systèmes Répartis Coopératifs (SRC), Université Pierre
@@ -22,6 +22,7 @@
 
 #include "config.h"
 #include "random.hh"
+#include <random>
 
 namespace spot
 {
@@ -36,28 +37,96 @@ namespace spot
   double
   drand()
   {
-    return
-      std::generate_canonical<double,
-			      std::numeric_limits<double>::digits>(gen);
+    return  gen() / (1.0 + gen.max());
   }
 
   int
   mrand(int max)
   {
-    std::uniform_int_distribution<> dis(0, max - 1);
-    return dis(gen);
+    return static_cast<int>(max * drand());
   }
 
   int
   rrand(int min, int max)
   {
-    std::uniform_int_distribution<> dis(min, max);
-    return dis(gen);
+    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
-  barand::rand()
+  prand(double p)
   {
-    return (*this)(gen);
+    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 d64ac194a..de121639c 100644
--- a/src/misc/random.hh
+++ b/src/misc/random.hh
@@ -1,5 +1,5 @@
 // -*- coding: utf-8 -*-
-// Copyright (C) 2013 Laboratoire de Recherche et Développement
+// Copyright (C) 2015 Laboratoire de Recherche et Développement
 // de l'Epita (LRDE).
 // Copyright (C) 2004  Laboratoire d'Informatique de Paris 6 (LIP6),
 // département Systèmes Répartis Coopératifs (SRC), Université Pierre
@@ -24,7 +24,8 @@
 # define SPOT_MISC_RANDOM_HH
 
 # include "common.hh"
-# include <random>
+# include <cmath>
+# include <vector>
 
 namespace spot
 {
@@ -55,18 +56,86 @@ namespace spot
   /// \see mrand, rrand, srand
   SPOT_API 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.
+  SPOT_API 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.
+  SPOT_API double bmrand();
+
   /// \brief Compute pseudo-random integer value between 0
   /// and \a n included, following a binomial distribution
   /// for probability \a p.
-  class SPOT_API barand : protected std::binomial_distribution<>
+  ///
+  /// \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) : binomial_distribution(n, p)
+    barand(int n, double p)
+      : n_(n), m_(n * p), s_(sqrt(n * p * (1 - p)))
     {
     }
 
-    int rand();
+    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>
+  SPOT_API int prand(double p);
+
+  /// \brief Shuffle the container using mrand function above.
+  /// This allows to get rid off shuffle or random_shuffle that use
+  /// uniform_distribution and RandomIterator that are not portables.
+  template<class iterator_type>
+  SPOT_API void mrandom_shuffle(iterator_type&& first, iterator_type&& last)
+  {
+    auto d = std::distance(first, last);
+    if (d > 1)
+      {
+	for (--last; first < last; ++first, --d)
+	  {
+	    auto i = mrand(d);
+	    std::swap(*first, *(first + i));
+	  }
+      }
+  }
+  /// @}
 }
 
 #endif // SPOT_MISC_RANDOM_HH
diff --git a/src/tgbaalgos/randomgraph.cc b/src/tgbaalgos/randomgraph.cc
index 52fafff6a..009cec846 100644
--- a/src/tgbaalgos/randomgraph.cc
+++ b/src/tgbaalgos/randomgraph.cc
@@ -157,7 +157,7 @@ namespace spot
     // We want to connect each node to a number of successors between
     // 1 and n.  If the probability to connect to each successor is d,
     // the number of connected successors follows a binomial distribution.
-    barand bin(n - 1, d);
+    barand<nrand> bin(n - 1, d);
 
     while (!nodes_to_process.empty())
       {
diff --git a/src/tgbaalgos/randomize.cc b/src/tgbaalgos/randomize.cc
index 675fdffd0..dad931904 100644
--- a/src/tgbaalgos/randomize.cc
+++ b/src/tgbaalgos/randomize.cc
@@ -19,6 +19,7 @@
 
 #include <algorithm>
 #include <numeric>
+#include <random>
 #include "randomize.hh"
 #include "misc/random.hh"
 
@@ -36,7 +37,7 @@ namespace spot
 	unsigned n = g.num_states();
 	std::vector<unsigned> nums(n);
 	std::iota(nums.begin(), nums.end(), 0);
-	std::random_shuffle(nums.begin(), nums.end(), spot::mrand);
+	mrandom_shuffle(nums.begin(), nums.end());
 	g.rename_states_(nums);
 	aut->set_init_state(nums[aut->get_init_state_number()]);
 
@@ -54,7 +55,7 @@ namespace spot
       {
 	g.remove_dead_transitions_();
 	auto& v = g.transition_vector();
-	std::random_shuffle(v.begin() + 1, v.end(), spot::mrand);
+	mrandom_shuffle(v.begin() + 1, v.end());
       }
 
     typedef tgba_digraph::graph_t::trans_storage_t tr_t;
diff --git a/src/tgbatest/randaut.test b/src/tgbatest/randaut.test
index 20df91e44..7247bf12d 100755
--- a/src/tgbatest/randaut.test
+++ b/src/tgbatest/randaut.test
@@ -58,11 +58,11 @@ test `expr $a + $b` = 100
 # not hesitate to adjust the expected values below.
 $randaut -n 5 --name='%F-%L-%s-%c-%e' -H a | grep '^name' >out
 cat >expected<<EOF
-name: "0-0-10-1-29"
-name: "0-1-10-1-30"
-name: "0-2-10-1-29"
-name: "0-3-10-3-24"
-name: "0-4-10-1-29"
+name: "0-0-10-1-30"
+name: "0-1-10-4-27"
+name: "0-2-10-6-20"
+name: "0-3-10-1-25"
+name: "0-4-10-2-20"
 EOF
 diff out expected
 
diff --git a/src/tgbatest/randomize.test b/src/tgbatest/randomize.test
index 201d07520..84aa716fd 100755
--- a/src/tgbatest/randomize.test
+++ b/src/tgbatest/randomize.test
@@ -86,32 +86,7 @@ grep "unknown argument for --randomize: 'f'" stderr
 cat >input <<EOF
 HOA: v1
 States: 5
-Start: 0
-AP: 4 "a" "b" "c" "d"
-acc-name: all
-Acceptance: 0 t
-properties: trans-labels explicit-labels state-acc
---BODY--
-State: 0 "s0"
-[0] 1
-[1] 2
-[2] 3
-[3] 4
-State: 1 "s1"
-[0] 1
-State: 2 "s2"
-[1] 2
-State: 3 "s3"
-[2] 3
-State: 4 "s4"
-[3] 4
---END--
-EOF
-$autfilt --randomize --seed=1 input -H > output
-cat >expected <<EOF
-HOA: v1
-States: 5
-Start: 1
+Start: 4
 AP: 4 "a" "b" "c" "d"
 acc-name: all
 Acceptance: 0 t
@@ -119,17 +94,42 @@ properties: trans-labels explicit-labels state-acc
 --BODY--
 State: 0 "s1"
 [0] 0
-State: 1 "s0"
-[1] 3
+State: 1 "s3"
+[2] 1
+State: 2 "s2"
+[1] 2
+State: 3 "s4"
+[3] 3
+State: 4 "s0"
+[3] 3
+[2] 1
 [0] 0
-[2] 2
-[3] 4
-State: 2 "s3"
-[2] 2
-State: 3 "s2"
-[1] 3
-State: 4 "s4"
-[3] 4
+[1] 2
+--END--
+EOF
+$autfilt --randomize --seed=1 input -H > output
+cat >expected <<EOF
+HOA: v1
+States: 5
+Start: 3
+AP: 4 "a" "b" "c" "d"
+acc-name: all
+Acceptance: 0 t
+properties: trans-labels explicit-labels state-acc
+--BODY--
+State: 0 "s4"
+[3] 0
+State: 1 "s2"
+[1] 1
+State: 2 "s1"
+[0] 2
+State: 3 "s0"
+[1] 1
+[3] 0
+[2] 4
+[0] 2
+State: 4 "s3"
+[2] 4
 --END--
 EOF
 
diff --git a/src/tgbatest/randpsl.test b/src/tgbatest/randpsl.test
index c8e183b60..916a684e9 100755
--- a/src/tgbatest/randpsl.test
+++ b/src/tgbatest/randpsl.test
@@ -25,7 +25,7 @@ set -e
 
 # Generate 50 random unique PSL formula that do not simplify to LTL
 # formulae, and that have a size of at lease 12.
-../../bin/randltl -n -1 --tree-size 30 --seed 0 --psl a b c |
+../../bin/randltl -n -1 --tree-size 30 --seed 12 --psl a b c |
 ../../bin/ltlfilt -r --size-min 12 --unique |
 ../../bin/ltlfilt -v --ltl -n 50 | tee formulas |
 ../../bin/ltlcross '../ltl2tgba -R3 -t %f >%T' '../ltl2tgba -x -R3 -t %f >%T' \
diff --git a/src/tgbatest/readsave.test b/src/tgbatest/readsave.test
index 98609e465..84b20350b 100755
--- a/src/tgbatest/readsave.test
+++ b/src/tgbatest/readsave.test
@@ -114,16 +114,16 @@ $randltl -n -1 a b |
     $autfilt -F- -F nonexistant --states=3 --edges=..10 --acc-sets=1.. \
 	     --name='%M, %S states' --stats='<%m>, %e, %a' -n 10 > output
 cat >expected <<EOF
-<F(b | GF!a), 3 states>, 6, 1
-<XFb, 3 states>, 4, 1
+<F(b | Ga), 3 states>, 5, 1
+<F(!b & G(!b | G!a)), 3 states>, 5, 1
 <XF!b, 3 states>, 4, 1
-<G!b & XF!a, 3 states>, 4, 1
-<F(b | GFa), 3 states>, 6, 1
-<F(Ga | XG(!a & Fb)), 3 states>, 6, 1
-<FG!b & F(b | GFb), 3 states>, 5, 1
-<Ga | G!a, 3 states>, 4, 1
-<a | G((!a & !b) | (a & b)), 3 states>, 4, 1
-<Fb U G!a, 3 states>, 7, 1
+<Gb | G!b, 3 states>, 4, 1
+<XFb, 3 states>, 4, 1
+<F(b W a), 3 states>, 6, 1
+<(a & !b & (b | (!b M F!a))) | (!a & (b | (!b & (b W Ga)))), 3 states>, 5, 1
+<(a & (a U !b)) | (!a & (!a R b)), 3 states>, 5, 1
+<a | G((a & GFa) | (!a & FG!a)), 3 states>, 4, 1
+<XXG(!a & (Fa W Gb)), 3 states>, 3, 1
 EOF
 diff output expected