Index: CMakeLists.txt
===================================================================
--- CMakeLists.txt	(revision 1346)
+++ CMakeLists.txt	(revision 1404)
@@ -87,35 +87,34 @@
 IF(LEMON_ENABLE_GLPK) 
   FIND_PACKAGE(GLPK 4.33)
+  IF(GLPK_FOUND)
+    SET(LEMON_HAVE_LP TRUE)
+    SET(LEMON_HAVE_MIP TRUE)
+    SET(LEMON_HAVE_GLPK TRUE)
+  ENDIF(GLPK_FOUND)
 ENDIF(LEMON_ENABLE_GLPK)
 IF(LEMON_ENABLE_ILOG)
   FIND_PACKAGE(ILOG)
+  IF(ILOG_FOUND)
+    SET(LEMON_HAVE_LP TRUE)
+    SET(LEMON_HAVE_MIP TRUE)
+    SET(LEMON_HAVE_CPLEX TRUE)
+  ENDIF(ILOG_FOUND)
 ENDIF(LEMON_ENABLE_ILOG)
 IF(LEMON_ENABLE_COIN)
   FIND_PACKAGE(COIN)
+  IF(COIN_FOUND)
+    SET(LEMON_HAVE_LP TRUE)
+    SET(LEMON_HAVE_MIP TRUE)
+    SET(LEMON_HAVE_CLP TRUE)
+    SET(LEMON_HAVE_CBC TRUE)
+  ENDIF(COIN_FOUND)
 ENDIF(LEMON_ENABLE_COIN)
 IF(LEMON_ENABLE_SOPLEX)
   FIND_PACKAGE(SOPLEX)
+  IF(SOPLEX_FOUND)
+    SET(LEMON_HAVE_LP TRUE)
+    SET(LEMON_HAVE_SOPLEX TRUE)
+  ENDIF(SOPLEX_FOUND)
 ENDIF(LEMON_ENABLE_SOPLEX)
-
-IF(GLPK_FOUND)
-  SET(LEMON_HAVE_LP TRUE)
-  SET(LEMON_HAVE_MIP TRUE)
-  SET(LEMON_HAVE_GLPK TRUE)
-ENDIF(GLPK_FOUND)
-IF(ILOG_FOUND)
-  SET(LEMON_HAVE_LP TRUE)
-  SET(LEMON_HAVE_MIP TRUE)
-  SET(LEMON_HAVE_CPLEX TRUE)
-ENDIF(ILOG_FOUND)
-IF(COIN_FOUND)
-  SET(LEMON_HAVE_LP TRUE)
-  SET(LEMON_HAVE_MIP TRUE)
-  SET(LEMON_HAVE_CLP TRUE)
-  SET(LEMON_HAVE_CBC TRUE)
-ENDIF(COIN_FOUND)
-IF(SOPLEX_FOUND)
-  SET(LEMON_HAVE_LP TRUE)
-  SET(LEMON_HAVE_SOPLEX TRUE)
-ENDIF(SOPLEX_FOUND)
 
 IF(ILOG_FOUND)
Index: cmake/FindCOIN.cmake
===================================================================
--- cmake/FindCOIN.cmake	(revision 1232)
+++ cmake/FindCOIN.cmake	(revision 1398)
@@ -66,4 +66,10 @@
 )
 
+FIND_LIBRARY(COIN_PTHREADS_LIBRARY
+  NAMES pthreads libpthreads
+  HINTS ${COIN_ROOT_DIR}/lib/coin
+  HINTS ${COIN_ROOT_DIR}/lib
+)
+
 INCLUDE(FindPackageHandleStandardArgs)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS(COIN DEFAULT_MSG
@@ -83,5 +89,5 @@
 IF(COIN_FOUND)
   SET(COIN_INCLUDE_DIRS ${COIN_INCLUDE_DIR})
-  SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY};${COIN_ZLIB_LIBRARY};${COIN_BZ2_LIBRARY}")
+  SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARY};${COIN_COIN_UTILS_LIBRARY}")
   IF(COIN_ZLIB_LIBRARY)
     SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_ZLIB_LIBRARY}")
@@ -90,5 +96,8 @@
     SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_BZ2_LIBRARY}")
   ENDIF(COIN_BZ2_LIBRARY)
-  SET(COIN_CBC_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_ZLIB_LIBRARY};${COIN_BZ2_LIBRARY};${COIN_CLP_LIBRARIES}")
+   IF(COIN_PTHREADS_LIBRARY)
+    SET(COIN_CLP_LIBRARIES "${COIN_CLP_LIBRARIES};${COIN_PTHREADS_LIBRARY}")
+  ENDIF(COIN_PTHREADS_LIBRARY)
+  SET(COIN_CBC_LIBRARIES "${COIN_CBC_LIBRARY};${COIN_CBC_SOLVER_LIBRARY};${COIN_CGL_LIBRARY};${COIN_OSI_LIBRARY};${COIN_OSI_CBC_LIBRARY};${COIN_OSI_CLP_LIBRARY};${COIN_CLP_LIBRARIES}")
   SET(COIN_LIBRARIES ${COIN_CBC_LIBRARIES})
 ENDIF(COIN_FOUND)
Index: lemon/adaptors.h
===================================================================
--- lemon/adaptors.h	(revision 1270)
+++ lemon/adaptors.h	(revision 1402)
@@ -3447,8 +3447,8 @@
     /// This map adaptor class adapts two node maps of the original digraph
     /// to get a node map of the split digraph.
-    /// Its value type is inherited from the first node map type (\c IN).
-    /// \tparam IN The type of the node map for the in-nodes.
-    /// \tparam OUT The type of the node map for the out-nodes.
-    template <typename IN, typename OUT>
+    /// Its value type is inherited from the first node map type (\c In).
+    /// \tparam In The type of the node map for the in-nodes.
+    /// \tparam Out The type of the node map for the out-nodes.
+    template <typename In, typename Out>
     class CombinedNodeMap {
     public:
@@ -3457,14 +3457,14 @@
       typedef Node Key;
       /// The value type of the map
-      typedef typename IN::Value Value;
-
-      typedef typename MapTraits<IN>::ReferenceMapTag ReferenceMapTag;
-      typedef typename MapTraits<IN>::ReturnValue ReturnValue;
-      typedef typename MapTraits<IN>::ConstReturnValue ConstReturnValue;
-      typedef typename MapTraits<IN>::ReturnValue Reference;
-      typedef typename MapTraits<IN>::ConstReturnValue ConstReference;
+      typedef typename In::Value Value;
+
+      typedef typename MapTraits<In>::ReferenceMapTag ReferenceMapTag;
+      typedef typename MapTraits<In>::ReturnValue ReturnValue;
+      typedef typename MapTraits<In>::ConstReturnValue ConstReturnValue;
+      typedef typename MapTraits<In>::ReturnValue Reference;
+      typedef typename MapTraits<In>::ConstReturnValue ConstReference;
 
       /// Constructor
-      CombinedNodeMap(IN& in_map, OUT& out_map)
+      CombinedNodeMap(In& in_map, Out& out_map)
         : _in_map(in_map), _out_map(out_map) {}
 
@@ -3498,6 +3498,6 @@
     private:
 
-      IN& _in_map;
-      OUT& _out_map;
+      In& _in_map;
+      Out& _out_map;
 
     };
@@ -3507,26 +3507,26 @@
     ///
     /// This function just returns a combined node map.
-    template <typename IN, typename OUT>
-    static CombinedNodeMap<IN, OUT>
-    combinedNodeMap(IN& in_map, OUT& out_map) {
-      return CombinedNodeMap<IN, OUT>(in_map, out_map);
-    }
-
-    template <typename IN, typename OUT>
-    static CombinedNodeMap<const IN, OUT>
-    combinedNodeMap(const IN& in_map, OUT& out_map) {
-      return CombinedNodeMap<const IN, OUT>(in_map, out_map);
-    }
-
-    template <typename IN, typename OUT>
-    static CombinedNodeMap<IN, const OUT>
-    combinedNodeMap(IN& in_map, const OUT& out_map) {
-      return CombinedNodeMap<IN, const OUT>(in_map, out_map);
-    }
-
-    template <typename IN, typename OUT>
-    static CombinedNodeMap<const IN, const OUT>
-    combinedNodeMap(const IN& in_map, const OUT& out_map) {
-      return CombinedNodeMap<const IN, const OUT>(in_map, out_map);
+    template <typename In, typename Out>
+    static CombinedNodeMap<In, Out>
+    combinedNodeMap(In& in_map, Out& out_map) {
+      return CombinedNodeMap<In, Out>(in_map, out_map);
+    }
+
+    template <typename In, typename Out>
+    static CombinedNodeMap<const In, Out>
+    combinedNodeMap(const In& in_map, Out& out_map) {
+      return CombinedNodeMap<const In, Out>(in_map, out_map);
+    }
+
+    template <typename In, typename Out>
+    static CombinedNodeMap<In, const Out>
+    combinedNodeMap(In& in_map, const Out& out_map) {
+      return CombinedNodeMap<In, const Out>(in_map, out_map);
+    }
+
+    template <typename In, typename Out>
+    static CombinedNodeMap<const In, const Out>
+    combinedNodeMap(const In& in_map, const Out& out_map) {
+      return CombinedNodeMap<const In, const Out>(in_map, out_map);
     }
 
Index: lemon/arg_parser.cc
===================================================================
--- lemon/arg_parser.cc	(revision 956)
+++ lemon/arg_parser.cc	(revision 1397)
@@ -222,7 +222,7 @@
   {
     Opts::iterator o = _opts.find(opt);
-    Opts::iterator s = _opts.find(syn);
     LEMON_ASSERT(o!=_opts.end(), "Unknown option: '"+opt+"'");
-    LEMON_ASSERT(s==_opts.end(), "Option already used: '"+syn+"'");
+    LEMON_ASSERT(_opts.find(syn)==_opts.end(),
+                 "Option already used: '"+syn+"'");
     ParData p;
     p.help=opt;
Index: lemon/lgf_reader.h
===================================================================
--- lemon/lgf_reader.h	(revision 1361)
+++ lemon/lgf_reader.h	(revision 1377)
@@ -452,11 +452,11 @@
   ///
   ///\code
-  /// DigraphReader<DGR>(digraph, std::cin).
-  ///   nodeMap("coordinates", coord_map).
-  ///   arcMap("capacity", cap_map).
-  ///   node("source", src).
-  ///   node("target", trg).
-  ///   attribute("caption", caption).
-  ///   run();
+  /// DigraphReader<DGR>(digraph, std::cin)
+  ///   .nodeMap("coordinates", coord_map)
+  ///   .arcMap("capacity", cap_map)
+  ///   .node("source", src)
+  ///   .node("target", trg)
+  ///   .attribute("caption", caption)
+  ///   .run();
   ///\endcode
   ///
@@ -1247,9 +1247,9 @@
   ///ListDigraph::ArcMap<int> cap(digraph);
   ///ListDigraph::Node src, trg;
-  ///digraphReader(digraph, std::cin).
-  ///  arcMap("capacity", cap).
-  ///  node("source", src).
-  ///  node("target", trg).
-  ///  run();
+  ///digraphReader(digraph, std::cin)
+  ///  .arcMap("capacity", cap)
+  ///  .node("source", src)
+  ///  .node("target", trg)
+  ///  .run();
   ///\endcode
   ///
@@ -2124,7 +2124,7 @@
   ///ListGraph graph;
   ///ListGraph::EdgeMap<int> weight(graph);
-  ///graphReader(graph, std::cin).
-  ///  edgeMap("weight", weight).
-  ///  run();
+  ///graphReader(graph, std::cin)
+  ///  .edgeMap("weight", weight)
+  ///  .run();
   ///\endcode
   ///
@@ -3192,7 +3192,7 @@
   ///ListBpGraph graph;
   ///ListBpGraph::EdgeMap<int> weight(graph);
-  ///bpGraphReader(graph, std::cin).
-  ///  edgeMap("weight", weight).
-  ///  run();
+  ///bpGraphReader(graph, std::cin)
+  ///  .edgeMap("weight", weight)
+  ///  .run();
   ///\endcode
   ///
Index: lemon/lgf_writer.h
===================================================================
--- lemon/lgf_writer.h	(revision 1270)
+++ lemon/lgf_writer.h	(revision 1377)
@@ -409,13 +409,13 @@
   ///
   ///\code
-  /// DigraphWriter<DGR>(digraph, std::cout).
-  ///   nodeMap("coordinates", coord_map).
-  ///   nodeMap("size", size).
-  ///   nodeMap("title", title).
-  ///   arcMap("capacity", cap_map).
-  ///   node("source", src).
-  ///   node("target", trg).
-  ///   attribute("caption", caption).
-  ///   run();
+  /// DigraphWriter<DGR>(digraph, std::cout)
+  ///   .nodeMap("coordinates", coord_map)
+  ///   .nodeMap("size", size)
+  ///   .nodeMap("title", title)
+  ///   .arcMap("capacity", cap_map)
+  ///   .node("source", src)
+  ///   .node("target", trg)
+  ///   .attribute("caption", caption)
+  ///   .run();
   ///\endcode
   ///
@@ -962,9 +962,9 @@
   ///ListDigraph::Node src, trg;
   ///  // Setting the capacity map and source and target nodes
-  ///digraphWriter(digraph, std::cout).
-  ///  arcMap("capacity", cap).
-  ///  node("source", src).
-  ///  node("target", trg).
-  ///  run();
+  ///digraphWriter(digraph, std::cout)
+  ///  .arcMap("capacity", cap)
+  ///  .node("source", src)
+  ///  .node("target", trg)
+  ///  .run();
   ///\endcode
   ///
@@ -1600,7 +1600,7 @@
   ///ListGraph::EdgeMap<int> weight(graph);
   ///  // Setting the weight map
-  ///graphWriter(graph, std::cout).
-  ///  edgeMap("weight", weight).
-  ///  run();
+  ///graphWriter(graph, std::cout)
+  ///  .edgeMap("weight", weight)
+  ///  .run();
   ///\endcode
   ///
@@ -2420,7 +2420,7 @@
   ///ListBpGraph::EdgeMap<int> weight(graph);
   ///  // Setting the weight map
-  ///bpGraphWriter(graph, std::cout).
-  ///  edgeMap("weight", weight).
-  ///  run();
+  ///bpGraphWriter(graph, std::cout)
+  ///  .edgeMap("weight", weight)
+  ///  .run();
   ///\endcode
   ///
Index: lemon/planarity.h
===================================================================
--- lemon/planarity.h	(revision 1270)
+++ lemon/planarity.h	(revision 1400)
@@ -2384,5 +2384,5 @@
       if (!pe.run()) return false;
 
-      run(pe);
+      run(pe.embeddingMap());
       return true;
     }
@@ -2398,4 +2398,13 @@
     void run(const EmbeddingMap& embedding) {
       typedef SmartEdgeSet<Graph> AuxGraph;
+
+      if (countNodes(_graph) < 3) {
+        int y = 0;
+        for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
+          _point_map[n].x = 0;
+          _point_map[n].y = y++;
+        }
+        return;
+      }
 
       if (3 * countNodes(_graph) - 6 == countEdges(_graph)) {
Index: lemon/preflow.h
===================================================================
--- lemon/preflow.h	(revision 1270)
+++ lemon/preflow.h	(revision 1385)
@@ -477,5 +477,5 @@
     /// flow to the given \c flowMap. The \c flowMap should contain a
     /// flow or at least a preflow, i.e. at each node excluding the
-    /// source node the incoming flow should greater or equal to the
+    /// source node the incoming flow should be greater or equal to the
     /// outgoing flow.
     /// \return \c false if the given \c flowMap is not a preflow.
@@ -496,5 +496,5 @@
           excess -= (*_flow)[e];
         }
-        if (excess < 0 && n != _source) return false;
+        if (_tolerance.negative(excess) && n != _source) return false;
         (*_excess)[n] = excess;
       }
@@ -640,5 +640,5 @@
           (*_excess)[n] = excess;
 
-          if (excess != 0) {
+          if (_tolerance.nonZero(excess)) {
             if (new_level + 1 < _level->maxLevel()) {
               _level->liftHighestActive(new_level + 1);
@@ -721,5 +721,5 @@
           (*_excess)[n] = excess;
 
-          if (excess != 0) {
+          if (_tolerance.nonZero(excess)) {
             if (new_level + 1 < _level->maxLevel()) {
               _level->liftActiveOn(level, new_level + 1);
@@ -792,5 +792,5 @@
         if (!reached[n]) {
           _level->dirtyTopButOne(n);
-        } else if ((*_excess)[n] > 0 && _target != n) {
+        } else if (_tolerance.positive((*_excess)[n]) && _target != n) {
           _level->activate(n);
         }
@@ -853,5 +853,5 @@
         (*_excess)[n] = excess;
 
-        if (excess != 0) {
+        if (_tolerance.nonZero(excess)) {
           if (new_level + 1 < _level->maxLevel()) {
             _level->liftHighestActive(new_level + 1);
Index: lemon/random.h
===================================================================
--- lemon/random.h	(revision 1343)
+++ lemon/random.h	(revision 1404)
@@ -111,5 +111,4 @@
       static const Word loMask = (1u << 31) - 1;
       static const Word hiMask = ~loMask;
-
 
       static Word tempering(Word rnd) {
@@ -243,5 +242,4 @@
 
     private:
-
 
       void fillState() {
@@ -272,5 +270,4 @@
       }
 
-
       Word *current;
       Word state[length];
@@ -297,5 +294,5 @@
     template <typename Result, typename Word,
               int rest = std::numeric_limits<Result>::digits, int shift = 0,
-              bool last = rest <= std::numeric_limits<Word>::digits>
+              bool last = (rest <= std::numeric_limits<Word>::digits)>
     struct IntConversion {
       static const int bits = std::numeric_limits<Word>::digits;
@@ -343,5 +340,5 @@
           num = rnd() & mask;
         } while (num > max);
-        return num;
+        return static_cast<Result>(num);
       }
     };
@@ -466,5 +463,544 @@
     };
 
-  }
+    /// \ingroup misc
+    ///
+    /// \brief Mersenne Twister random number generator
+    ///
+    /// The Mersenne Twister is a twisted generalized feedback
+    /// shift-register generator of Matsumoto and Nishimura. The period
+    /// of this generator is \f$ 2^{19937} - 1\f$ and it is
+    /// equi-distributed in 623 dimensions for 32-bit numbers. The time
+    /// performance of this generator is comparable to the commonly used
+    /// generators.
+    ///
+    /// This is a template implementation of both 32-bit and
+    /// 64-bit architecture optimized versions. The generators differ
+    /// sligthly in the initialization and generation phase so they
+    /// produce two completly different sequences.
+    ///
+    /// \alert Do not use this class directly, but instead one of \ref
+    /// Random, \ref Random32 or \ref Random64.
+    ///
+    /// The generator gives back random numbers of serveral types. To
+    /// get a random number from a range of a floating point type, you
+    /// can use one form of the \c operator() or the \c real() member
+    /// function. If you want to get random number from the {0, 1, ...,
+    /// n-1} integer range, use the \c operator[] or the \c integer()
+    /// method. And to get random number from the whole range of an
+    /// integer type, you can use the argumentless \c integer() or
+    /// \c uinteger() functions. Finally, you can get random bool with
+    /// equal chance of true and false or with given probability of true
+    /// result using the \c boolean() member functions.
+    ///
+    /// Various non-uniform distributions are also supported: normal (Gauss),
+    /// exponential, gamma, Poisson, etc.; and a few two-dimensional
+    /// distributions, too.
+    ///
+    ///\code
+    /// // The commented code is identical to the other
+    /// double a = rnd();                     // [0.0, 1.0)
+    /// // double a = rnd.real();             // [0.0, 1.0)
+    /// double b = rnd(100.0);                // [0.0, 100.0)
+    /// // double b = rnd.real(100.0);        // [0.0, 100.0)
+    /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
+    /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
+    /// int d = rnd[100000];                  // 0..99999
+    /// // int d = rnd.integer(100000);       // 0..99999
+    /// int e = rnd[6] + 1;                   // 1..6
+    /// // int e = rnd.integer(1, 1 + 6);     // 1..6
+    /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
+    /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
+    /// bool g = rnd.boolean();               // P(g = true) = 0.5
+    /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
+    ///\endcode
+    ///
+    /// LEMON provides a global instance of the random number generator:
+    /// \ref lemon::rnd "rnd". In most cases, it is a good practice
+    /// to use this global generator to get random numbers.
+    ///
+    /// \sa \ref Random, \ref Random32 or \ref Random64.
+    template<class Word>
+    class Random {
+    private:
+
+      _random_bits::RandomCore<Word> core;
+      _random_bits::BoolProducer<Word> bool_producer;
+
+
+    public:
+
+      ///\name Initialization
+      ///
+      /// @{
+
+      /// \brief Default constructor
+      ///
+      /// Constructor with constant seeding.
+      Random() { core.initState(); }
+
+      /// \brief Constructor with seed
+      ///
+      /// Constructor with seed. The current number type will be converted
+      /// to the architecture word type.
+      template <typename Number>
+      Random(Number seed) {
+        _random_bits::Initializer<Number, Word>::init(core, seed);
+      }
+
+      /// \brief Constructor with array seeding
+      ///
+      /// Constructor with array seeding. The given range should contain
+      /// any number type and the numbers will be converted to the
+      /// architecture word type.
+      template <typename Iterator>
+      Random(Iterator begin, Iterator end) {
+        typedef typename std::iterator_traits<Iterator>::value_type Number;
+        _random_bits::Initializer<Number, Word>::init(core, begin, end);
+      }
+
+      /// \brief Copy constructor
+      ///
+      /// Copy constructor. The generated sequence will be identical to
+      /// the other sequence. It can be used to save the current state
+      /// of the generator and later use it to generate the same
+      /// sequence.
+      Random(const Random& other) {
+        core.copyState(other.core);
+      }
+
+      /// \brief Assign operator
+      ///
+      /// Assign operator. The generated sequence will be identical to
+      /// the other sequence. It can be used to save the current state
+      /// of the generator and later use it to generate the same
+      /// sequence.
+      Random& operator=(const Random& other) {
+        if (&other != this) {
+          core.copyState(other.core);
+        }
+        return *this;
+      }
+
+      /// \brief Seeding random sequence
+      ///
+      /// Seeding the random sequence. The current number type will be
+      /// converted to the architecture word type.
+      template <typename Number>
+      void seed(Number seed) {
+        _random_bits::Initializer<Number, Word>::init(core, seed);
+      }
+
+      /// \brief Seeding random sequence
+      ///
+      /// Seeding the random sequence. The given range should contain
+      /// any number type and the numbers will be converted to the
+      /// architecture word type.
+      template <typename Iterator>
+      void seed(Iterator begin, Iterator end) {
+        typedef typename std::iterator_traits<Iterator>::value_type Number;
+        _random_bits::Initializer<Number, Word>::init(core, begin, end);
+      }
+
+      /// \brief Seeding from file or from process id and time
+      ///
+      /// By default, this function calls the \c seedFromFile() member
+      /// function with the <tt>/dev/urandom</tt> file. If it does not success,
+      /// it uses the \c seedFromTime().
+      /// \return Currently always \c true.
+      bool seed() {
+#ifndef LEMON_WIN32
+        if (seedFromFile("/dev/urandom", 0)) return true;
+#endif
+        if (seedFromTime()) return true;
+        return false;
+      }
+
+      /// \brief Seeding from file
+      ///
+      /// Seeding the random sequence from file. The linux kernel has two
+      /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
+      /// could give good seed values for pseudo random generators (The
+      /// difference between two devices is that the <tt>random</tt> may
+      /// block the reading operation while the kernel can give good
+      /// source of randomness, while the <tt>urandom</tt> does not
+      /// block the input, but it could give back bytes with worse
+      /// entropy).
+      /// \param file The source file
+      /// \param offset The offset, from the file read.
+      /// \return \c true when the seeding successes.
+#ifndef LEMON_WIN32
+      bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
+#else
+        bool seedFromFile(const std::string& file = "", int offset = 0)
+#endif
+      {
+        std::ifstream rs(file.c_str());
+        const int size = 4;
+        Word buf[size];
+        if (offset != 0 && !rs.seekg(offset)) return false;
+        if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
+        seed(buf, buf + size);
+        return true;
+      }
+
+      /// \brief Seeding from process id and time
+      ///
+      /// Seeding from process id and time. This function uses the
+      /// current process id and the current time for initialize the
+      /// random sequence.
+      /// \return Currently always \c true.
+      bool seedFromTime() {
+#ifndef LEMON_WIN32
+        timeval tv;
+        gettimeofday(&tv, 0);
+        seed(getpid() + tv.tv_sec + tv.tv_usec);
+#else
+        seed(bits::getWinRndSeed());
+#endif
+        return true;
+      }
+
+      /// @}
+
+      ///\name Uniform Distributions
+      ///
+      /// @{
+
+      /// \brief Returns a random real number from the range [0, 1)
+      ///
+      /// It returns a random real number from the range [0, 1). The
+      /// default Number type is \c double.
+      template <typename Number>
+      Number real() {
+        return _random_bits::RealConversion<Number, Word>::convert(core);
+      }
+
+      double real() {
+        return real<double>();
+      }
+
+      /// \brief Returns a random real number from the range [0, 1)
+      ///
+      /// It returns a random double from the range [0, 1).
+      double operator()() {
+        return real<double>();
+      }
+
+      /// \brief Returns a random real number from the range [0, b)
+      ///
+      /// It returns a random real number from the range [0, b).
+      double operator()(double b) {
+        return real<double>() * b;
+      }
+
+      /// \brief Returns a random real number from the range [a, b)
+      ///
+      /// It returns a random real number from the range [a, b).
+      double operator()(double a, double b) {
+        return real<double>() * (b - a) + a;
+      }
+
+      /// \brief Returns a random integer from a range
+      ///
+      /// It returns a random integer from the range {0, 1, ..., b - 1}.
+      template <typename Number>
+      Number integer(Number b) {
+        return _random_bits::Mapping<Number, Word>::map(core, b);
+      }
+
+      /// \brief Returns a random integer from a range
+      ///
+      /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
+      template <typename Number>
+      Number integer(Number a, Number b) {
+        return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
+      }
+
+      /// \brief Returns a random integer from a range
+      ///
+      /// It returns a random integer from the range {0, 1, ..., b - 1}.
+      template <typename Number>
+      Number operator[](Number b) {
+        return _random_bits::Mapping<Number, Word>::map(core, b);
+      }
+
+      /// \brief Returns a random non-negative integer
+      ///
+      /// It returns a random non-negative integer uniformly from the
+      /// whole range of the current \c Number type. The default result
+      /// type of this function is <tt>unsigned int</tt>.
+      template <typename Number>
+      Number uinteger() {
+        return _random_bits::IntConversion<Number, Word>::convert(core);
+      }
+
+      unsigned int uinteger() {
+        return uinteger<unsigned int>();
+      }
+
+      /// \brief Returns a random integer
+      ///
+      /// It returns a random integer uniformly from the whole range of
+      /// the current \c Number type. The default result type of this
+      /// function is \c int.
+      template <typename Number>
+      Number integer() {
+        static const int nb = std::numeric_limits<Number>::digits +
+          (std::numeric_limits<Number>::is_signed ? 1 : 0);
+        return _random_bits::IntConversion<Number, Word, nb>::convert(core);
+      }
+
+      int integer() {
+        return integer<int>();
+      }
+
+      /// \brief Returns a random bool
+      ///
+      /// It returns a random bool. The generator holds a buffer for
+      /// random bits. Every time when it become empty the generator makes
+      /// a new random word and fill the buffer up.
+      bool boolean() {
+        return bool_producer.convert(core);
+      }
+
+      /// @}
+
+      ///\name Non-uniform Distributions
+      ///
+      ///@{
+
+      /// \brief Returns a random bool with given probability of true result.
+      ///
+      /// It returns a random bool with given probability of true result.
+      bool boolean(double p) {
+        return operator()() < p;
+      }
+
+      /// Standard normal (Gauss) distribution
+
+      /// Standard normal (Gauss) distribution.
+      /// \note The Cartesian form of the Box-Muller
+      /// transformation is used to generate a random normal distribution.
+      double gauss()
+      {
+        double V1,V2,S;
+        do {
+          V1=2*real<double>()-1;
+          V2=2*real<double>()-1;
+          S=V1*V1+V2*V2;
+        } while(S>=1);
+        return std::sqrt(-2*std::log(S)/S)*V1;
+      }
+      /// Normal (Gauss) distribution with given mean and standard deviation
+
+      /// Normal (Gauss) distribution with given mean and standard deviation.
+      /// \sa gauss()
+      double gauss(double mean,double std_dev)
+      {
+        return gauss()*std_dev+mean;
+      }
+
+      /// Lognormal distribution
+
+      /// Lognormal distribution. The parameters are the mean and the standard
+      /// deviation of <tt>exp(X)</tt>.
+      ///
+      double lognormal(double n_mean,double n_std_dev)
+      {
+        return std::exp(gauss(n_mean,n_std_dev));
+      }
+      /// Lognormal distribution
+
+      /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
+      /// the mean and the standard deviation of <tt>exp(X)</tt>.
+      ///
+      double lognormal(const std::pair<double,double> &params)
+      {
+        return std::exp(gauss(params.first,params.second));
+      }
+      /// Compute the lognormal parameters from mean and standard deviation
+
+      /// This function computes the lognormal parameters from mean and
+      /// standard deviation. The return value can direcly be passed to
+      /// lognormal().
+      std::pair<double,double> lognormalParamsFromMD(double mean,
+                                                     double std_dev)
+      {
+        double fr=std_dev/mean;
+        fr*=fr;
+        double lg=std::log(1+fr);
+        return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
+      }
+      /// Lognormal distribution with given mean and standard deviation
+
+      /// Lognormal distribution with given mean and standard deviation.
+      ///
+      double lognormalMD(double mean,double std_dev)
+      {
+        return lognormal(lognormalParamsFromMD(mean,std_dev));
+      }
+
+      /// Exponential distribution with given mean
+
+      /// This function generates an exponential distribution random number
+      /// with mean <tt>1/lambda</tt>.
+      ///
+      double exponential(double lambda=1.0)
+      {
+        return -std::log(1.0-real<double>())/lambda;
+      }
+
+      /// Gamma distribution with given integer shape
+
+      /// This function generates a gamma distribution random number.
+      ///
+      ///\param k shape parameter (<tt>k>0</tt> integer)
+      double gamma(int k)
+      {
+        double s = 0;
+        for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
+        return s;
+      }
+
+      /// Gamma distribution with given shape and scale parameter
+
+      /// This function generates a gamma distribution random number.
+      ///
+      ///\param k shape parameter (<tt>k>0</tt>)
+      ///\param theta scale parameter
+      ///
+      double gamma(double k,double theta=1.0)
+      {
+        double xi,nu;
+        const double delta = k-std::floor(k);
+        const double v0=E/(E-delta);
+        do {
+          double V0=1.0-real<double>();
+          double V1=1.0-real<double>();
+          double V2=1.0-real<double>();
+          if(V2<=v0)
+            {
+              xi=std::pow(V1,1.0/delta);
+              nu=V0*std::pow(xi,delta-1.0);
+            }
+          else
+            {
+              xi=1.0-std::log(V1);
+              nu=V0*std::exp(-xi);
+            }
+        } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
+        return theta*(xi+gamma(int(std::floor(k))));
+      }
+
+      /// Weibull distribution
+
+      /// This function generates a Weibull distribution random number.
+      ///
+      ///\param k shape parameter (<tt>k>0</tt>)
+      ///\param lambda scale parameter (<tt>lambda>0</tt>)
+      ///
+      double weibull(double k,double lambda)
+      {
+        return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
+      }
+
+      /// Pareto distribution
+
+      /// This function generates a Pareto distribution random number.
+      ///
+      ///\param k shape parameter (<tt>k>0</tt>)
+      ///\param x_min location parameter (<tt>x_min>0</tt>)
+      ///
+      double pareto(double k,double x_min)
+      {
+        return exponential(gamma(k,1.0/x_min))+x_min;
+      }
+
+      /// Poisson distribution
+
+      /// This function generates a Poisson distribution random number with
+      /// parameter \c lambda.
+      ///
+      /// The probability mass function of this distribusion is
+      /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
+      /// \note The algorithm is taken from the book of Donald E. Knuth titled
+      /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
+      /// return value.
+
+      int poisson(double lambda)
+      {
+        const double l = std::exp(-lambda);
+        int k=0;
+        double p = 1.0;
+        do {
+          k++;
+          p*=real<double>();
+        } while (p>=l);
+        return k-1;
+      }
+
+      ///@}
+
+      ///\name Two-Dimensional Distributions
+      ///
+      ///@{
+
+      /// Uniform distribution on the full unit circle
+
+      /// Uniform distribution on the full unit circle.
+      ///
+      dim2::Point<double> disc()
+      {
+        double V1,V2;
+        do {
+          V1=2*real<double>()-1;
+          V2=2*real<double>()-1;
+
+        } while(V1*V1+V2*V2>=1);
+        return dim2::Point<double>(V1,V2);
+      }
+      /// A kind of two-dimensional normal (Gauss) distribution
+
+      /// This function provides a turning symmetric two-dimensional distribution.
+      /// Both coordinates are of standard normal distribution, but they are not
+      /// independent.
+      ///
+      /// \note The coordinates are the two random variables provided by
+      /// the Box-Muller method.
+      dim2::Point<double> gauss2()
+      {
+        double V1,V2,S;
+        do {
+          V1=2*real<double>()-1;
+          V2=2*real<double>()-1;
+          S=V1*V1+V2*V2;
+        } while(S>=1);
+        double W=std::sqrt(-2*std::log(S)/S);
+        return dim2::Point<double>(W*V1,W*V2);
+      }
+      /// A kind of two-dimensional exponential distribution
+
+      /// This function provides a turning symmetric two-dimensional distribution.
+      /// The x-coordinate is of conditionally exponential distribution
+      /// with the condition that x is positive and y=0. If x is negative and
+      /// y=0 then, -x is of exponential distribution. The same is true for the
+      /// y-coordinate.
+      dim2::Point<double> exponential2()
+      {
+        double V1,V2,S;
+        do {
+          V1=2*real<double>()-1;
+          V2=2*real<double>()-1;
+          S=V1*V1+V2*V2;
+        } while(S>=1);
+        double W=-std::log(S)/S;
+        return dim2::Point<double>(W*V1,W*V2);
+      }
+
+      ///@}
+    };
+
+
+  };
 
   /// \ingroup misc
@@ -472,536 +1008,48 @@
   /// \brief Mersenne Twister random number generator
   ///
-  /// The Mersenne Twister is a twisted generalized feedback
-  /// shift-register generator of Matsumoto and Nishimura. The period
-  /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
-  /// equi-distributed in 623 dimensions for 32-bit numbers. The time
-  /// performance of this generator is comparable to the commonly used
-  /// generators.
+  /// This class implements either the 32-bit or the 64-bit version of
+  /// the Mersenne Twister random number generator algorithm
+  /// depending on the system architecture.
+  /// 
+  /// For the API description, see its base class
+  /// \ref _random_bits::Random.
   ///
-  /// This implementation is specialized for both 32-bit and 64-bit
-  /// architectures. The generators differ sligthly in the
-  /// initialization and generation phase so they produce two
-  /// completly different sequences.
+  /// \sa \ref _random_bits::Random
+  typedef _random_bits::Random<unsigned long> Random;
+
+  /// \ingroup misc
   ///
-  /// The generator gives back random numbers of serveral types. To
-  /// get a random number from a range of a floating point type you
-  /// can use one form of the \c operator() or the \c real() member
-  /// function. If you want to get random number from the {0, 1, ...,
-  /// n-1} integer range use the \c operator[] or the \c integer()
-  /// method. And to get random number from the whole range of an
-  /// integer type you can use the argumentless \c integer() or \c
-  /// uinteger() functions. After all you can get random bool with
-  /// equal chance of true and false or given probability of true
-  /// result with the \c boolean() member functions.
+  /// \brief Mersenne Twister random number generator (32-bit version)
   ///
-  ///\code
-  /// // The commented code is identical to the other
-  /// double a = rnd();                     // [0.0, 1.0)
-  /// // double a = rnd.real();             // [0.0, 1.0)
-  /// double b = rnd(100.0);                // [0.0, 100.0)
-  /// // double b = rnd.real(100.0);        // [0.0, 100.0)
-  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
-  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
-  /// int d = rnd[100000];                  // 0..99999
-  /// // int d = rnd.integer(100000);       // 0..99999
-  /// int e = rnd[6] + 1;                   // 1..6
-  /// // int e = rnd.integer(1, 1 + 6);     // 1..6
-  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
-  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
-  /// bool g = rnd.boolean();               // P(g = true) = 0.5
-  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
-  ///\endcode
+  /// This class implements the 32-bit version of the Mersenne Twister
+  /// random number generator algorithm. It is recommended to be used
+  /// when someone wants to make sure that the \e same pseudo random
+  /// sequence will be generated on every platfrom.
   ///
-  /// LEMON provides a global instance of the random number
-  /// generator which name is \ref lemon::rnd "rnd". Usually it is a
-  /// good programming convenience to use this global generator to get
-  /// random numbers.
-  class Random {
-  private:
-
-    // Architecture word
-    typedef unsigned long Word;
-
-    _random_bits::RandomCore<Word> core;
-    _random_bits::BoolProducer<Word> bool_producer;
-
-
-  public:
-
-    ///\name Initialization
-    ///
-    /// @{
-
-    /// \brief Default constructor
-    ///
-    /// Constructor with constant seeding.
-    Random() { core.initState(); }
-
-    /// \brief Constructor with seed
-    ///
-    /// Constructor with seed. The current number type will be converted
-    /// to the architecture word type.
-    template <typename Number>
-    Random(Number seed) {
-      _random_bits::Initializer<Number, Word>::init(core, seed);
-    }
-
-    /// \brief Constructor with array seeding
-    ///
-    /// Constructor with array seeding. The given range should contain
-    /// any number type and the numbers will be converted to the
-    /// architecture word type.
-    template <typename Iterator>
-    Random(Iterator begin, Iterator end) {
-      typedef typename std::iterator_traits<Iterator>::value_type Number;
-      _random_bits::Initializer<Number, Word>::init(core, begin, end);
-    }
-
-    /// \brief Copy constructor
-    ///
-    /// Copy constructor. The generated sequence will be identical to
-    /// the other sequence. It can be used to save the current state
-    /// of the generator and later use it to generate the same
-    /// sequence.
-    Random(const Random& other) {
-      core.copyState(other.core);
-    }
-
-    /// \brief Assign operator
-    ///
-    /// Assign operator. The generated sequence will be identical to
-    /// the other sequence. It can be used to save the current state
-    /// of the generator and later use it to generate the same
-    /// sequence.
-    Random& operator=(const Random& other) {
-      if (&other != this) {
-        core.copyState(other.core);
-      }
-      return *this;
-    }
-
-    /// \brief Seeding random sequence
-    ///
-    /// Seeding the random sequence. The current number type will be
-    /// converted to the architecture word type.
-    template <typename Number>
-    void seed(Number seed) {
-      _random_bits::Initializer<Number, Word>::init(core, seed);
-    }
-
-    /// \brief Seeding random sequence
-    ///
-    /// Seeding the random sequence. The given range should contain
-    /// any number type and the numbers will be converted to the
-    /// architecture word type.
-    template <typename Iterator>
-    void seed(Iterator begin, Iterator end) {
-      typedef typename std::iterator_traits<Iterator>::value_type Number;
-      _random_bits::Initializer<Number, Word>::init(core, begin, end);
-    }
-
-    /// \brief Seeding from file or from process id and time
-    ///
-    /// By default, this function calls the \c seedFromFile() member
-    /// function with the <tt>/dev/urandom</tt> file. If it does not success,
-    /// it uses the \c seedFromTime().
-    /// \return Currently always \c true.
-    bool seed() {
-#ifndef LEMON_WIN32
-      if (seedFromFile("/dev/urandom", 0)) return true;
+  /// For the API description, see its base class
+  /// \ref _random_bits::Random.
+  ///
+  /// \sa \ref _random_bits::Random
+  typedef _random_bits::Random<unsigned int> Random32;
+
+  /// \ingroup misc
+  ///
+  /// \brief Mersenne Twister random number generator (64-bit version)
+  ///
+  /// This class implements the 64-bit version of the Mersenne Twister
+  /// random number generator algorithm. (Even though it runs
+  /// on 32-bit architectures, too.) It is recommended to be used when
+  /// someone wants to make sure that the \e same pseudo random sequence
+  /// will be generated on every platfrom.
+  ///
+  /// For the API description, see its base class
+  /// \ref _random_bits::Random.
+  ///
+  /// \sa \ref _random_bits::Random
+  typedef _random_bits::Random<unsigned long long> Random64;
+
+  extern Random rnd;
+  
+}
+
 #endif
-      if (seedFromTime()) return true;
-      return false;
-    }
-
-    /// \brief Seeding from file
-    ///
-    /// Seeding the random sequence from file. The linux kernel has two
-    /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which
-    /// could give good seed values for pseudo random generators (The
-    /// difference between two devices is that the <tt>random</tt> may
-    /// block the reading operation while the kernel can give good
-    /// source of randomness, while the <tt>urandom</tt> does not
-    /// block the input, but it could give back bytes with worse
-    /// entropy).
-    /// \param file The source file
-    /// \param offset The offset, from the file read.
-    /// \return \c true when the seeding successes.
-#ifndef LEMON_WIN32
-    bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
-#else
-    bool seedFromFile(const std::string& file = "", int offset = 0)
-#endif
-    {
-      std::ifstream rs(file.c_str());
-      const int size = 4;
-      Word buf[size];
-      if (offset != 0 && !rs.seekg(offset)) return false;
-      if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;
-      seed(buf, buf + size);
-      return true;
-    }
-
-    /// \brief Seding from process id and time
-    ///
-    /// Seding from process id and time. This function uses the
-    /// current process id and the current time for initialize the
-    /// random sequence.
-    /// \return Currently always \c true.
-    bool seedFromTime() {
-#ifndef LEMON_WIN32
-      timeval tv;
-      gettimeofday(&tv, 0);
-      seed(getpid() + tv.tv_sec + tv.tv_usec);
-#else
-      seed(bits::getWinRndSeed());
-#endif
-      return true;
-    }
-
-    /// @}
-
-    ///\name Uniform Distributions
-    ///
-    /// @{
-
-    /// \brief Returns a random real number from the range [0, 1)
-    ///
-    /// It returns a random real number from the range [0, 1). The
-    /// default Number type is \c double.
-    template <typename Number>
-    Number real() {
-      return _random_bits::RealConversion<Number, Word>::convert(core);
-    }
-
-    double real() {
-      return real<double>();
-    }
-
-    /// \brief Returns a random real number from the range [0, 1)
-    ///
-    /// It returns a random double from the range [0, 1).
-    double operator()() {
-      return real<double>();
-    }
-
-    /// \brief Returns a random real number from the range [0, b)
-    ///
-    /// It returns a random real number from the range [0, b).
-    double operator()(double b) {
-      return real<double>() * b;
-    }
-
-    /// \brief Returns a random real number from the range [a, b)
-    ///
-    /// It returns a random real number from the range [a, b).
-    double operator()(double a, double b) {
-      return real<double>() * (b - a) + a;
-    }
-
-    /// \brief Returns a random integer from a range
-    ///
-    /// It returns a random integer from the range {0, 1, ..., b - 1}.
-    template <typename Number>
-    Number integer(Number b) {
-      return _random_bits::Mapping<Number, Word>::map(core, b);
-    }
-
-    /// \brief Returns a random integer from a range
-    ///
-    /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
-    template <typename Number>
-    Number integer(Number a, Number b) {
-      return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
-    }
-
-    /// \brief Returns a random integer from a range
-    ///
-    /// It returns a random integer from the range {0, 1, ..., b - 1}.
-    template <typename Number>
-    Number operator[](Number b) {
-      return _random_bits::Mapping<Number, Word>::map(core, b);
-    }
-
-    /// \brief Returns a random non-negative integer
-    ///
-    /// It returns a random non-negative integer uniformly from the
-    /// whole range of the current \c Number type. The default result
-    /// type of this function is <tt>unsigned int</tt>.
-    template <typename Number>
-    Number uinteger() {
-      return _random_bits::IntConversion<Number, Word>::convert(core);
-    }
-
-    unsigned int uinteger() {
-      return uinteger<unsigned int>();
-    }
-
-    /// \brief Returns a random integer
-    ///
-    /// It returns a random integer uniformly from the whole range of
-    /// the current \c Number type. The default result type of this
-    /// function is \c int.
-    template <typename Number>
-    Number integer() {
-      static const int nb = std::numeric_limits<Number>::digits +
-        (std::numeric_limits<Number>::is_signed ? 1 : 0);
-      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
-    }
-
-    int integer() {
-      return integer<int>();
-    }
-
-    /// \brief Returns a random bool
-    ///
-    /// It returns a random bool. The generator holds a buffer for
-    /// random bits. Every time when it become empty the generator makes
-    /// a new random word and fill the buffer up.
-    bool boolean() {
-      return bool_producer.convert(core);
-    }
-
-    /// @}
-
-    ///\name Non-uniform Distributions
-    ///
-    ///@{
-
-    /// \brief Returns a random bool with given probability of true result.
-    ///
-    /// It returns a random bool with given probability of true result.
-    bool boolean(double p) {
-      return operator()() < p;
-    }
-
-    /// Standard normal (Gauss) distribution
-
-    /// Standard normal (Gauss) distribution.
-    /// \note The Cartesian form of the Box-Muller
-    /// transformation is used to generate a random normal distribution.
-    double gauss()
-    {
-      double V1,V2,S;
-      do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        S=V1*V1+V2*V2;
-      } while(S>=1);
-      return std::sqrt(-2*std::log(S)/S)*V1;
-    }
-    /// Normal (Gauss) distribution with given mean and standard deviation
-
-    /// Normal (Gauss) distribution with given mean and standard deviation.
-    /// \sa gauss()
-    double gauss(double mean,double std_dev)
-    {
-      return gauss()*std_dev+mean;
-    }
-
-    /// Lognormal distribution
-
-    /// Lognormal distribution. The parameters are the mean and the standard
-    /// deviation of <tt>exp(X)</tt>.
-    ///
-    double lognormal(double n_mean,double n_std_dev)
-    {
-      return std::exp(gauss(n_mean,n_std_dev));
-    }
-    /// Lognormal distribution
-
-    /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of
-    /// the mean and the standard deviation of <tt>exp(X)</tt>.
-    ///
-    double lognormal(const std::pair<double,double> &params)
-    {
-      return std::exp(gauss(params.first,params.second));
-    }
-    /// Compute the lognormal parameters from mean and standard deviation
-
-    /// This function computes the lognormal parameters from mean and
-    /// standard deviation. The return value can direcly be passed to
-    /// lognormal().
-    std::pair<double,double> lognormalParamsFromMD(double mean,
-                                                   double std_dev)
-    {
-      double fr=std_dev/mean;
-      fr*=fr;
-      double lg=std::log(1+fr);
-      return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));
-    }
-    /// Lognormal distribution with given mean and standard deviation
-
-    /// Lognormal distribution with given mean and standard deviation.
-    ///
-    double lognormalMD(double mean,double std_dev)
-    {
-      return lognormal(lognormalParamsFromMD(mean,std_dev));
-    }
-
-    /// Exponential distribution with given mean
-
-    /// This function generates an exponential distribution random number
-    /// with mean <tt>1/lambda</tt>.
-    ///
-    double exponential(double lambda=1.0)
-    {
-      return -std::log(1.0-real<double>())/lambda;
-    }
-
-    /// Gamma distribution with given integer shape
-
-    /// This function generates a gamma distribution random number.
-    ///
-    ///\param k shape parameter (<tt>k>0</tt> integer)
-    double gamma(int k)
-    {
-      double s = 0;
-      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
-      return s;
-    }
-
-    /// Gamma distribution with given shape and scale parameter
-
-    /// This function generates a gamma distribution random number.
-    ///
-    ///\param k shape parameter (<tt>k>0</tt>)
-    ///\param theta scale parameter
-    ///
-    double gamma(double k,double theta=1.0)
-    {
-      double xi,nu;
-      const double delta = k-std::floor(k);
-      const double v0=E/(E-delta);
-      do {
-        double V0=1.0-real<double>();
-        double V1=1.0-real<double>();
-        double V2=1.0-real<double>();
-        if(V2<=v0)
-          {
-            xi=std::pow(V1,1.0/delta);
-            nu=V0*std::pow(xi,delta-1.0);
-          }
-        else
-          {
-            xi=1.0-std::log(V1);
-            nu=V0*std::exp(-xi);
-          }
-      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
-      return theta*(xi+gamma(int(std::floor(k))));
-    }
-
-    /// Weibull distribution
-
-    /// This function generates a Weibull distribution random number.
-    ///
-    ///\param k shape parameter (<tt>k>0</tt>)
-    ///\param lambda scale parameter (<tt>lambda>0</tt>)
-    ///
-    double weibull(double k,double lambda)
-    {
-      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
-    }
-
-    /// Pareto distribution
-
-    /// This function generates a Pareto distribution random number.
-    ///
-    ///\param k shape parameter (<tt>k>0</tt>)
-    ///\param x_min location parameter (<tt>x_min>0</tt>)
-    ///
-    double pareto(double k,double x_min)
-    {
-      return exponential(gamma(k,1.0/x_min))+x_min;
-    }
-
-    /// Poisson distribution
-
-    /// This function generates a Poisson distribution random number with
-    /// parameter \c lambda.
-    ///
-    /// The probability mass function of this distribusion is
-    /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
-    /// \note The algorithm is taken from the book of Donald E. Knuth titled
-    /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
-    /// return value.
-
-    int poisson(double lambda)
-    {
-      const double l = std::exp(-lambda);
-      int k=0;
-      double p = 1.0;
-      do {
-        k++;
-        p*=real<double>();
-      } while (p>=l);
-      return k-1;
-    }
-
-    ///@}
-
-    ///\name Two Dimensional Distributions
-    ///
-    ///@{
-
-    /// Uniform distribution on the full unit circle
-
-    /// Uniform distribution on the full unit circle.
-    ///
-    dim2::Point<double> disc()
-    {
-      double V1,V2;
-      do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-
-      } while(V1*V1+V2*V2>=1);
-      return dim2::Point<double>(V1,V2);
-    }
-    /// A kind of two dimensional normal (Gauss) distribution
-
-    /// This function provides a turning symmetric two-dimensional distribution.
-    /// Both coordinates are of standard normal distribution, but they are not
-    /// independent.
-    ///
-    /// \note The coordinates are the two random variables provided by
-    /// the Box-Muller method.
-    dim2::Point<double> gauss2()
-    {
-      double V1,V2,S;
-      do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        S=V1*V1+V2*V2;
-      } while(S>=1);
-      double W=std::sqrt(-2*std::log(S)/S);
-      return dim2::Point<double>(W*V1,W*V2);
-    }
-    /// A kind of two dimensional exponential distribution
-
-    /// This function provides a turning symmetric two-dimensional distribution.
-    /// The x-coordinate is of conditionally exponential distribution
-    /// with the condition that x is positive and y=0. If x is negative and
-    /// y=0 then, -x is of exponential distribution. The same is true for the
-    /// y-coordinate.
-    dim2::Point<double> exponential2()
-    {
-      double V1,V2,S;
-      do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        S=V1*V1+V2*V2;
-      } while(S>=1);
-      double W=-std::log(S)/S;
-      return dim2::Point<double>(W*V1,W*V2);
-    }
-
-    ///@}
-  };
-
-
-  extern Random rnd;
-
-}
-
-#endif
Index: test/bellman_ford_test.cc
===================================================================
--- test/bellman_ford_test.cc	(revision 1337)
+++ test/bellman_ford_test.cc	(revision 1378)
@@ -101,4 +101,5 @@
     pp = const_bf_test.negativeCycle();
 
+#ifdef LEMON_CXX11
     for (BF::ActiveIt it(const_bf_test); it != INVALID; ++it) {}
     for (auto n: const_bf_test.activeNodes()) { ::lemon::ignore_unused_variable_warning(n); }
@@ -106,4 +107,5 @@
       ::lemon::ignore_unused_variable_warning(n);
     }
+#endif
   }
   {
Index: test/max_flow_test.cc
===================================================================
--- test/max_flow_test.cc	(revision 1270)
+++ test/max_flow_test.cc	(revision 1396)
@@ -27,4 +27,5 @@
 #include <lemon/lgf_reader.h>
 #include <lemon/elevator.h>
+#include <lemon/tolerance.h>
 
 using namespace lemon;
@@ -66,4 +67,35 @@
   "target 8\n";
 
+char test_lgf_float[] =
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "8\n"
+  "9\n"
+  "@arcs\n"
+  "      capacity\n"
+  "0 1 0.1\n"
+  "0 2 0.1\n"
+  "0 3 0.1\n"
+  "1 4 0.1\n"
+  "2 4 0.1\n"
+  "3 4 0.1\n"
+  "4 5 0.3\n"
+  "5 6 0.1\n"
+  "5 7 0.1\n"
+  "5 8 0.1\n"
+  "6 9 0.1\n"
+  "7 9 0.1\n"
+  "8 9 0.1\n"
+  "@attributes\n"
+  "source 0\n"
+  "target 9\n";
 
 // Checks the general interface of a max flow algorithm
@@ -166,6 +198,4 @@
   typedef concepts::Digraph Digraph;
   typedef concepts::ReadMap<Digraph::Arc, Value> CapMap;
-  typedef Elevator<Digraph, Digraph::Node> Elev;
-  typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;
 
   Digraph g;
@@ -185,11 +215,11 @@
 
 template <typename T>
-T cutValue (const SmartDigraph& g,
-              const SmartDigraph::NodeMap<bool>& cut,
-              const SmartDigraph::ArcMap<T>& cap) {
-
-  T c=0;
-  for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {
-    if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];
+T cutValue(const SmartDigraph& g,
+           const SmartDigraph::NodeMap<bool>& cut,
+           const SmartDigraph::ArcMap<T>& cap) {
+
+  T c = 0;
+  for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
+    if (cut[g.source(e)] && !cut[g.target(e)]) c += cap[e];
   }
   return c;
@@ -200,8 +230,9 @@
                const SmartDigraph::ArcMap<T>& flow,
                const SmartDigraph::ArcMap<T>& cap,
-               SmartDigraph::Node s, SmartDigraph::Node t) {
+               SmartDigraph::Node s, SmartDigraph::Node t,
+               const Tolerance<T>& tol) {
 
   for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
-    if (flow[e] < 0 || flow[e] > cap[e]) return false;
+    if (tol.negative(flow[e]) || tol.less(cap[e], flow[e])) return false;
   }
 
@@ -215,26 +246,27 @@
       sum -= flow[e];
     }
-    if (sum != 0) return false;
+    if (tol.nonZero(sum)) return false;
   }
   return true;
 }
 
-void initFlowTest()
+void checkInitPreflow()
 {
   DIGRAPH_TYPEDEFS(SmartDigraph);
 
   SmartDigraph g;
-  SmartDigraph::ArcMap<int> cap(g),iflow(g);
-  Node s=g.addNode(); Node t=g.addNode();
-  Node n1=g.addNode(); Node n2=g.addNode();
+  SmartDigraph::ArcMap<int> cap(g), iflow(g);
+  Node s = g.addNode(); Node t = g.addNode();
+  Node n1 = g.addNode(); Node n2 = g.addNode();
   Arc a;
-  a=g.addArc(s,n1); cap[a]=20; iflow[a]=20;
-  a=g.addArc(n1,n2); cap[a]=10; iflow[a]=0;
-  a=g.addArc(n2,t); cap[a]=20; iflow[a]=0;
-
-  Preflow<SmartDigraph> pre(g,cap,s,t);
+  a = g.addArc(s, n1); cap[a] = 20; iflow[a] = 20;
+  a = g.addArc(n1, n2); cap[a] = 10; iflow[a] = 0;
+  a = g.addArc(n2, t); cap[a] = 20; iflow[a] = 0;
+
+  Preflow<SmartDigraph> pre(g, cap, s, t);
   pre.init(iflow);
   pre.startFirstPhase();
-  check(pre.flowValue() == 10, "The incorrect max flow value.");
+
+  check(pre.flowValue() == 10, "Incorrect max flow value.");
   check(pre.minCut(s), "Wrong min cut (Node s).");
   check(pre.minCut(n1), "Wrong min cut (Node n1).");
@@ -244,5 +276,5 @@
 
 template <typename MF, typename SF>
-void checkMaxFlowAlg() {
+void checkMaxFlowAlg(const char *input_lgf,  typename MF::Value expected) {
   typedef SmartDigraph Digraph;
   DIGRAPH_TYPEDEFS(Digraph);
@@ -253,12 +285,14 @@
   typedef BoolNodeMap CutMap;
 
+  Tolerance<Value> tol;
+
   Digraph g;
   Node s, t;
   CapMap cap(g);
-  std::istringstream input(test_lgf);
-  DigraphReader<Digraph>(g,input)
+  std::istringstream input(input_lgf);
+  DigraphReader<Digraph>(g, input)
       .arcMap("capacity", cap)
-      .node("source",s)
-      .node("target",t)
+      .node("source", s)
+      .node("target", t)
       .run();
 
@@ -266,5 +300,7 @@
   max_flow.run();
 
-  check(checkFlow(g, max_flow.flowMap(), cap, s, t),
+  check(!tol.different(expected, max_flow.flowValue()),
+        "Incorrect max flow value.");
+  check(checkFlow(g, max_flow.flowMap(), cap, s, t, tol),
         "The flow is not feasible.");
 
@@ -273,13 +309,10 @@
   Value min_cut_value = cutValue(g, min_cut, cap);
 
-  check(max_flow.flowValue() == min_cut_value,
-        "The max flow value is not equal to the min cut value.");
+  check(!tol.different(expected, min_cut_value),
+        "Incorrect min cut value.");
 
   FlowMap flow(g);
-  for (ArcIt e(g); e != INVALID; ++e) flow[e] = max_flow.flowMap()[e];
-
-  Value flow_value = max_flow.flowValue();
-
-  for (ArcIt e(g); e != INVALID; ++e) cap[e] = 2 * cap[e];
+  for (ArcIt e(g); e != INVALID; ++e) flow[e] = 13 * max_flow.flowMap()[e];
+  for (ArcIt e(g); e != INVALID; ++e) cap[e] = 17 * cap[e];
   max_flow.init(flow);
 
@@ -290,11 +323,12 @@
   min_cut_value = cutValue(g, min_cut1, cap);
 
-  check(max_flow.flowValue() == min_cut_value &&
-        min_cut_value == 2 * flow_value,
-        "The max flow value or the min cut value is wrong.");
+  check(!tol.different(17 * expected, max_flow.flowValue()),
+        "Incorrect max flow value.");
+  check(!tol.different(17 * expected, min_cut_value),
+        "Incorrect min cut value.");
 
   SF::startSecondPhase(max_flow);       // start second phase of the algorithm
 
-  check(checkFlow(g, max_flow.flowMap(), cap, s, t),
+  check(checkFlow(g, max_flow.flowMap(), cap, s, t, tol),
         "The flow is not feasible.");
 
@@ -303,8 +337,8 @@
   min_cut_value = cutValue(g, min_cut2, cap);
 
-  check(max_flow.flowValue() == min_cut_value &&
-        min_cut_value == 2 * flow_value,
-        "The max flow value or the min cut value was not doubled");
-
+  check(!tol.different(17 * expected, max_flow.flowValue()),
+        "Incorrect max flow value.");
+  check(!tol.different(17 * expected, min_cut_value),
+        "Incorrect min cut value.");
 
   max_flow.flowMap(flow);
@@ -325,7 +359,7 @@
   CutMap min_cut3(g);
   max_flow.minCutMap(min_cut3);
-  min_cut_value=cutValue(g, min_cut3, cap);
-
-  check(max_flow.flowValue() == min_cut_value,
+  min_cut_value = cutValue(g, min_cut3, cap);
+
+  check(!tol.different(max_flow.flowValue(), min_cut_value),
         "The max flow value or the min cut value is wrong.");
 }
@@ -380,15 +414,26 @@
   typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<int> > PType1;
   typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<float> > PType2;
-  checkMaxFlowAlg<PType1, PreflowStartFunctions<PType1> >();
-  checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >();
-  initFlowTest();
+  typedef Preflow<SmartDigraph, SmartDigraph::ArcMap<double> > PType3;
+
+  checkMaxFlowAlg<PType1, PreflowStartFunctions<PType1> >(test_lgf, 13);
+  checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(test_lgf, 13);
+  checkMaxFlowAlg<PType3, PreflowStartFunctions<PType3> >(test_lgf, 13);
+
+  checkMaxFlowAlg<PType2, PreflowStartFunctions<PType2> >(test_lgf_float, 0.3f);
+  checkMaxFlowAlg<PType3, PreflowStartFunctions<PType3> >(test_lgf_float, 0.3);
+
+  checkInitPreflow();
 
   // Check EdmondsKarp
   typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<int> > EKType1;
   typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<float> > EKType2;
-  checkMaxFlowAlg<EKType1, GeneralStartFunctions<EKType1> >();
-  checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >();
-
-  initFlowTest();
+  typedef EdmondsKarp<SmartDigraph, SmartDigraph::ArcMap<double> > EKType3;
+
+  checkMaxFlowAlg<EKType1, GeneralStartFunctions<EKType1> >(test_lgf, 13);
+  checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(test_lgf, 13);
+  checkMaxFlowAlg<EKType3, GeneralStartFunctions<EKType3> >(test_lgf, 13);
+
+  checkMaxFlowAlg<EKType2, GeneralStartFunctions<EKType2> >(test_lgf_float, 0.3f);
+  checkMaxFlowAlg<EKType3, GeneralStartFunctions<EKType3> >(test_lgf_float, 0.3);
 
   return 0;
Index: test/planarity_test.cc
===================================================================
--- test/planarity_test.cc	(revision 862)
+++ test/planarity_test.cc	(revision 1400)
@@ -31,6 +31,36 @@
 using namespace lemon::dim2;
 
-const int lgfn = 4;
+const int lgfn = 8;
 const std::string lgf[lgfn] = {
+  "@nodes\n"
+  "label\n"
+  "@edges\n"
+  "     label\n",
+
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "@edges\n"
+  "     label\n",
+
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "@edges\n"
+  "     label\n"
+  "0 1  0\n",
+
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "@edges\n"
+  "     label\n"
+  "0 1  0\n"
+  "1 2  1\n"
+  "2 0  2\n",
+
   "@nodes\n"
   "label\n"
@@ -137,6 +167,9 @@
     }
   }
-  check(face_num + countNodes(graph) - countConnectedComponents(graph) ==
-        countEdges(graph) + 1, "Euler test does not passed");
+
+  if (face_num != 0) {
+    check(face_num + countNodes(graph) - countConnectedComponents(graph) ==
+          countEdges(graph) + 1, "Euler test does not passed");
+  }
 }
 
@@ -246,11 +279,27 @@
       checkEmbedding(graph, pe);
 
-      PlanarDrawing<Graph> pd(graph);
-      pd.run(pe.embeddingMap());
-      checkDrawing(graph, pd);
-
-      PlanarColoring<Graph> pc(graph);
-      pc.runFiveColoring(pe.embeddingMap());
-      checkColoring(graph, pc, 5);
+      {
+        PlanarDrawing<Graph> pd(graph);
+        pd.run(pe.embeddingMap());
+        checkDrawing(graph, pd);
+      }
+
+      {
+        PlanarDrawing<Graph> pd(graph);
+        pd.run();
+        checkDrawing(graph, pd);
+      }
+
+      {
+        PlanarColoring<Graph> pc(graph);
+        pc.runFiveColoring(pe.embeddingMap());
+        checkColoring(graph, pc, 5);
+      }
+
+      {
+        PlanarColoring<Graph> pc(graph);
+        pc.runFiveColoring();
+        checkColoring(graph, pc, 5);
+      }
 
     } else {
Index: test/random_test.cc
===================================================================
--- test/random_test.cc	(revision 463)
+++ test/random_test.cc	(revision 1379)
@@ -22,4 +22,31 @@
 int seed_array[] = {1, 2};
 
+int rnd_seq32[] = {
+2732, 43567, 42613, 52416, 45891, 21243, 30403, 32103, 
+62501, 33003, 12172, 5192, 32511, 50057, 43723, 7813, 
+23720, 35343, 6637, 30280, 44566, 31019, 18898, 33867, 
+5994, 1688, 11513, 59011, 48056, 25544, 39168, 25365, 
+17530, 8366, 27063, 49861, 55169, 63848, 11863, 49608
+};
+int rnd_seq64[] = {
+56382, 63883, 59577, 64750, 9644, 59886, 57647, 18152, 
+28520, 64078, 17818, 49294, 26424, 26697, 53684, 19209, 
+35404, 12121, 12837, 11827, 32156, 58333, 62553, 7907, 
+64427, 39399, 21971, 48789, 46981, 15716, 53335, 65256, 
+12999, 15308, 10906, 42162, 47587, 43006, 53921, 18716
+};
+
+void seq_test() {
+  for(int i=0;i<5;i++) {
+    lemon::Random32 r32(i);
+    lemon::Random64 r64(i);
+    for(int j=0;j<8;j++) {
+      check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence");
+      check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence");
+    }
+  }
+}
+
+
 int main()
 {
@@ -37,4 +64,5 @@
                   (sizeof(seed_array) / sizeof(seed_array[0])));
 
+  seq_test();
   return 0;
 }
Index: tools/dimacs-solver.cc
===================================================================
--- tools/dimacs-solver.cc	(revision 1271)
+++ tools/dimacs-solver.cc	(revision 1386)
@@ -223,4 +223,5 @@
         throw IoError("Cannot open the file for writing", ap.files()[1]);
       }
+      // fall through
     case 1:
       input.open(ap.files()[0].c_str());
@@ -228,4 +229,5 @@
         throw IoError("File cannot be found", ap.files()[0]);
       }
+      // fall through
     case 0:
       break;
@@ -252,4 +254,5 @@
         case DimacsDescriptor::SP:
           std::cout << "sp";
+          break;
         case DimacsDescriptor::MAT:
           std::cout << "mat";
Index: tools/dimacs-to-lgf.cc
===================================================================
--- tools/dimacs-to-lgf.cc	(revision 631)
+++ tools/dimacs-to-lgf.cc	(revision 1397)
@@ -74,4 +74,5 @@
         throw IoError("Cannot open the file for writing", ap.files()[1]);
       }
+      // fall through
     case 1:
       input.open(ap.files()[0].c_str());
@@ -79,4 +80,5 @@
         throw IoError("File cannot be found", ap.files()[0]);
       }
+      // fall through
     case 0:
       break;
