New random interface
authordeba
Sat, 14 Oct 2006 15:26:05 +0000
changeset 224216523135943d
parent 2241 37e0966e43b6
child 2243 5deb7b22a0ec
New random interface
Switching to the new interface
benchmark/edge_lookup.cc
benchmark/radix_sort-bench.cc
benchmark/swap_bipartite_bench.cc
demo/descriptor_map_demo.cc
demo/simann_maxcut_demo.cc
demo/topology_demo.cc
lemon/hypercube_graph.h
lemon/random.h
lemon/simann.h
test/all_pairs_shortest_path_test.cc
test/arborescence_test.cc
test/graph_utils_test.cc
test/graph_utils_test.h
test/heap_test.h
test/matrix_maps_test.cc
test/radix_sort_test.cc
test/simann_test.cc
test/test_tools.h
     1.1 --- a/benchmark/edge_lookup.cc	Fri Oct 13 15:10:50 2006 +0000
     1.2 +++ b/benchmark/edge_lookup.cc	Sat Oct 14 15:26:05 2006 +0000
     1.3 @@ -2,6 +2,7 @@
     1.4  #include<lemon/smart_graph.h>
     1.5  #include<vector>
     1.6  #include<lemon/time_measure.h>
     1.7 +#include<lemon/random.h>
     1.8  
     1.9  using namespace lemon;
    1.10  
    1.11 @@ -100,7 +101,7 @@
    1.12    
    1.13    std::vector<Node> v;
    1.14    for(int i=0;i<N;i++) v.push_back(g.addNode());
    1.15 -  for(int i=0;i<M;i++) g.addEdge(v[int(drand48()*N)],v[int(drand48()*N)]);
    1.16 +  for(int i=0;i<M;i++) g.addEdge(v[rnd[N]],v[rnd[N]]);
    1.17  
    1.18  //   {
    1.19  //     Edge e;
     2.1 --- a/benchmark/radix_sort-bench.cc	Fri Oct 13 15:10:50 2006 +0000
     2.2 +++ b/benchmark/radix_sort-bench.cc	Sat Oct 14 15:26:05 2006 +0000
     2.3 @@ -26,6 +26,8 @@
     2.4  
     2.5  #include <lemon/radix_sort.h>
     2.6  
     2.7 +#include <lemon/random.h>
     2.8 +
     2.9  #include <vector>
    2.10  #include <algorithm>
    2.11  #include <cmath>
    2.12 @@ -37,7 +39,7 @@
    2.13    int n = 10000000;
    2.14    vector<int> data(n);
    2.15    for (int i = 0; i < n; ++i) {
    2.16 -    data[i] = (int)(1000 * (rand() / (RAND_MAX + 1.0))) - 500;
    2.17 +    data[i] = rnd[1000] - 500;
    2.18    }
    2.19    radixSort(data.begin(), data.end());
    2.20  }
    2.21 @@ -47,7 +49,7 @@
    2.22    int n = 10000000;
    2.23    vector<int> data(n);
    2.24    for (int i = 0; i < n; ++i) {
    2.25 -    data[i] = (int)(1000 * (rand() / (RAND_MAX + 1.0))) - 500;
    2.26 +    data[i] = rnd[1000] - 500;
    2.27    }
    2.28    counterSort(data.begin(), data.end());
    2.29  }
    2.30 @@ -56,7 +58,7 @@
    2.31    int n = 10000000;
    2.32    vector<int> data(n);
    2.33    for (int i = 0; i < n; ++i) {
    2.34 -    data[i] = (int)(1000 * (rand() / (RAND_MAX + 1.0))) - 500;
    2.35 +    data[i] = rnd[1000] - 500;
    2.36    }
    2.37    sort(data.begin(), data.end());
    2.38  }
    2.39 @@ -65,7 +67,7 @@
    2.40    int n = 10000000;
    2.41    vector<int> data(n);
    2.42    for (int i = 0; i < n; ++i) {
    2.43 -    data[i] = (int)(1000 * (rand() / (RAND_MAX + 1.0))) - 500;
    2.44 +    data[i] = rnd[1000] - 500;
    2.45    }
    2.46    stable_sort(data.begin(), data.end());
    2.47  }
     3.1 --- a/benchmark/swap_bipartite_bench.cc	Fri Oct 13 15:10:50 2006 +0000
     3.2 +++ b/benchmark/swap_bipartite_bench.cc	Sat Oct 14 15:26:05 2006 +0000
     3.3 @@ -12,6 +12,7 @@
     3.4  #include <lemon/graph_to_eps.h>
     3.5  
     3.6  #include <lemon/time_measure.h>
     3.7 +#include <lemon/random.h>
     3.8  
     3.9  using namespace std;
    3.10  using namespace lemon;
    3.11 @@ -20,17 +21,6 @@
    3.12  typedef ListBpUGraph LGraph;
    3.13  BPUGRAPH_TYPEDEFS(Graph);
    3.14  
    3.15 -int _urandom_init() {
    3.16 -  int seed = time(0);
    3.17 -  srand(seed);
    3.18 -  return seed;
    3.19 -}
    3.20 -
    3.21 -int urandom(int n) {
    3.22 -  static int seed = _urandom_init();
    3.23 -  ignore_unused_variable_warning(seed);
    3.24 -  return (int)(rand() / (1.0 + RAND_MAX) * n);
    3.25 -}
    3.26  
    3.27  int main() {
    3.28  
    3.29 @@ -66,8 +56,8 @@
    3.30        }
    3.31        for (int i = 0; i < e; ++i) {
    3.32          int a,b;
    3.33 -	Node aNode = aNodes[a=urandom(n)];
    3.34 -        Node bNode = bNodes[b=urandom(m)];
    3.35 +	Node aNode = aNodes[a=rnd[n]];
    3.36 +        Node bNode = bNodes[b=rnd[m]];
    3.37          graph.addEdge(aNode, bNode);
    3.38  	LGraph::Node laNode = laNodes[a];
    3.39          LGraph::Node lbNode = lbNodes[b];
     4.1 --- a/demo/descriptor_map_demo.cc	Fri Oct 13 15:10:50 2006 +0000
     4.2 +++ b/demo/descriptor_map_demo.cc	Sat Oct 14 15:26:05 2006 +0000
     4.3 @@ -32,11 +32,12 @@
     4.4  #include <lemon/dim2.h>
     4.5  #include <lemon/graph_to_eps.h>
     4.6  
     4.7 +#include <lemon/random.h>
     4.8 +
     4.9  #include <iostream>
    4.10  
    4.11 -#include <cstdlib>
    4.12  #include <cmath>
    4.13 -#include <ctime>
    4.14 +
    4.15  
    4.16  using namespace lemon;
    4.17  
    4.18 @@ -75,7 +76,6 @@
    4.19  };
    4.20  
    4.21  int main() {
    4.22 -  std::srand(std::time(0));
    4.23    typedef ListGraph Graph;
    4.24    typedef Graph::Node Node;
    4.25    typedef Graph::Edge Edge;
    4.26 @@ -105,8 +105,8 @@
    4.27    // it holds reference to it. 
    4.28    const int EDGE = (int)(NODE * std::log(double(NODE)));
    4.29    for (int i = 0; i < EDGE; ++i) {
    4.30 -    int si = (int)(std::rand() / (RAND_MAX + 1.0) * NODE);
    4.31 -    int ti = (int)(std::rand() / (RAND_MAX + 1.0) * NODE);
    4.32 +    int si = rnd[NODE];
    4.33 +    int ti = rnd[NODE];
    4.34        
    4.35      graph.addEdge(nodeInv[si], nodeInv[ti]);
    4.36    }
     5.1 --- a/demo/simann_maxcut_demo.cc	Fri Oct 13 15:10:50 2006 +0000
     5.2 +++ b/demo/simann_maxcut_demo.cc	Sat Oct 14 15:26:05 2006 +0000
     5.3 @@ -54,7 +54,7 @@
     5.4      Entity(Graph& _g, Graph::EdgeMap<int>& _w) : g(_g), w(_w), a(_g) {}
     5.5      double mutate() {
     5.6        static const int node_num = countNodes(g);
     5.7 -      int i = 1 + (int) (node_num * (rand() / (RAND_MAX + 1.0)));
     5.8 +      int i = 1 + rnd[node_num];
     5.9        NodeIt n(g);
    5.10        int j = 1;
    5.11        while (j < i) {
    5.12 @@ -91,7 +91,7 @@
    5.13        for (NodeIt n(g); n != INVALID; ++n)
    5.14          a[n] = false;
    5.15        for (NodeIt n(g); n != INVALID; ++n)
    5.16 -        if (rand() < 0.5) a[n] = true;
    5.17 +        if (rnd.boolean(0.5)) a[n] = true;
    5.18        sum = 0;
    5.19        for (EdgeIt e(g); e != INVALID; ++e)
    5.20          if (a[g.source(e)] != a[g.target(e)])
     6.1 --- a/demo/topology_demo.cc	Fri Oct 13 15:10:50 2006 +0000
     6.2 +++ b/demo/topology_demo.cc	Sat Oct 14 15:26:05 2006 +0000
     6.3 @@ -184,8 +184,6 @@
     6.4  } 
     6.5  
     6.6  int main() {
     6.7 -  srand(time(0));
     6.8 -
     6.9    drawConnectedComponents();
    6.10    drawStronglyConnectedComponents();
    6.11    drawNodeBiconnectedComponents();
     7.1 --- a/lemon/hypercube_graph.h	Fri Oct 13 15:10:50 2006 +0000
     7.2 +++ b/lemon/hypercube_graph.h	Sat Oct 14 15:26:05 2006 +0000
     7.3 @@ -260,8 +260,8 @@
     7.4      /// HyperCubeGraph graph(DIM);
     7.5      /// dim2::Point<double> base[DIM];
     7.6      /// for (int k = 0; k < DIM; ++k) {
     7.7 -    ///   base[k].x = random.getReal();
     7.8 -    ///   base[k].y = random.getReal();
     7.9 +    ///   base[k].x = rnd();
    7.10 +    ///   base[k].y = rnd();
    7.11      /// } 
    7.12      /// HyperCubeGraph::HyperMap<dim2::Point<double> > 
    7.13      ///   pos(graph, base, base + DIM, dim2::Point<double>(0.0, 0.0));
     8.1 --- a/lemon/random.h	Fri Oct 13 15:10:50 2006 +0000
     8.2 +++ b/lemon/random.h	Sat Oct 14 15:26:05 2006 +0000
     8.3 @@ -20,6 +20,8 @@
     8.4  #define LEMON_RANDOM_H
     8.5  
     8.6  #include <algorithm>
     8.7 +#include <iterator>
     8.8 +#include <vector>
     8.9  
    8.10  #include <ctime>
    8.11  #include <cmath>
    8.12 @@ -33,46 +35,460 @@
    8.13  
    8.14  namespace lemon {
    8.15  
    8.16 -#if WORD_BIT == 32
    8.17 +  namespace _random_bits {
    8.18 +    
    8.19 +    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
    8.20 +    struct RandomTraits {};
    8.21 +
    8.22 +    template <typename _Word>
    8.23 +    struct RandomTraits<_Word, 32> {
    8.24 +
    8.25 +      typedef _Word Word;
    8.26 +      static const int bits = 32;
    8.27 +
    8.28 +      static const int length = 624;
    8.29 +      static const int shift = 397;
    8.30 +      
    8.31 +      static const Word mul = 0x6c078965u;
    8.32 +      static const Word arrayInit = 0x012BD6AAu;
    8.33 +      static const Word arrayMul1 = 0x0019660Du;
    8.34 +      static const Word arrayMul2 = 0x5D588B65u;
    8.35 +
    8.36 +      static const Word mask = 0x9908B0DFu;
    8.37 +      static const Word loMask = (1u << 31) - 1;
    8.38 +      static const Word hiMask = ~loMask;
    8.39 +
    8.40 +
    8.41 +      static Word tempering(Word rnd) {
    8.42 +        rnd ^= (rnd >> 11);
    8.43 +        rnd ^= (rnd << 7) & 0x9D2C5680u;
    8.44 +        rnd ^= (rnd << 15) & 0xEFC60000u;
    8.45 +        rnd ^= (rnd >> 18);
    8.46 +        return rnd;
    8.47 +      }
    8.48 +
    8.49 +    };
    8.50 +
    8.51 +    template <typename _Word>
    8.52 +    struct RandomTraits<_Word, 64> {
    8.53 +
    8.54 +      typedef _Word Word;
    8.55 +      static const int bits = 64;
    8.56 +
    8.57 +      static const int length = 312;
    8.58 +      static const int shift = 156;
    8.59 +
    8.60 +      static const Word mul = (Word)0x5851F42Du << 32 | (Word)0x4C957F2Du;
    8.61 +      static const Word arrayInit = (Word)0x00000000u << 32 |(Word)0x012BD6AAu;
    8.62 +      static const Word arrayMul1 = (Word)0x369DEA0Fu << 32 |(Word)0x31A53F85u;
    8.63 +      static const Word arrayMul2 = (Word)0x27BB2EE6u << 32 |(Word)0x87B0B0FDu;
    8.64 +
    8.65 +      static const Word mask = (Word)0xB5026F5Au << 32 | (Word)0xA96619E9u;
    8.66 +      static const Word loMask = ((Word)1u << 31) - 1;
    8.67 +      static const Word hiMask = ~loMask;
    8.68 +
    8.69 +      static Word tempering(Word rnd) {
    8.70 +        rnd ^= (rnd >> 29) & ((Word)0x55555555u << 32 | (Word)0x55555555u);
    8.71 +        rnd ^= (rnd << 17) & ((Word)0x71D67FFFu << 32 | (Word)0xEDA60000u);
    8.72 +        rnd ^= (rnd << 37) & ((Word)0xFFF7EEE0u << 32 | (Word)0x00000000u);
    8.73 +        rnd ^= (rnd >> 43);
    8.74 +        return rnd;
    8.75 +      }
    8.76 +
    8.77 +    };
    8.78 +
    8.79 +    template <typename _Word>
    8.80 +    class RandomCore {
    8.81 +    public:
    8.82 +
    8.83 +      typedef _Word Word;
    8.84 +
    8.85 +    private:
    8.86 +
    8.87 +      static const int bits = RandomTraits<Word>::bits;
    8.88 +
    8.89 +      static const int length = RandomTraits<Word>::length;
    8.90 +      static const int shift = RandomTraits<Word>::shift;
    8.91 +
    8.92 +    public:
    8.93 +
    8.94 +      void initState() {
    8.95 +        static const Word seedArray[4] = {
    8.96 +          0x12345u, 0x23456u, 0x34567u, 0x45678u
    8.97 +        };
    8.98 +    
    8.99 +        initState(seedArray, seedArray + 4);
   8.100 +      }
   8.101 +
   8.102 +      void initState(Word seed) {
   8.103 +
   8.104 +        static const Word mul = RandomTraits<Word>::mul;
   8.105 +
   8.106 +        current = state; 
   8.107 +
   8.108 +        Word *curr = state + length - 1;
   8.109 +        curr[0] = seed; --curr;
   8.110 +        for (int i = 1; i < length; ++i) {
   8.111 +          curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
   8.112 +          --curr;
   8.113 +        }
   8.114 +      }
   8.115 +
   8.116 +      template <typename Iterator>
   8.117 +      void initState(Iterator begin, Iterator end) {
   8.118 +
   8.119 +        static const Word init = RandomTraits<Word>::arrayInit;
   8.120 +        static const Word mul1 = RandomTraits<Word>::arrayMul1;
   8.121 +        static const Word mul2 = RandomTraits<Word>::arrayMul2;
   8.122 +
   8.123 +
   8.124 +        Word *curr = state + length - 1; --curr;
   8.125 +        Iterator it = begin; int cnt = 0;
   8.126 +        int num;
   8.127 +
   8.128 +        initState(init);
   8.129 +
   8.130 +        num = length > end - begin ? length : end - begin;
   8.131 +        while (num--) {
   8.132 +          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 
   8.133 +            + *it + cnt;
   8.134 +          ++it; ++cnt;
   8.135 +          if (it == end) {
   8.136 +            it = begin; cnt = 0;
   8.137 +          }
   8.138 +          if (curr == state) {
   8.139 +            curr = state + length - 1; curr[0] = state[0];
   8.140 +          }
   8.141 +          --curr;
   8.142 +        }
   8.143 +
   8.144 +        num = length - 1; cnt = length - (curr - state) - 1;
   8.145 +        while (num--) {
   8.146 +          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
   8.147 +            - cnt;
   8.148 +          --curr; ++cnt;
   8.149 +          if (curr == state) {
   8.150 +            curr = state + length - 1; curr[0] = state[0]; --curr;
   8.151 +            cnt = 1;
   8.152 +          }
   8.153 +        }
   8.154 +        
   8.155 +        state[length - 1] = (Word)1 << (bits - 1);
   8.156 +      }
   8.157 +      
   8.158 +      void copyState(const RandomCore& other) {
   8.159 +        std::copy(other.state, other.state + length, state);
   8.160 +        current = state + (other.current - other.state);
   8.161 +      }
   8.162 +
   8.163 +      Word operator()() {
   8.164 +        if (current == state) fillState();
   8.165 +        --current;
   8.166 +        Word rnd = *current;
   8.167 +        return RandomTraits<Word>::tempering(rnd);
   8.168 +      }
   8.169 +
   8.170 +    private:
   8.171 +
   8.172 +  
   8.173 +      void fillState() {
   8.174 +        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
   8.175 +        static const Word loMask = RandomTraits<Word>::loMask;
   8.176 +        static const Word hiMask = RandomTraits<Word>::hiMask;
   8.177 +
   8.178 +        current = state + length; 
   8.179 +
   8.180 +        register Word *curr = state + length - 1;
   8.181 +        register long num;
   8.182 +      
   8.183 +        num = length - shift;
   8.184 +        while (num--) {
   8.185 +          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   8.186 +            curr[- shift] ^ mask[curr[-1] & 1ul];
   8.187 +          --curr;
   8.188 +        }
   8.189 +        num = shift - 1;
   8.190 +        while (num--) {
   8.191 +          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   8.192 +            curr[length - shift] ^ mask[curr[-1] & 1ul];
   8.193 +          --curr;
   8.194 +        }
   8.195 +        curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   8.196 +          curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   8.197 +
   8.198 +      }
   8.199 +
   8.200 +  
   8.201 +      Word *current;
   8.202 +      Word state[length];
   8.203 +      
   8.204 +    };
   8.205 +
   8.206 +
   8.207 +    template <typename Result, 
   8.208 +              int shift = (std::numeric_limits<Result>::digits + 1) / 2>
   8.209 +    struct Masker {
   8.210 +      static Result mask(const Result& result) {
   8.211 +        return Masker<Result, (shift + 1) / 2>::
   8.212 +          mask((Result)(result | (result >> shift)));
   8.213 +      }
   8.214 +    };
   8.215 +    
   8.216 +    template <typename Result>
   8.217 +    struct Masker<Result, 1> {
   8.218 +      static Result mask(const Result& result) {
   8.219 +        return (Result)(result | (result >> 1));
   8.220 +      }
   8.221 +    };
   8.222 +
   8.223 +    template <typename Result, typename Word, 
   8.224 +              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   8.225 +              bool last = rest <= std::numeric_limits<Word>::digits>
   8.226 +    struct IntConversion {
   8.227 +      static const int bits = std::numeric_limits<Word>::digits;
   8.228 +    
   8.229 +      static Result convert(RandomCore<Word>& rnd) {
   8.230 +        return (Result)(rnd() >> (bits - rest)) << shift;
   8.231 +      }
   8.232 +      
   8.233 +    }; 
   8.234 +
   8.235 +    template <typename Result, typename Word, int rest, int shift> 
   8.236 +    struct IntConversion<Result, Word, rest, shift, false> {
   8.237 +      static const int bits = std::numeric_limits<Word>::digits;
   8.238 +
   8.239 +      static Result convert(RandomCore<Word>& rnd) {
   8.240 +        return ((Result)rnd() << shift) | 
   8.241 +          IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
   8.242 +      }
   8.243 +    };
   8.244 +
   8.245 +
   8.246 +    template <typename Result, typename Word,
   8.247 +              bool one_word = std::numeric_limits<Word>::digits < 
   8.248 +                              std::numeric_limits<Result>::digits>
   8.249 +    struct Mapping {
   8.250 +      static Result map(RandomCore<Word>& rnd, const Result& bound) {
   8.251 +        Word max = (Word)(bound - 1);
   8.252 +        Result mask = Masker<Result>::mask(bound - 1);
   8.253 +        Result num;
   8.254 +        do {
   8.255 +          num = IntConversion<Result, Word>::convert(rnd) & mask; 
   8.256 +        } while (num > max);
   8.257 +        return num;
   8.258 +      }
   8.259 +    };
   8.260 +
   8.261 +    template <typename Result, typename Word>
   8.262 +    struct Mapping<Result, Word, false> {
   8.263 +      static Result map(RandomCore<Word>& rnd, const Result& bound) {
   8.264 +        Word max = (Word)(bound - 1);
   8.265 +        Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
   8.266 +          ::mask(max);
   8.267 +        Word num;
   8.268 +        do {
   8.269 +          num = rnd() & mask;
   8.270 +        } while (num > max);
   8.271 +        return num;
   8.272 +      }
   8.273 +    };
   8.274 +
   8.275 +    template <typename Result, int exp, bool pos = (exp >= 0)>
   8.276 +    struct ShiftMultiplier {
   8.277 +      static const Result multiplier() {
   8.278 +        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   8.279 +        res *= res;
   8.280 +        if ((exp & 1) == 1) res *= (Result)2.0;
   8.281 +        return res; 
   8.282 +      }
   8.283 +    };
   8.284 +
   8.285 +    template <typename Result, int exp>
   8.286 +    struct ShiftMultiplier<Result, exp, false> {
   8.287 +      static const Result multiplier() {
   8.288 +        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   8.289 +        res *= res;
   8.290 +        if ((exp & 1) == 1) res *= (Result)0.5;
   8.291 +        return res; 
   8.292 +      }
   8.293 +    };
   8.294 +
   8.295 +    template <typename Result>
   8.296 +    struct ShiftMultiplier<Result, 0, true> {
   8.297 +      static const Result multiplier() {
   8.298 +        return (Result)1.0; 
   8.299 +      }
   8.300 +    };
   8.301 +
   8.302 +    template <typename Result>
   8.303 +    struct ShiftMultiplier<Result, -20, true> {
   8.304 +      static const Result multiplier() {
   8.305 +        return (Result)(1.0/1048576.0); 
   8.306 +      }
   8.307 +    };
   8.308 +    
   8.309 +    template <typename Result>
   8.310 +    struct ShiftMultiplier<Result, -32, true> {
   8.311 +      static const Result multiplier() {
   8.312 +        return (Result)(1.0/424967296.0); 
   8.313 +      }
   8.314 +    };
   8.315 +
   8.316 +    template <typename Result>
   8.317 +    struct ShiftMultiplier<Result, -53, true> {
   8.318 +      static const Result multiplier() {
   8.319 +        return (Result)(1.0/9007199254740992.0); 
   8.320 +      }
   8.321 +    };
   8.322 +
   8.323 +    template <typename Result>
   8.324 +    struct ShiftMultiplier<Result, -64, true> {
   8.325 +      static const Result multiplier() {
   8.326 +        return (Result)(1.0/18446744073709551616.0); 
   8.327 +      }
   8.328 +    };
   8.329 +
   8.330 +    template <typename Result, int exp>
   8.331 +    struct Shifting {
   8.332 +      static Result shift(const Result& result) {
   8.333 +        return result * ShiftMultiplier<Result, exp>::multiplier();
   8.334 +      }
   8.335 +    };
   8.336 +
   8.337 +    template <typename Result, typename Word,
   8.338 +              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   8.339 +              bool last = rest <= std::numeric_limits<Word>::digits>
   8.340 +    struct RealConversion{ 
   8.341 +      static const int bits = std::numeric_limits<Word>::digits;
   8.342 +
   8.343 +      static Result convert(RandomCore<Word>& rnd) {
   8.344 +        return Shifting<Result, - shift - rest>::
   8.345 +          shift((Result)(rnd() >> (bits - rest)));
   8.346 +      }
   8.347 +    };
   8.348 +
   8.349 +    template <typename Result, typename Word, int rest, int shift>
   8.350 +    struct RealConversion<Result, Word, rest, shift, false> { 
   8.351 +      static const int bits = std::numeric_limits<Word>::digits;
   8.352 +
   8.353 +      static Result convert(RandomCore<Word>& rnd) {
   8.354 +        return Shifting<Result, - shift - bits>::shift((Result)rnd()) +
   8.355 +          RealConversion<Result, Word, rest-bits, shift + bits>::convert(rnd);
   8.356 +      }
   8.357 +    };
   8.358 +
   8.359 +    template <typename Result, typename Word>
   8.360 +    struct Initializer {
   8.361 +
   8.362 +      template <typename Iterator>
   8.363 +      static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
   8.364 +        std::vector<Word> ws;
   8.365 +        for (Iterator it = begin; it != end; ++it) {
   8.366 +          ws.push_back((Word)*it);
   8.367 +        }
   8.368 +        rnd.initState(ws.begin(), ws.end());
   8.369 +      }
   8.370 +
   8.371 +      template <typename Iterator>
   8.372 +      static void init(RandomCore<Word>& rnd, Result seed) {
   8.373 +        rnd.initState(seed);
   8.374 +      }
   8.375 +    };
   8.376 +
   8.377 +    template <typename Word>
   8.378 +    struct BoolConversion {
   8.379 +      static bool convert(RandomCore<Word>& rnd) {
   8.380 +        return (rnd() & 1) == 1;
   8.381 +      }
   8.382 +    };
   8.383 +
   8.384 +  }
   8.385  
   8.386    /// \ingroup misc
   8.387    ///
   8.388    /// \brief Mersenne Twister random number generator
   8.389    ///
   8.390    /// The Mersenne Twister is a twisted generalized feedback
   8.391 -  /// shift-register generator of Matsumoto and Nishimura. The period of
   8.392 -  /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
   8.393 -  /// 623 dimensions. The time performance of this generator is comparable
   8.394 -  /// to the commonly used generators.
   8.395 +  /// shift-register generator of Matsumoto and Nishimura. The period
   8.396 +  /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
   8.397 +  /// equi-distributed in 623 dimensions for 32-bit numbers. The time
   8.398 +  /// performance of this generator is comparable to the commonly used
   8.399 +  /// generators.
   8.400 +  ///
   8.401 +  /// This implementation is specialized for both 32-bit and 64-bit
   8.402 +  /// architectures. The generators differ sligthly in the
   8.403 +  /// initialization and generation phase so they produce two
   8.404 +  /// completly different sequences.
   8.405 +  ///
   8.406 +  /// The generator gives back random numbers of serveral types. To
   8.407 +  /// get a random number from a range of a floating point type you
   8.408 +  /// can use one form of the \c operator(). If you want to get random
   8.409 +  /// number from a the {0, 1, ..., n-1} integer range use the \c
   8.410 +  /// operator[]. And to get random number from the whole range of an
   8.411 +  /// integer type you can use the \c integer() or \c uinteger()
   8.412 +  /// functions. After all you can get random bool with equal chance
   8.413 +  /// or given probability with the \c boolean() member function.
   8.414 +  ///
   8.415 +  ///\code
   8.416 +  /// int a = rnd[100000];                  // 0..99999
   8.417 +  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
   8.418 +  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
   8.419 +  /// double d = rnd();                     // [0.0, 1.0)
   8.420 +  /// double e = rnd(100.0);                // [0.0, 100.0)
   8.421 +  /// double f = rnd(1.0, 2.0);             // [1.0, 2.0)
   8.422 +  /// bool g = rnd.boolean();               // P(g = true) = 0.5
   8.423 +  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
   8.424 +  ///\endcode
   8.425 +  ///
   8.426 +  /// The lemon provides a global instance of the random number generator
   8.427 +  /// which name is \c rnd. Usually it is a good programming convenience
   8.428 +  /// to use this global generator to get random numbers.
   8.429    ///
   8.430    /// \author Balazs Dezso
   8.431    class Random {
   8.432 +  private:
   8.433  
   8.434 -    static const int length = 624;
   8.435 -    static const int shift = 397;
   8.436 +#if WORD_BIT == 32
   8.437 +    typedef unsigned int Word;
   8.438 +#elif WORD_BIT == 64
   8.439 +    typedef unsigned long Word;
   8.440 +#endif
   8.441 +    
   8.442 +    _random_bits::RandomCore<Word> core;
   8.443  
   8.444    public:
   8.445  
   8.446 -    static const unsigned long min = 0ul;
   8.447 -    static const unsigned long max = ~0ul;
   8.448 -  
   8.449      /// \brief Constructor
   8.450      ///
   8.451 -    /// Constructor with time dependent seeding.
   8.452 -    Random() { initState(std::time(0)); }
   8.453 +    /// Constructor with constant seeding.
   8.454 +    Random() { core.initState(); }
   8.455  
   8.456      /// \brief Constructor
   8.457      ///
   8.458 -    /// Constructor
   8.459 -    Random(unsigned long seed) { initState(seed); }
   8.460 +    /// Constructor with seed. The current number type will be converted
   8.461 +    /// to the architecture word type.
   8.462 +    template <typename Number>
   8.463 +    Random(Number seed) { 
   8.464 +      _random_bits::Initializer<Number, Word>::init(core, seed);
   8.465 +    }
   8.466 +
   8.467 +    /// \brief Constructor
   8.468 +    ///
   8.469 +    /// Constructor with array seeding. The given range should contain
   8.470 +    /// any number type and the numbers will be converted to the
   8.471 +    /// architecture word type.
   8.472 +    template <typename Iterator>
   8.473 +    Random(Iterator begin, Iterator end) { 
   8.474 +      typedef typename std::iterator_traits<Iterator>::value_type Number;
   8.475 +      _random_bits::Initializer<Number, Word>::init(core, begin, end);
   8.476 +    }
   8.477  
   8.478      /// \brief Copy constructor
   8.479      ///
   8.480      /// Copy constructor. The generated sequence will be identical to
   8.481      /// the other sequence.
   8.482 -    Random(const Random& other) { 
   8.483 -      std::copy(other.state, other.state + length, state);
   8.484 -      current = state + (other.current - other.state);
   8.485 +    Random(const Random& other) {
   8.486 +      core.copyState(other.core);
   8.487      }
   8.488  
   8.489      /// \brief Assign operator
   8.490 @@ -81,429 +497,94 @@
   8.491      /// the other sequence.
   8.492      Random& operator=(const Random& other) {
   8.493        if (&other != this) {
   8.494 -        std::copy(other.state, other.state + length, state);
   8.495 -        current = state + (other.current - other.state);
   8.496 +        core.copyState(other.core);
   8.497        }
   8.498        return *this;
   8.499      }
   8.500  
   8.501 -    /// \brief Returns an unsigned random number
   8.502 +    /// \brief Returns a random real number
   8.503      ///
   8.504 -    /// It returns an unsigned integer random number from the range 
   8.505 -    /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$.
   8.506 -    unsigned long getUnsigned() {
   8.507 -      if (current == state) fillState();
   8.508 -      --current;
   8.509 -      unsigned long rnd = *current;
   8.510 -      rnd ^= (rnd >> 11);
   8.511 -      rnd ^= (rnd << 7) & 0x9D2C5680ul;
   8.512 -      rnd ^= (rnd << 15) & 0xEFC60000ul;
   8.513 -      rnd ^= (rnd >> 18);
   8.514 -      return rnd;
   8.515 +    /// It returns a random real number from the range [0, 1). The default
   8.516 +    /// result type of this function is double.
   8.517 +    template <typename Number>
   8.518 +    Number operator()() {
   8.519 +      return _random_bits::RealConversion<Number, Word>::convert(core);
   8.520      }
   8.521  
   8.522 -    /// \brief Returns a random number
   8.523 +    double operator()() {
   8.524 +      return operator()<double>();
   8.525 +    }
   8.526 +
   8.527 +    /// \brief Returns a random real number
   8.528      ///
   8.529 -    /// It returns an integer random number from the range 
   8.530 -    /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$.
   8.531 -    long getInt() {
   8.532 -      return (long)getUnsigned();
   8.533 +    /// It returns a random real number from the range [0, b).
   8.534 +    template <typename Number>
   8.535 +    Number operator()(Number b) { 
   8.536 +      return operator()<Number>() * b; 
   8.537 +    }
   8.538 +
   8.539 +    /// \brief Returns a random real number
   8.540 +    ///
   8.541 +    /// It returns a random real number from the range [a, b).
   8.542 +    template <typename Number>
   8.543 +    Number operator()(Number a, Number b) { 
   8.544 +      return operator()<Number>() * (b - a) + a; 
   8.545 +    }
   8.546 +
   8.547 +    /// \brief Returns a random integer from a range
   8.548 +    ///
   8.549 +    /// It returns a random integer from the range {0, 1, ..., bound - 1}.
   8.550 +    template <typename Number>
   8.551 +    Number operator[](const Number& bound) {
   8.552 +      return _random_bits::Mapping<Number, Word>::map(core, bound);
   8.553 +    }
   8.554 +
   8.555 +    /// \brief Returns a random non-negative integer
   8.556 +    ///
   8.557 +    /// It returns a random non-negative integer uniformly from the
   8.558 +    /// whole range of the current \c Number type.  The default result
   8.559 +    /// type of this function is unsigned int.
   8.560 +    template <typename Number>
   8.561 +    Number uinteger() {
   8.562 +      return _random_bits::IntConversion<Number, Word>::convert(core);
   8.563 +    }
   8.564 +
   8.565 +    unsigned int uinteger() {
   8.566 +      return uinteger<unsigned int>();
   8.567 +    }
   8.568 +
   8.569 +    /// \brief Returns a random integer
   8.570 +    ///
   8.571 +    /// It returns a random integer uniformly from the whole range of
   8.572 +    /// the current \c Number type. The default result type of this
   8.573 +    /// function is int.
   8.574 +    template <typename Number>
   8.575 +    Number integer() {
   8.576 +      static const int nb = std::numeric_limits<Number>::digits + 
   8.577 +        (std::numeric_limits<Number>::is_signed ? 1 : 0);
   8.578 +      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
   8.579 +    }
   8.580 +
   8.581 +    int integer() {
   8.582 +      return integer<int>();
   8.583      }
   8.584      
   8.585 -    /// \brief Returns an unsigned integer random number
   8.586 +    /// \brief Returns a random bool
   8.587      ///
   8.588 -    /// It returns an unsigned integer random number from the range 
   8.589 -    /// \f$ \{0, 1, \dots, 2^{31} - 1\} \f$.
   8.590 -    long getNatural() {
   8.591 -      return (long)(getUnsigned() >> 1);
   8.592 +    /// It returns a random bool
   8.593 +    bool boolean() {
   8.594 +      return _random_bits::BoolConversion<Word>::convert(core);
   8.595      }
   8.596  
   8.597      /// \brief Returns a random bool
   8.598      ///
   8.599 -    /// It returns a random bool.
   8.600 -    bool getBool() {
   8.601 -      return (bool)(getUnsigned() & 1);
   8.602 +    /// It returns a random bool with given probability of true result
   8.603 +    bool boolean(double p) {
   8.604 +      return operator()() < p;
   8.605      }
   8.606 -
   8.607 -    /// \brief Returns a real random number
   8.608 -    ///
   8.609 -    /// It returns a real random number from the range 
   8.610 -    /// \f$ [0, 1) \f$. The double will have 32 significant bits.
   8.611 -    double getReal() {
   8.612 -      return std::ldexp((double)getUnsigned(), -32);
   8.613 -    }
   8.614 -
   8.615 -    /// \brief Returns a real random number
   8.616 -    ///
   8.617 -    /// It returns a real random number from the range 
   8.618 -    /// \f$ [0, 1) \f$. The double will have 53 significant bits,
   8.619 -    /// but it is slower than the \c getReal().
   8.620 -    double getPrecReal() {
   8.621 -      return std::ldexp((double)(getUnsigned() >> 5), -27) + 
   8.622 -        std::ldexp((double)(getUnsigned() >> 6), -53);
   8.623 -    }
   8.624 -
   8.625 -    /// \brief Returns an unsigned random number from a given range
   8.626 -    ///
   8.627 -    /// It returns an unsigned integer random number from the range 
   8.628 -    /// \f$ \{0, 1, \dots, n - 1\} \f$.
   8.629 -    unsigned long getUnsigned(unsigned long n) {
   8.630 -      unsigned long mask = n - 1, rnd;
   8.631 -      mask |= mask >> 1;
   8.632 -      mask |= mask >> 2;
   8.633 -      mask |= mask >> 4;
   8.634 -      mask |= mask >> 8;
   8.635 -      mask |= mask >> 16;
   8.636 -      do rnd = getUnsigned() & mask; while (rnd >= n);
   8.637 -      return rnd;
   8.638 -    }
   8.639 -
   8.640 -    /// \brief Returns a random number from a given range
   8.641 -    ///
   8.642 -    /// It returns an unsigned integer random number from the range 
   8.643 -    /// \f$ \{0, 1, \dots, n - 1\} \f$.
   8.644 -    long getInt(long n) {
   8.645 -      long mask = n - 1, rnd;
   8.646 -      mask |= mask >> 1;
   8.647 -      mask |= mask >> 2;
   8.648 -      mask |= mask >> 4;
   8.649 -      mask |= mask >> 8;
   8.650 -      mask |= mask >> 16;
   8.651 -      do rnd = getUnsigned() & mask; while (rnd >= n);
   8.652 -      return rnd;
   8.653 -    }
   8.654 -
   8.655 -  private:
   8.656 -
   8.657 -    void initState(unsigned long seed) {
   8.658 -      static const unsigned long mul = 0x6c078965ul;
   8.659 -
   8.660 -      current = state; 
   8.661 -
   8.662 -      unsigned long *curr = state + length - 1;
   8.663 -      curr[0] = seed; --curr;
   8.664 -      for (int i = 1; i < length; ++i) {
   8.665 -        curr[0] = (mul * ( curr[1] ^ (curr[1] >> 30) ) + i);
   8.666 -        --curr;
   8.667 -      }
   8.668 -    }
   8.669 -  
   8.670 -    void fillState() {
   8.671 -      static const unsigned long mask[2] = { 0x0ul, 0x9908B0DFul };
   8.672 -      static const unsigned long loMask = (1ul << 31) - 1;
   8.673 -      static const unsigned long hiMask = ~loMask;
   8.674 -
   8.675 -      current = state + length; 
   8.676 -
   8.677 -      register unsigned long *curr = state + length - 1;
   8.678 -      register long num;
   8.679 -      
   8.680 -      num = length - shift;
   8.681 -      while (num--) {
   8.682 -        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   8.683 -          curr[- shift] ^ mask[curr[-1] & 1ul];
   8.684 -        --curr;
   8.685 -      }
   8.686 -      num = shift - 1;
   8.687 -      while (num--) {
   8.688 -        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   8.689 -          curr[length - shift] ^ mask[curr[-1] & 1ul];
   8.690 -        --curr;
   8.691 -      }
   8.692 -      curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   8.693 -        curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   8.694 -
   8.695 -    }
   8.696 -  
   8.697 -    unsigned long *current;
   8.698 -    unsigned long state[length];
   8.699      
   8.700    };
   8.701  
   8.702 -#elif WORD_BIT == 64
   8.703 -
   8.704 -  /// \ingroup misc
   8.705 -  ///
   8.706 -  /// \brief Mersenne Twister random number generator
   8.707 -  ///
   8.708 -  /// The Mersenne Twister is a twisted generalized feedback
   8.709 -  /// shift-register generator of Matsumoto and Nishimura. The period of
   8.710 -  /// this generator is \f$ 2^{19937} - 1 \f$ and it is equi-distributed in
   8.711 -  /// 311 dimensions. The time performance of this generator is comparable
   8.712 -  /// to the commonly used generators.
   8.713 -  class Random {
   8.714 -
   8.715 -    static const int length = 312;
   8.716 -    static const int shift = 156;
   8.717 -
   8.718 -  public:
   8.719 -
   8.720 -    static const unsigned long min = 0ul;
   8.721 -    static const unsigned long max = ~0ul;
   8.722 -  
   8.723 -    /// \brief Constructor
   8.724 -    ///
   8.725 -    /// Constructor with time dependent seeding.
   8.726 -    Random() { initState(std::time(0)); }
   8.727 -
   8.728 -    /// \brief Constructor
   8.729 -    ///
   8.730 -    /// Constructor
   8.731 -    Random(unsigned long seed) { initState(seed); }
   8.732 -
   8.733 -    /// \brief Copy constructor
   8.734 -    ///
   8.735 -    /// Copy constructor. The generated sequence will be identical to
   8.736 -    /// the other sequence.
   8.737 -    Random(const Random& other) { 
   8.738 -      std::copy(other.state, other.state + length, state);
   8.739 -      current = state + (other.current - other.state);
   8.740 -    }
   8.741 -
   8.742 -    /// \brief Assign operator
   8.743 -    ///
   8.744 -    /// Assign operator. The generated sequence will be identical to
   8.745 -    /// the other sequence.
   8.746 -    Random& operator=(const Random& other) {
   8.747 -      if (&other != this) {
   8.748 -        std::copy(other.state, other.state + length, state);
   8.749 -        current = state + (other.current - other.state);
   8.750 -      }
   8.751 -      return *this;
   8.752 -    }
   8.753 -
   8.754 -    /// \brief Returns an unsigned random number
   8.755 -    ///
   8.756 -    /// It returns an unsigned integer random number from the range 
   8.757 -    /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$.
   8.758 -    unsigned long getUnsigned() {
   8.759 -      if (current == state) fillState();
   8.760 -      --current;
   8.761 -      unsigned long rnd = *current;
   8.762 -      rnd ^= (rnd >> 29) & 0x5555555555555555ul;
   8.763 -      rnd ^= (rnd << 17) & 0x71D67FFFEDA60000ul;
   8.764 -      rnd ^= (rnd << 37) & 0xFFF7EEE000000000ul;
   8.765 -      rnd ^= (rnd >> 43);
   8.766 -      return rnd;
   8.767 -    }
   8.768 -
   8.769 -    /// \brief Returns a random number
   8.770 -    ///
   8.771 -    /// It returns an integer random number from the range 
   8.772 -    /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$.
   8.773 -    long getInt() {
   8.774 -      return (long)getUnsigned();
   8.775 -    }
   8.776 -    
   8.777 -    /// \brief Returns an unsigned integer random number
   8.778 -    ///
   8.779 -    /// It returns an unsigned integer random number from the range 
   8.780 -    /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$.
   8.781 -    long getNatural() {
   8.782 -      return (long)(getUnsigned() >> 1);
   8.783 -    }
   8.784 -
   8.785 -    /// \brief Returns a random bool
   8.786 -    ///
   8.787 -    /// It returns a random bool.
   8.788 -    bool getBool() {
   8.789 -      return (bool)(getUnsigned() & 1);
   8.790 -    }
   8.791 -
   8.792 -    /// \brief Returns a real random number
   8.793 -    ///
   8.794 -    /// It returns a real random number from the range 
   8.795 -    /// \f$ [0, 1) \f$.
   8.796 -    double getReal() {
   8.797 -      return std::ldexp((double)(getUnsigned() >> 11), -53);
   8.798 -    }
   8.799 -
   8.800 -    /// \brief Returns a real random number
   8.801 -    ///
   8.802 -    /// It returns a real random number from the range 
   8.803 -    /// \f$ [0, 1) \f$. This function is identical to the \c getReal().
   8.804 -    double getPrecReal() {
   8.805 -      return getReal();
   8.806 -    }
   8.807 -
   8.808 -    /// \brief Returns an unsigned random number from a given range
   8.809 -    ///
   8.810 -    /// It returns an unsigned integer random number from the range 
   8.811 -    /// \f$ \{0, 1, \dots, n - 1\} \f$.
   8.812 -    unsigned long getUnsigned(unsigned long n) {
   8.813 -      unsigned long mask = n - 1, rnd;
   8.814 -      mask |= mask >> 1;
   8.815 -      mask |= mask >> 2;
   8.816 -      mask |= mask >> 4;
   8.817 -      mask |= mask >> 8;
   8.818 -      mask |= mask >> 16;
   8.819 -      mask |= mask >> 32;
   8.820 -      do rnd = getUnsigned() & mask; while (rnd >= n);
   8.821 -      return rnd;
   8.822 -    }
   8.823 -
   8.824 -    /// \brief Returns a random number from a given range
   8.825 -    ///
   8.826 -    /// It returns an unsigned integer random number from the range 
   8.827 -    /// \f$ \{0, 1, \dots, n - 1\} \f$.
   8.828 -    long getInt(long n) {
   8.829 -      long mask = n - 1, rnd;
   8.830 -      mask |= mask >> 1;
   8.831 -      mask |= mask >> 2;
   8.832 -      mask |= mask >> 4;
   8.833 -      mask |= mask >> 8;
   8.834 -      mask |= mask >> 16;
   8.835 -      mask |= mask >> 32;
   8.836 -      do rnd = getUnsigned() & mask; while (rnd >= n);
   8.837 -      return rnd;
   8.838 -    }
   8.839 -
   8.840 -  private:
   8.841 -
   8.842 -    void initState(unsigned long seed) {
   8.843 -
   8.844 -      static const unsigned long mul = 0x5851F42D4C957F2Dul;
   8.845 -
   8.846 -      current = state; 
   8.847 -
   8.848 -      unsigned long *curr = state + length - 1;
   8.849 -      curr[0] = seed; --curr;
   8.850 -      for (int i = 1; i < length; ++i) {
   8.851 -        curr[0] = (mul * ( curr[1] ^ (curr[1] >> 62) ) + i);
   8.852 -        --curr;
   8.853 -      }
   8.854 -    }
   8.855 -  
   8.856 -    void fillState() {
   8.857 -      static const unsigned long mask[2] = { 0x0ul, 0xB5026F5AA96619E9ul };
   8.858 -      static const unsigned long loMask = (1ul << 31) - 1;
   8.859 -      static const unsigned long hiMask = ~loMask;
   8.860 -
   8.861 -      current = state + length; 
   8.862 -
   8.863 -      register unsigned long *curr = state + length - 1;
   8.864 -      register int num;
   8.865 -      
   8.866 -      num = length - shift;
   8.867 -      while (num--) {
   8.868 -        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   8.869 -          curr[- shift] ^ mask[curr[-1] & 1ul];
   8.870 -        --curr;
   8.871 -      }
   8.872 -      num = shift - 1;
   8.873 -      while (num--) {
   8.874 -        curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   8.875 -          curr[length - shift] ^ mask[curr[-1] & 1ul];
   8.876 -        --curr;
   8.877 -      }
   8.878 -      curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   8.879 -        curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   8.880 -
   8.881 -    }
   8.882 -  
   8.883 -    unsigned long *current;
   8.884 -    unsigned long state[length];
   8.885 -    
   8.886 -  };
   8.887 -
   8.888 -#else
   8.889 -
   8.890 -  /// \ingroup misc
   8.891 -  ///
   8.892 -  /// \brief Mersenne Twister random number generator
   8.893 -  ///
   8.894 -  /// The Mersenne Twister is a twisted generalized feedback
   8.895 -  /// shift-register generator of Matsumoto and Nishimura. There is
   8.896 -  /// two different implementation in the Lemon library, one for
   8.897 -  /// 32-bit processors and one for 64-bit processors. The period of
   8.898 -  /// the generated sequence is \f$ 2^{19937} - 1 \f$, the generated
   8.899 -  /// sequence of 32-bit random numbers is equi-distributed in 623
   8.900 -  /// dimensions. The time performance of this generator is comparable
   8.901 -  /// to the commonly used generators.
   8.902 -  class Random {
   8.903 -  public:
   8.904 -
   8.905 -    static const unsigned long min = 0ul;
   8.906 -    static const unsigned long max = ~0ul;
   8.907 -  
   8.908 -    /// \brief Constructor
   8.909 -    ///
   8.910 -    /// Constructor with time dependent seeding.
   8.911 -    Random() {}
   8.912 -
   8.913 -    /// \brief Constructor
   8.914 -    ///
   8.915 -    /// Constructor
   8.916 -    Random(unsigned long seed) {}
   8.917 -
   8.918 -    /// \brief Copy constructor
   8.919 -    ///
   8.920 -    /// Copy constructor. The generated sequence will be identical to
   8.921 -    /// the other sequence.
   8.922 -    Random(const Random& other) {}
   8.923 -
   8.924 -    /// \brief Assign operator
   8.925 -    ///
   8.926 -    /// Assign operator. The generated sequence will be identical to
   8.927 -    /// the other sequence.
   8.928 -    Random& operator=(const Random& other) { return *this; }
   8.929 -
   8.930 -    /// \brief Returns an unsigned random number
   8.931 -    ///
   8.932 -    /// It returns an unsigned integer random number from the range 
   8.933 -    /// \f$ \{0, 1, \dots, 2^{64} - 1\} \f$ for 64-bit processors and from
   8.934 -    /// \f$ \{0, 1, \dots, 2^{32} - 1\} \f$ for 32-bit processors.
   8.935 -    unsigned long getUnsigned() { return 0ul; }
   8.936 -
   8.937 -    /// \brief Returns a random number
   8.938 -    ///
   8.939 -    /// It returns an integer random number from the range 
   8.940 -    /// \f$ \{-2^{63}, \dots, 2^{63} - 1\} \f$ for 64-bit processors and from
   8.941 -    /// \f$ \{-2^{31}, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
   8.942 -    long getInt() { return 0l; }
   8.943 -    
   8.944 -    /// \brief Returns an unsigned integer random number
   8.945 -    ///
   8.946 -    /// It returns an unsigned integer random number from the range 
   8.947 -    /// \f$ \{0, 1, \dots, 2^{63} - 1\} \f$ for 64-bit processors and
   8.948 -    /// from \f$ \{0, 1, \dots, 2^{31} - 1\} \f$ for 32-bit processors.
   8.949 -    long getNatural() { return 0l; }
   8.950 -
   8.951 -    /// \brief Returns a random bool
   8.952 -    ///
   8.953 -    /// It returns a random bool.
   8.954 -    bool getBool() { return false; }
   8.955 -
   8.956 -    /// \brief Returns a real random number
   8.957 -    ///
   8.958 -    /// It returns a real random number from the range 
   8.959 -    /// \f$ [0, 1) \f$. For 32-bit processors the generated random
   8.960 -    /// number will just have 32 significant bits.
   8.961 -    double getReal() { return 0.0; }
   8.962 -
   8.963 -    /// \brief Returns a real random number
   8.964 -    ///
   8.965 -    /// It returns a real random number from the range 
   8.966 -    /// \f$ [0, 1) \f$. This function returns random numbers with 53
   8.967 -    /// significant bits for 32-bit processors. For 64-bit processors
   8.968 -    /// it is identical to the \c getReal().
   8.969 -    double getPrecReal() { return 0.0; }
   8.970 -
   8.971 -    /// \brief Returns an unsigned random number from a given range
   8.972 -    ///
   8.973 -    /// It returns an unsigned integer random number from the range 
   8.974 -    /// \f$ \{0, 1, \dots, n - 1\} \f$.
   8.975 -    unsigned long getUnsigned(unsigned long n) { return 0; }
   8.976 -
   8.977 -    /// \brief Returns a random number from a given range
   8.978 -    ///
   8.979 -    /// It returns an unsigned integer random number from the range 
   8.980 -    /// \f$ \{0, 1, \dots, n - 1\} \f$.
   8.981 -    long getInt(long n) { return 0; }
   8.982 -
   8.983 -  };
   8.984 -
   8.985 -#endif
   8.986  
   8.987    /// \brief Global random number generator instance
   8.988    ///
     9.1 --- a/lemon/simann.h	Fri Oct 13 15:10:50 2006 +0000
     9.2 +++ b/lemon/simann.h	Sat Oct 14 15:26:05 2006 +0000
     9.3 @@ -257,7 +257,7 @@
     9.4      /// \brief Decides whether to accept the current solution or not.
     9.5      bool accept() {
     9.6        double cost_diff = simann->getCurrCost() - simann->getPrevCost();
     9.7 -      return (rnd.getReal() <= exp(-(cost_diff / temp)));
     9.8 +      return (rnd() <= exp(-(cost_diff / temp)));
     9.9      }
    9.10      /// \brief Destructor.
    9.11      virtual ~SimpleController() {}
    9.12 @@ -365,7 +365,7 @@
    9.13        }
    9.14        else {
    9.15          double cost_diff = simann->getCurrCost() - simann->getPrevCost();
    9.16 -        return (rnd.getReal() <= exp(-(cost_diff / temp)));
    9.17 +        return (rnd() <= exp(-(cost_diff / temp)));
    9.18        }
    9.19      }
    9.20      /// \brief Destructor.
    10.1 --- a/test/all_pairs_shortest_path_test.cc	Fri Oct 13 15:10:50 2006 +0000
    10.2 +++ b/test/all_pairs_shortest_path_test.cc	Sat Oct 14 15:26:05 2006 +0000
    10.3 @@ -36,19 +36,18 @@
    10.4  using namespace std;
    10.5  
    10.6  int main(int argc, const char *argv[]) {
    10.7 -  srand(time(0));
    10.8    typedef SmartGraph Graph;
    10.9    typedef Graph::Node Node;
   10.10    typedef Graph::Edge Edge;
   10.11    typedef Graph::NodeIt NodeIt;
   10.12    typedef Graph::EdgeIt EdgeIt;
   10.13  
   10.14 -  typedef Graph::EdgeMap<double> LengthMap;
   10.15 -  typedef Graph::NodeMap<double> DistMap;
   10.16 +  typedef Graph::EdgeMap<int> LengthMap;
   10.17 +  typedef Graph::NodeMap<int> DistMap;
   10.18  
   10.19    const int n = argc > 1 ? atoi(argv[1]) : 20;
   10.20    const int e = argc > 2 ? atoi(argv[2]) : (int)(n * log((double)n));
   10.21 -  const double m = argc > 3 ? (double)atoi(argv[3]) : 100.0;
   10.22 +  const int m = argc > 3 ? atoi(argv[3]) : 100;
   10.23  
   10.24    Graph graph;
   10.25    LengthMap length(graph);
   10.26 @@ -58,15 +57,14 @@
   10.27    for (int i = 0; i < n; ++i) {
   10.28      Node node = graph.addNode();
   10.29      nodes.push_back(node);
   10.30 -    shift[node] = m * (double)rand() / (RAND_MAX + 1.0);
   10.31 +    shift[node] = rnd[m];
   10.32    }
   10.33  
   10.34    for (int i = 0; i < e; ++i) {
   10.35 -    int s = (int)(n * (double)rand() / (RAND_MAX + 1.0));
   10.36 -    int t = (int)(n * (double)rand() / (RAND_MAX + 1.0));
   10.37 -    double c = m * (double)rand() / (RAND_MAX + 1.0);
   10.38 +    int s = rnd[n];
   10.39 +    int t = rnd[n];
   10.40      Edge edge = graph.addEdge(nodes[s], nodes[t]);
   10.41 -    length[edge] = c - shift[nodes[s]] + shift[nodes[t]];
   10.42 +    length[edge] = rnd[m] - shift[nodes[s]] + shift[nodes[t]];
   10.43    }
   10.44  
   10.45    Johnson<Graph, LengthMap> johnson(graph, length);
   10.46 @@ -76,8 +74,8 @@
   10.47      cout << "Johnson: " << timer << endl;
   10.48    }
   10.49  
   10.50 -  typedef FibHeap<Node, double, Graph::NodeMap<int> > DoubleFibHeap;
   10.51 -  Johnson<Graph, LengthMap>::DefStandardHeap<DoubleFibHeap>
   10.52 +  typedef FibHeap<Node, int, Graph::NodeMap<int> > IntFibHeap;
   10.53 +  Johnson<Graph, LengthMap>::DefStandardHeap<IntFibHeap>
   10.54      ::Create fibJohnson(graph, length);
   10.55    {
   10.56      Timer timer;
    11.1 --- a/test/arborescence_test.cc	Fri Oct 13 15:10:50 2006 +0000
    11.2 +++ b/test/arborescence_test.cc	Sat Oct 14 15:26:05 2006 +0000
    11.3 @@ -47,7 +47,6 @@
    11.4  
    11.5  
    11.6  int main() {
    11.7 -  srand(time(0));
    11.8    typedef SmartGraph Graph;
    11.9    GRAPH_TYPEDEFS(Graph);
   11.10  
    12.1 --- a/test/graph_utils_test.cc	Fri Oct 13 15:10:50 2006 +0000
    12.2 +++ b/test/graph_utils_test.cc	Sat Oct 14 15:26:05 2006 +0000
    12.3 @@ -116,7 +116,7 @@
    12.4      std::vector<ListGraph::Edge> edges(edgeNum);
    12.5      for (int i = 0; i < edgeNum; ++i) {
    12.6        edges[i] = 
    12.7 -	graph.addEdge(nodes[urandom(nodeNum)], nodes[urandom(nodeNum)]);
    12.8 +	graph.addEdge(nodes[rnd[nodeNum]], nodes[rnd[nodeNum]]);
    12.9      }
   12.10      for (int i = 0; i < nodeNum; ++i) {
   12.11        check(inDeg[nodes[i]] == countInEdges(graph, nodes[i]), 
    13.1 --- a/test/graph_utils_test.h	Fri Oct 13 15:10:50 2006 +0000
    13.2 +++ b/test/graph_utils_test.h	Sat Oct 14 15:26:05 2006 +0000
    13.3 @@ -56,8 +56,8 @@
    13.4      DescriptorMap<Graph, Node> nodes(graph);
    13.5      typename DescriptorMap<Graph, Node>::InverseMap invNodes(nodes);
    13.6      for (int i = 0; i < 100; ++i) {
    13.7 -      int src = rnd.getInt(invNodes.size());
    13.8 -      int trg = rnd.getInt(invNodes.size());
    13.9 +      int src = rnd[invNodes.size()];
   13.10 +      int trg = rnd[invNodes.size()];
   13.11        graph.addEdge(invNodes[src], invNodes[trg]);
   13.12      }
   13.13      typename Graph::template EdgeMap<bool> found(graph, false);
    14.1 --- a/test/heap_test.h	Fri Oct 13 15:10:50 2006 +0000
    14.2 +++ b/test/heap_test.h	Sat Oct 14 15:26:05 2006 +0000
    14.3 @@ -45,7 +45,7 @@
    14.4    std::vector<int> v(n);
    14.5  
    14.6    for (int i = 0; i < n; ++i) {
    14.7 -    v[i] = rnd.getInt(1000);
    14.8 +    v[i] = rnd[1000];
    14.9      heap.push(i, v[i]);
   14.10    }
   14.11    std::sort(v.begin(), v.end());
   14.12 @@ -65,11 +65,11 @@
   14.13    std::vector<int> v(n);
   14.14  
   14.15    for (int i = 0; i < n; ++i) {
   14.16 -    v[i] = rnd.getInt(1000);
   14.17 +    v[i] = rnd[1000];
   14.18      heap.push(i, v[i]);
   14.19    }
   14.20    for (int i = 0; i < n; ++i) {
   14.21 -    v[i] += rnd.getInt(1000);
   14.22 +    v[i] += rnd[1000];
   14.23      heap.increase(i, v[i]);
   14.24    }
   14.25    std::sort(v.begin(), v.end());
    15.1 --- a/test/matrix_maps_test.cc	Fri Oct 13 15:10:50 2006 +0000
    15.2 +++ b/test/matrix_maps_test.cc	Sat Oct 14 15:26:05 2006 +0000
    15.3 @@ -67,7 +67,7 @@
    15.4      }
    15.5      for (Graph::NodeIt it(graph); it != INVALID; ++it) {
    15.6        for (Graph::NodeIt jt(graph); jt != INVALID; ++jt) {
    15.7 -	int val = urandom(100);
    15.8 +	int val = rnd[100];
    15.9  	matrix.set(it, jt, val);
   15.10  	check(matrix(it, jt) == val, "Wrong assign");
   15.11  	check(matrix(it, jt) == matrixRowMap(matrix, it)[jt], "Wrong rowMap");
   15.12 @@ -108,7 +108,7 @@
   15.13      }
   15.14      for (Graph::NodeIt it(graph); it != INVALID; ++it) {
   15.15        for (Graph::NodeIt jt(graph); jt != INVALID; ++jt) {
   15.16 -	int val = urandom(100);
   15.17 +	int val = rnd[100];
   15.18  	matrix.set(it, jt, val);
   15.19  	check(matrix(it, jt) == val, "Wrong assign");
   15.20  	check(matrix(jt, it) == val, "Wrong assign");
   15.21 @@ -154,7 +154,7 @@
   15.22      }
   15.23      for (Graph::NodeIt it(graph1); it != INVALID; ++it) {
   15.24        for (Graph::EdgeIt jt(graph2); jt != INVALID; ++jt) {
   15.25 -	int val = urandom(100);
   15.26 +	int val = rnd[100];
   15.27  	matrix.set(it, jt, val);
   15.28  	check(matrix(it, jt) == val, "Wrong assign");
   15.29  	check(matrix(it, jt) == matrixRowMap(matrix, it)[jt], "Wrong rowMap");
   15.30 @@ -198,7 +198,7 @@
   15.31      }
   15.32      for (Graph::NodeIt it(graph); it != INVALID; ++it) {
   15.33        for (Graph::NodeIt jt(graph); jt != INVALID; ++jt) {
   15.34 -	int val = urandom(100);
   15.35 +	int val = rnd[100];
   15.36  	matrix.set(it, jt, val);
   15.37  	check(matrix(it, jt) == val, "Wrong assign");
   15.38  	check(matrix(it, jt) == matrixRowMap(matrix, it)[jt], "Wrong rowMap");
    16.1 --- a/test/radix_sort_test.cc	Fri Oct 13 15:10:50 2006 +0000
    16.2 +++ b/test/radix_sort_test.cc	Sat Oct 14 15:26:05 2006 +0000
    16.3 @@ -37,7 +37,7 @@
    16.4      int n = 10000;
    16.5      vector<int> data1(n), data2(n);
    16.6      for (int i = 0; i < n; ++i) {
    16.7 -      data1[i] = data2[i] = (int)(1000 * (rand() / (RAND_MAX + 1.0))) - 500;
    16.8 +      data1[i] = data2[i] = rnd[1000] - 500;
    16.9      }
   16.10      radixSort(data1.begin(), data1.end());
   16.11      sort(data2.begin(), data2.end());
   16.12 @@ -48,7 +48,7 @@
   16.13      int n = 10000;
   16.14      vector<unsigned char> data1(n), data2(n);
   16.15      for (int i = 0; i < n; ++i) {
   16.16 -      data1[i] = data2[i] = (int)(200 * (rand() / (RAND_MAX + 1.0)));
   16.17 +      data1[i] = data2[i] = rnd[(unsigned char)200];
   16.18      }
   16.19      radixSort(data1.begin(), data1.end());
   16.20      sort(data2.begin(), data2.end());
   16.21 @@ -64,7 +64,7 @@
   16.22      int n = 10000;
   16.23      vector<int> data1(n), data2(n);
   16.24      for (int i = 0; i < n; ++i) {
   16.25 -      data1[i] = data2[i] = (int)(1000 * (rand() / (RAND_MAX + 1.0))) - 500;
   16.26 +      data1[i] = data2[i] = rnd[1000] - 500;
   16.27      }
   16.28      counterSort(data1.begin(), data1.end());
   16.29      sort(data2.begin(), data2.end());
   16.30 @@ -75,7 +75,7 @@
   16.31      int n = 10000;
   16.32      vector<unsigned char> data1(n), data2(n);
   16.33      for (int i = 0; i < n; ++i) {
   16.34 -      data1[i] = data2[i] = (int)(200 * (rand() / (RAND_MAX + 1.0)));
   16.35 +      data1[i] = data2[i] = rnd[(unsigned char)200];
   16.36      }
   16.37      counterSort(data1.begin(), data1.end());
   16.38      sort(data2.begin(), data2.end());
   16.39 @@ -106,8 +106,8 @@
   16.40    }
   16.41    vector<Edge> edges;
   16.42    for (int i = 0; i < e; ++i) {
   16.43 -    int s = (int)(n * (double)rand() / (RAND_MAX + 1.0));
   16.44 -    int t = (int)(n * (double)rand() / (RAND_MAX + 1.0));
   16.45 +    int s = rnd[n];
   16.46 +    int t = rnd[n];
   16.47      edges.push_back(graph.addEdge(nodes[s], nodes[t]));
   16.48    }
   16.49  
   16.50 @@ -145,8 +145,8 @@
   16.51    }
   16.52    vector<Edge> edges;
   16.53    for (int i = 0; i < e; ++i) {
   16.54 -    int s = (int)(n * (double)rand() / (RAND_MAX + 1.0));
   16.55 -    int t = (int)(n * (double)rand() / (RAND_MAX + 1.0));
   16.56 +    int s = rnd[n];
   16.57 +    int t = rnd[n];
   16.58      edges.push_back(graph.addEdge(nodes[s], nodes[t]));
   16.59    }
   16.60  
    17.1 --- a/test/simann_test.cc	Fri Oct 13 15:10:50 2006 +0000
    17.2 +++ b/test/simann_test.cc	Sat Oct 14 15:26:05 2006 +0000
    17.3 @@ -23,10 +23,10 @@
    17.4  class MyEntity : public EntityBase {
    17.5  public:
    17.6    double d, prev_d;
    17.7 -  MyEntity() : d(100000.0) { srand48(time(0)); }
    17.8 +  MyEntity() : d(100000.0) {}
    17.9    double mutate() {
   17.10      prev_d = d;
   17.11 -    if (drand48() < 0.8) { d += 1.0; }
   17.12 +    if (rnd.boolean(0.8)) { d += 1.0; }
   17.13      else { d -= 1.0; }
   17.14      return d;
   17.15    }
    18.1 --- a/test/test_tools.h	Fri Oct 13 15:10:50 2006 +0000
    18.2 +++ b/test/test_tools.h	Sat Oct 14 15:26:05 2006 +0000
    18.3 @@ -174,12 +174,8 @@
    18.4     n.chords.push_back(G.addEdge(n.outer[i],n.inner[i]));
    18.5     n.outcir.push_back(G.addEdge(n.outer[i],n.outer[(i+1)%5]));
    18.6     n.incir.push_back(G.addEdge(n.inner[i],n.inner[(i+2)%5]));
    18.7 -  }
    18.8 + }
    18.9   return n;
   18.10  }
   18.11  
   18.12 -int urandom(int n) {
   18.13 -  return rnd.getInt(n);
   18.14 -}
   18.15 -
   18.16  #endif