lemon/random.h
changeset 22 45f8b617339e
parent 12 435bbc8127b3
child 23 0ba375bf5dae
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/random.h	Thu Jan 03 11:13:29 2008 +0100
     1.3 @@ -0,0 +1,872 @@
     1.4 +/* -*- C++ -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library
     1.7 + *
     1.8 + * Copyright (C) 2003-2007
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +/*
    1.23 + * This file contains the reimplemented version of the Mersenne Twister
    1.24 + * Generator of Matsumoto and Nishimura.
    1.25 + *
    1.26 + * See the appropriate copyright notice below.
    1.27 + * 
    1.28 + * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
    1.29 + * All rights reserved.                          
    1.30 + *
    1.31 + * Redistribution and use in source and binary forms, with or without
    1.32 + * modification, are permitted provided that the following conditions
    1.33 + * are met:
    1.34 + *
    1.35 + * 1. Redistributions of source code must retain the above copyright
    1.36 + *    notice, this list of conditions and the following disclaimer.
    1.37 + *
    1.38 + * 2. Redistributions in binary form must reproduce the above copyright
    1.39 + *    notice, this list of conditions and the following disclaimer in the
    1.40 + *    documentation and/or other materials provided with the distribution.
    1.41 + *
    1.42 + * 3. The names of its contributors may not be used to endorse or promote 
    1.43 + *    products derived from this software without specific prior written 
    1.44 + *    permission.
    1.45 + *
    1.46 + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    1.47 + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    1.48 + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
    1.49 + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
    1.50 + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
    1.51 + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    1.52 + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
    1.53 + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
    1.54 + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
    1.55 + * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    1.56 + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
    1.57 + * OF THE POSSIBILITY OF SUCH DAMAGE.
    1.58 + *
    1.59 + *
    1.60 + * Any feedback is very welcome.
    1.61 + * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
    1.62 + * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
    1.63 + */
    1.64 +
    1.65 +#ifndef LEMON_RANDOM_H
    1.66 +#define LEMON_RANDOM_H
    1.67 +
    1.68 +#include <algorithm>
    1.69 +#include <iterator>
    1.70 +#include <vector>
    1.71 +
    1.72 +#include <ctime>
    1.73 +#include <cmath>
    1.74 +
    1.75 +#include <lemon/dim2.h>
    1.76 +///\ingroup misc
    1.77 +///\file
    1.78 +///\brief Mersenne Twister random number generator
    1.79 +
    1.80 +namespace lemon {
    1.81 +
    1.82 +  namespace _random_bits {
    1.83 +    
    1.84 +    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
    1.85 +    struct RandomTraits {};
    1.86 +
    1.87 +    template <typename _Word>
    1.88 +    struct RandomTraits<_Word, 32> {
    1.89 +
    1.90 +      typedef _Word Word;
    1.91 +      static const int bits = 32;
    1.92 +
    1.93 +      static const int length = 624;
    1.94 +      static const int shift = 397;
    1.95 +      
    1.96 +      static const Word mul = 0x6c078965u;
    1.97 +      static const Word arrayInit = 0x012BD6AAu;
    1.98 +      static const Word arrayMul1 = 0x0019660Du;
    1.99 +      static const Word arrayMul2 = 0x5D588B65u;
   1.100 +
   1.101 +      static const Word mask = 0x9908B0DFu;
   1.102 +      static const Word loMask = (1u << 31) - 1;
   1.103 +      static const Word hiMask = ~loMask;
   1.104 +
   1.105 +
   1.106 +      static Word tempering(Word rnd) {
   1.107 +        rnd ^= (rnd >> 11);
   1.108 +        rnd ^= (rnd << 7) & 0x9D2C5680u;
   1.109 +        rnd ^= (rnd << 15) & 0xEFC60000u;
   1.110 +        rnd ^= (rnd >> 18);
   1.111 +        return rnd;
   1.112 +      }
   1.113 +
   1.114 +    };
   1.115 +
   1.116 +    template <typename _Word>
   1.117 +    struct RandomTraits<_Word, 64> {
   1.118 +
   1.119 +      typedef _Word Word;
   1.120 +      static const int bits = 64;
   1.121 +
   1.122 +      static const int length = 312;
   1.123 +      static const int shift = 156;
   1.124 +
   1.125 +      static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du);
   1.126 +      static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu);
   1.127 +      static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u);
   1.128 +      static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu);
   1.129 +
   1.130 +      static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u);
   1.131 +      static const Word loMask = (Word(1u) << 31) - 1;
   1.132 +      static const Word hiMask = ~loMask;
   1.133 +
   1.134 +      static Word tempering(Word rnd) {
   1.135 +        rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u));
   1.136 +        rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u));
   1.137 +        rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u));
   1.138 +        rnd ^= (rnd >> 43);
   1.139 +        return rnd;
   1.140 +      }
   1.141 +
   1.142 +    };
   1.143 +
   1.144 +    template <typename _Word>
   1.145 +    class RandomCore {
   1.146 +    public:
   1.147 +
   1.148 +      typedef _Word Word;
   1.149 +
   1.150 +    private:
   1.151 +
   1.152 +      static const int bits = RandomTraits<Word>::bits;
   1.153 +
   1.154 +      static const int length = RandomTraits<Word>::length;
   1.155 +      static const int shift = RandomTraits<Word>::shift;
   1.156 +
   1.157 +    public:
   1.158 +
   1.159 +      void initState() {
   1.160 +        static const Word seedArray[4] = {
   1.161 +          0x12345u, 0x23456u, 0x34567u, 0x45678u
   1.162 +        };
   1.163 +    
   1.164 +        initState(seedArray, seedArray + 4);
   1.165 +      }
   1.166 +
   1.167 +      void initState(Word seed) {
   1.168 +
   1.169 +        static const Word mul = RandomTraits<Word>::mul;
   1.170 +
   1.171 +        current = state; 
   1.172 +
   1.173 +        Word *curr = state + length - 1;
   1.174 +        curr[0] = seed; --curr;
   1.175 +        for (int i = 1; i < length; ++i) {
   1.176 +          curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
   1.177 +          --curr;
   1.178 +        }
   1.179 +      }
   1.180 +
   1.181 +      template <typename Iterator>
   1.182 +      void initState(Iterator begin, Iterator end) {
   1.183 +
   1.184 +        static const Word init = RandomTraits<Word>::arrayInit;
   1.185 +        static const Word mul1 = RandomTraits<Word>::arrayMul1;
   1.186 +        static const Word mul2 = RandomTraits<Word>::arrayMul2;
   1.187 +
   1.188 +
   1.189 +        Word *curr = state + length - 1; --curr;
   1.190 +        Iterator it = begin; int cnt = 0;
   1.191 +        int num;
   1.192 +
   1.193 +        initState(init);
   1.194 +
   1.195 +        num = length > end - begin ? length : end - begin;
   1.196 +        while (num--) {
   1.197 +          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 
   1.198 +            + *it + cnt;
   1.199 +          ++it; ++cnt;
   1.200 +          if (it == end) {
   1.201 +            it = begin; cnt = 0;
   1.202 +          }
   1.203 +          if (curr == state) {
   1.204 +            curr = state + length - 1; curr[0] = state[0];
   1.205 +          }
   1.206 +          --curr;
   1.207 +        }
   1.208 +
   1.209 +        num = length - 1; cnt = length - (curr - state) - 1;
   1.210 +        while (num--) {
   1.211 +          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
   1.212 +            - cnt;
   1.213 +          --curr; ++cnt;
   1.214 +          if (curr == state) {
   1.215 +            curr = state + length - 1; curr[0] = state[0]; --curr;
   1.216 +            cnt = 1;
   1.217 +          }
   1.218 +        }
   1.219 +        
   1.220 +        state[length - 1] = Word(1) << (bits - 1);
   1.221 +      }
   1.222 +      
   1.223 +      void copyState(const RandomCore& other) {
   1.224 +        std::copy(other.state, other.state + length, state);
   1.225 +        current = state + (other.current - other.state);
   1.226 +      }
   1.227 +
   1.228 +      Word operator()() {
   1.229 +        if (current == state) fillState();
   1.230 +        --current;
   1.231 +        Word rnd = *current;
   1.232 +        return RandomTraits<Word>::tempering(rnd);
   1.233 +      }
   1.234 +
   1.235 +    private:
   1.236 +
   1.237 +  
   1.238 +      void fillState() {
   1.239 +        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
   1.240 +        static const Word loMask = RandomTraits<Word>::loMask;
   1.241 +        static const Word hiMask = RandomTraits<Word>::hiMask;
   1.242 +
   1.243 +        current = state + length; 
   1.244 +
   1.245 +        register Word *curr = state + length - 1;
   1.246 +        register long num;
   1.247 +      
   1.248 +        num = length - shift;
   1.249 +        while (num--) {
   1.250 +          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   1.251 +            curr[- shift] ^ mask[curr[-1] & 1ul];
   1.252 +          --curr;
   1.253 +        }
   1.254 +        num = shift - 1;
   1.255 +        while (num--) {
   1.256 +          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   1.257 +            curr[length - shift] ^ mask[curr[-1] & 1ul];
   1.258 +          --curr;
   1.259 +        }
   1.260 +        curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   1.261 +          curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   1.262 +
   1.263 +      }
   1.264 +
   1.265 +  
   1.266 +      Word *current;
   1.267 +      Word state[length];
   1.268 +      
   1.269 +    };
   1.270 +
   1.271 +
   1.272 +    template <typename Result, 
   1.273 +              int shift = (std::numeric_limits<Result>::digits + 1) / 2>
   1.274 +    struct Masker {
   1.275 +      static Result mask(const Result& result) {
   1.276 +        return Masker<Result, (shift + 1) / 2>::
   1.277 +          mask(static_cast<Result>(result | (result >> shift)));
   1.278 +      }
   1.279 +    };
   1.280 +    
   1.281 +    template <typename Result>
   1.282 +    struct Masker<Result, 1> {
   1.283 +      static Result mask(const Result& result) {
   1.284 +        return static_cast<Result>(result | (result >> 1));
   1.285 +      }
   1.286 +    };
   1.287 +
   1.288 +    template <typename Result, typename Word, 
   1.289 +              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   1.290 +              bool last = rest <= std::numeric_limits<Word>::digits>
   1.291 +    struct IntConversion {
   1.292 +      static const int bits = std::numeric_limits<Word>::digits;
   1.293 +    
   1.294 +      static Result convert(RandomCore<Word>& rnd) {
   1.295 +        return static_cast<Result>(rnd() >> (bits - rest)) << shift;
   1.296 +      }
   1.297 +      
   1.298 +    }; 
   1.299 +
   1.300 +    template <typename Result, typename Word, int rest, int shift> 
   1.301 +    struct IntConversion<Result, Word, rest, shift, false> {
   1.302 +      static const int bits = std::numeric_limits<Word>::digits;
   1.303 +
   1.304 +      static Result convert(RandomCore<Word>& rnd) {
   1.305 +        return (static_cast<Result>(rnd()) << shift) | 
   1.306 +          IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
   1.307 +      }
   1.308 +    };
   1.309 +
   1.310 +
   1.311 +    template <typename Result, typename Word,
   1.312 +              bool one_word = (std::numeric_limits<Word>::digits < 
   1.313 +			       std::numeric_limits<Result>::digits) >
   1.314 +    struct Mapping {
   1.315 +      static Result map(RandomCore<Word>& rnd, const Result& bound) {
   1.316 +        Word max = Word(bound - 1);
   1.317 +        Result mask = Masker<Result>::mask(bound - 1);
   1.318 +        Result num;
   1.319 +        do {
   1.320 +          num = IntConversion<Result, Word>::convert(rnd) & mask; 
   1.321 +        } while (num > max);
   1.322 +        return num;
   1.323 +      }
   1.324 +    };
   1.325 +
   1.326 +    template <typename Result, typename Word>
   1.327 +    struct Mapping<Result, Word, false> {
   1.328 +      static Result map(RandomCore<Word>& rnd, const Result& bound) {
   1.329 +        Word max = Word(bound - 1);
   1.330 +        Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
   1.331 +          ::mask(max);
   1.332 +        Word num;
   1.333 +        do {
   1.334 +          num = rnd() & mask;
   1.335 +        } while (num > max);
   1.336 +        return num;
   1.337 +      }
   1.338 +    };
   1.339 +
   1.340 +    template <typename Result, int exp, bool pos = (exp >= 0)>
   1.341 +    struct ShiftMultiplier {
   1.342 +      static const Result multiplier() {
   1.343 +        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   1.344 +        res *= res;
   1.345 +        if ((exp & 1) == 1) res *= static_cast<Result>(2.0);
   1.346 +        return res; 
   1.347 +      }
   1.348 +    };
   1.349 +
   1.350 +    template <typename Result, int exp>
   1.351 +    struct ShiftMultiplier<Result, exp, false> {
   1.352 +      static const Result multiplier() {
   1.353 +        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   1.354 +        res *= res;
   1.355 +        if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
   1.356 +        return res; 
   1.357 +      }
   1.358 +    };
   1.359 +
   1.360 +    template <typename Result>
   1.361 +    struct ShiftMultiplier<Result, 0, true> {
   1.362 +      static const Result multiplier() {
   1.363 +        return static_cast<Result>(1.0); 
   1.364 +      }
   1.365 +    };
   1.366 +
   1.367 +    template <typename Result>
   1.368 +    struct ShiftMultiplier<Result, -20, true> {
   1.369 +      static const Result multiplier() {
   1.370 +        return static_cast<Result>(1.0/1048576.0); 
   1.371 +      }
   1.372 +    };
   1.373 +    
   1.374 +    template <typename Result>
   1.375 +    struct ShiftMultiplier<Result, -32, true> {
   1.376 +      static const Result multiplier() {
   1.377 +        return static_cast<Result>(1.0/424967296.0); 
   1.378 +      }
   1.379 +    };
   1.380 +
   1.381 +    template <typename Result>
   1.382 +    struct ShiftMultiplier<Result, -53, true> {
   1.383 +      static const Result multiplier() {
   1.384 +        return static_cast<Result>(1.0/9007199254740992.0); 
   1.385 +      }
   1.386 +    };
   1.387 +
   1.388 +    template <typename Result>
   1.389 +    struct ShiftMultiplier<Result, -64, true> {
   1.390 +      static const Result multiplier() {
   1.391 +        return static_cast<Result>(1.0/18446744073709551616.0); 
   1.392 +      }
   1.393 +    };
   1.394 +
   1.395 +    template <typename Result, int exp>
   1.396 +    struct Shifting {
   1.397 +      static Result shift(const Result& result) {
   1.398 +        return result * ShiftMultiplier<Result, exp>::multiplier();
   1.399 +      }
   1.400 +    };
   1.401 +
   1.402 +    template <typename Result, typename Word,
   1.403 +              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   1.404 +              bool last = rest <= std::numeric_limits<Word>::digits>
   1.405 +    struct RealConversion{ 
   1.406 +      static const int bits = std::numeric_limits<Word>::digits;
   1.407 +
   1.408 +      static Result convert(RandomCore<Word>& rnd) {
   1.409 +        return Shifting<Result, - shift - rest>::
   1.410 +          shift(static_cast<Result>(rnd() >> (bits - rest)));
   1.411 +      }
   1.412 +    };
   1.413 +
   1.414 +    template <typename Result, typename Word, int rest, int shift>
   1.415 +    struct RealConversion<Result, Word, rest, shift, false> { 
   1.416 +      static const int bits = std::numeric_limits<Word>::digits;
   1.417 +
   1.418 +      static Result convert(RandomCore<Word>& rnd) {
   1.419 +        return Shifting<Result, - shift - bits>::
   1.420 +          shift(static_cast<Result>(rnd())) +
   1.421 +          RealConversion<Result, Word, rest-bits, shift + bits>::
   1.422 +          convert(rnd);
   1.423 +      }
   1.424 +    };
   1.425 +
   1.426 +    template <typename Result, typename Word>
   1.427 +    struct Initializer {
   1.428 +
   1.429 +      template <typename Iterator>
   1.430 +      static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
   1.431 +        std::vector<Word> ws;
   1.432 +        for (Iterator it = begin; it != end; ++it) {
   1.433 +          ws.push_back(Word(*it));
   1.434 +        }
   1.435 +        rnd.initState(ws.begin(), ws.end());
   1.436 +      }
   1.437 +
   1.438 +      static void init(RandomCore<Word>& rnd, Result seed) {
   1.439 +        rnd.initState(seed);
   1.440 +      }
   1.441 +    };
   1.442 +
   1.443 +    template <typename Word>
   1.444 +    struct BoolConversion {
   1.445 +      static bool convert(RandomCore<Word>& rnd) {
   1.446 +        return (rnd() & 1) == 1;
   1.447 +      }
   1.448 +    };
   1.449 +
   1.450 +    template <typename Word>
   1.451 +    struct BoolProducer {
   1.452 +      Word buffer;
   1.453 +      int num;
   1.454 +      
   1.455 +      BoolProducer() : num(0) {}
   1.456 +
   1.457 +      bool convert(RandomCore<Word>& rnd) {
   1.458 +        if (num == 0) {
   1.459 +          buffer = rnd();
   1.460 +          num = RandomTraits<Word>::bits;
   1.461 +        }
   1.462 +        bool r = (buffer & 1);
   1.463 +        buffer >>= 1;
   1.464 +        --num;
   1.465 +        return r;
   1.466 +      }
   1.467 +    };
   1.468 +
   1.469 +  }
   1.470 +
   1.471 +  /// \ingroup misc
   1.472 +  ///
   1.473 +  /// \brief Mersenne Twister random number generator
   1.474 +  ///
   1.475 +  /// The Mersenne Twister is a twisted generalized feedback
   1.476 +  /// shift-register generator of Matsumoto and Nishimura. The period
   1.477 +  /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
   1.478 +  /// equi-distributed in 623 dimensions for 32-bit numbers. The time
   1.479 +  /// performance of this generator is comparable to the commonly used
   1.480 +  /// generators.
   1.481 +  ///
   1.482 +  /// This implementation is specialized for both 32-bit and 64-bit
   1.483 +  /// architectures. The generators differ sligthly in the
   1.484 +  /// initialization and generation phase so they produce two
   1.485 +  /// completly different sequences.
   1.486 +  ///
   1.487 +  /// The generator gives back random numbers of serveral types. To
   1.488 +  /// get a random number from a range of a floating point type you
   1.489 +  /// can use one form of the \c operator() or the \c real() member
   1.490 +  /// function. If you want to get random number from the {0, 1, ...,
   1.491 +  /// n-1} integer range use the \c operator[] or the \c integer()
   1.492 +  /// method. And to get random number from the whole range of an
   1.493 +  /// integer type you can use the argumentless \c integer() or \c
   1.494 +  /// uinteger() functions. After all you can get random bool with
   1.495 +  /// equal chance of true and false or given probability of true
   1.496 +  /// result with the \c boolean() member functions.
   1.497 +  ///
   1.498 +  ///\code
   1.499 +  /// // The commented code is identical to the other
   1.500 +  /// double a = rnd();                     // [0.0, 1.0)
   1.501 +  /// // double a = rnd.real();             // [0.0, 1.0)
   1.502 +  /// double b = rnd(100.0);                // [0.0, 100.0)
   1.503 +  /// // double b = rnd.real(100.0);        // [0.0, 100.0)
   1.504 +  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
   1.505 +  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
   1.506 +  /// int d = rnd[100000];                  // 0..99999
   1.507 +  /// // int d = rnd.integer(100000);       // 0..99999
   1.508 +  /// int e = rnd[6] + 1;                   // 1..6
   1.509 +  /// // int e = rnd.integer(1, 1 + 6);     // 1..6
   1.510 +  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
   1.511 +  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
   1.512 +  /// bool g = rnd.boolean();               // P(g = true) = 0.5
   1.513 +  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
   1.514 +  ///\endcode
   1.515 +  ///
   1.516 +  /// The lemon provides a global instance of the random number
   1.517 +  /// generator which name is \ref lemon::rnd "rnd". Usually it is a
   1.518 +  /// good programming convenience to use this global generator to get
   1.519 +  /// random numbers.
   1.520 +  class Random {
   1.521 +  private:
   1.522 +
   1.523 +    // Architecture word
   1.524 +    typedef unsigned long Word;
   1.525 +    
   1.526 +    _random_bits::RandomCore<Word> core;
   1.527 +    _random_bits::BoolProducer<Word> bool_producer;
   1.528 +    
   1.529 +
   1.530 +  public:
   1.531 +
   1.532 +    /// \brief Constructor
   1.533 +    ///
   1.534 +    /// Constructor with constant seeding.
   1.535 +    Random() { core.initState(); }
   1.536 +
   1.537 +    /// \brief Constructor
   1.538 +    ///
   1.539 +    /// Constructor with seed. The current number type will be converted
   1.540 +    /// to the architecture word type.
   1.541 +    template <typename Number>
   1.542 +    Random(Number seed) { 
   1.543 +      _random_bits::Initializer<Number, Word>::init(core, seed);
   1.544 +    }
   1.545 +
   1.546 +    /// \brief Constructor
   1.547 +    ///
   1.548 +    /// Constructor with array seeding. The given range should contain
   1.549 +    /// any number type and the numbers will be converted to the
   1.550 +    /// architecture word type.
   1.551 +    template <typename Iterator>
   1.552 +    Random(Iterator begin, Iterator end) { 
   1.553 +      typedef typename std::iterator_traits<Iterator>::value_type Number;
   1.554 +      _random_bits::Initializer<Number, Word>::init(core, begin, end);
   1.555 +    }
   1.556 +
   1.557 +    /// \brief Copy constructor
   1.558 +    ///
   1.559 +    /// Copy constructor. The generated sequence will be identical to
   1.560 +    /// the other sequence. It can be used to save the current state
   1.561 +    /// of the generator and later use it to generate the same
   1.562 +    /// sequence.
   1.563 +    Random(const Random& other) {
   1.564 +      core.copyState(other.core);
   1.565 +    }
   1.566 +
   1.567 +    /// \brief Assign operator
   1.568 +    ///
   1.569 +    /// Assign operator. The generated sequence will be identical to
   1.570 +    /// the other sequence. It can be used to save the current state
   1.571 +    /// of the generator and later use it to generate the same
   1.572 +    /// sequence.
   1.573 +    Random& operator=(const Random& other) {
   1.574 +      if (&other != this) {
   1.575 +        core.copyState(other.core);
   1.576 +      }
   1.577 +      return *this;
   1.578 +    }
   1.579 +
   1.580 +    /// \brief Returns a random real number from the range [0, 1)
   1.581 +    ///
   1.582 +    /// It returns a random real number from the range [0, 1). The
   1.583 +    /// default Number type is double.
   1.584 +    template <typename Number>
   1.585 +    Number real() {
   1.586 +      return _random_bits::RealConversion<Number, Word>::convert(core);
   1.587 +    }
   1.588 +
   1.589 +    double real() {
   1.590 +      return real<double>();
   1.591 +    }
   1.592 +
   1.593 +    /// \brief Returns a random real number the range [0, b)
   1.594 +    ///
   1.595 +    /// It returns a random real number from the range [0, b).
   1.596 +    template <typename Number>
   1.597 +    Number real(Number b) { 
   1.598 +      return real<Number>() * b; 
   1.599 +    }
   1.600 +
   1.601 +    /// \brief Returns a random real number from the range [a, b)
   1.602 +    ///
   1.603 +    /// It returns a random real number from the range [a, b).
   1.604 +    template <typename Number>
   1.605 +    Number real(Number a, Number b) { 
   1.606 +      return real<Number>() * (b - a) + a; 
   1.607 +    }
   1.608 +
   1.609 +    /// \brief Returns a random real number from the range [0, 1)
   1.610 +    ///
   1.611 +    /// It returns a random double from the range [0, 1).
   1.612 +    double operator()() {
   1.613 +      return real<double>();
   1.614 +    }
   1.615 +
   1.616 +    /// \brief Returns a random real number from the range [0, b)
   1.617 +    ///
   1.618 +    /// It returns a random real number from the range [0, b).
   1.619 +    template <typename Number>
   1.620 +    Number operator()(Number b) { 
   1.621 +      return real<Number>() * b; 
   1.622 +    }
   1.623 +
   1.624 +    /// \brief Returns a random real number from the range [a, b)
   1.625 +    ///
   1.626 +    /// It returns a random real number from the range [a, b).
   1.627 +    template <typename Number>
   1.628 +    Number operator()(Number a, Number b) { 
   1.629 +      return real<Number>() * (b - a) + a; 
   1.630 +    }
   1.631 +
   1.632 +    /// \brief Returns a random integer from a range
   1.633 +    ///
   1.634 +    /// It returns a random integer from the range {0, 1, ..., b - 1}.
   1.635 +    template <typename Number>
   1.636 +    Number integer(Number b) {
   1.637 +      return _random_bits::Mapping<Number, Word>::map(core, b);
   1.638 +    }
   1.639 +
   1.640 +    /// \brief Returns a random integer from a range
   1.641 +    ///
   1.642 +    /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
   1.643 +    template <typename Number>
   1.644 +    Number integer(Number a, Number b) {
   1.645 +      return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
   1.646 +    }
   1.647 +
   1.648 +    /// \brief Returns a random integer from a range
   1.649 +    ///
   1.650 +    /// It returns a random integer from the range {0, 1, ..., b - 1}.
   1.651 +    template <typename Number>
   1.652 +    Number operator[](Number b) {
   1.653 +      return _random_bits::Mapping<Number, Word>::map(core, b);
   1.654 +    }
   1.655 +
   1.656 +    /// \brief Returns a random non-negative integer
   1.657 +    ///
   1.658 +    /// It returns a random non-negative integer uniformly from the
   1.659 +    /// whole range of the current \c Number type.  The default result
   1.660 +    /// type of this function is unsigned int.
   1.661 +    template <typename Number>
   1.662 +    Number uinteger() {
   1.663 +      return _random_bits::IntConversion<Number, Word>::convert(core);
   1.664 +    }
   1.665 +
   1.666 +    unsigned int uinteger() {
   1.667 +      return uinteger<unsigned int>();
   1.668 +    }
   1.669 +
   1.670 +    /// \brief Returns a random integer
   1.671 +    ///
   1.672 +    /// It returns a random integer uniformly from the whole range of
   1.673 +    /// the current \c Number type. The default result type of this
   1.674 +    /// function is int.
   1.675 +    template <typename Number>
   1.676 +    Number integer() {
   1.677 +      static const int nb = std::numeric_limits<Number>::digits + 
   1.678 +        (std::numeric_limits<Number>::is_signed ? 1 : 0);
   1.679 +      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
   1.680 +    }
   1.681 +
   1.682 +    int integer() {
   1.683 +      return integer<int>();
   1.684 +    }
   1.685 +    
   1.686 +    /// \brief Returns a random bool
   1.687 +    ///
   1.688 +    /// It returns a random bool. The generator holds a buffer for
   1.689 +    /// random bits. Every time when it become empty the generator makes
   1.690 +    /// a new random word and fill the buffer up.
   1.691 +    bool boolean() {
   1.692 +      return bool_producer.convert(core);
   1.693 +    }
   1.694 +
   1.695 +    ///\name Nonuniform distributions
   1.696 +    ///
   1.697 +    
   1.698 +    ///@{
   1.699 +    
   1.700 +    /// \brief Returns a random bool
   1.701 +    ///
   1.702 +    /// It returns a random bool with given probability of true result
   1.703 +    bool boolean(double p) {
   1.704 +      return operator()() < p;
   1.705 +    }
   1.706 +
   1.707 +    /// Standard Gauss distribution
   1.708 +
   1.709 +    /// Standard Gauss distribution.
   1.710 +    /// \note The Cartesian form of the Box-Muller
   1.711 +    /// transformation is used to generate a random normal distribution.
   1.712 +    /// \todo Consider using the "ziggurat" method instead.
   1.713 +    double gauss() 
   1.714 +    {
   1.715 +      double V1,V2,S;
   1.716 +      do {
   1.717 +	V1=2*real<double>()-1;
   1.718 +	V2=2*real<double>()-1;
   1.719 +	S=V1*V1+V2*V2;
   1.720 +      } while(S>=1);
   1.721 +      return std::sqrt(-2*std::log(S)/S)*V1;
   1.722 +    }
   1.723 +    /// Gauss distribution with given mean and standard deviation
   1.724 +
   1.725 +    /// Gauss distribution with given mean and standard deviation
   1.726 +    /// \sa gauss()
   1.727 +    double gauss(double mean,double std_dev)
   1.728 +    {
   1.729 +      return gauss()*std_dev+mean;
   1.730 +    }
   1.731 +
   1.732 +    /// Exponential distribution with given mean
   1.733 +
   1.734 +    /// This function generates an exponential distribution random number
   1.735 +    /// with mean <tt>1/lambda</tt>.
   1.736 +    ///
   1.737 +    double exponential(double lambda=1.0)
   1.738 +    {
   1.739 +      return -std::log(1.0-real<double>())/lambda;
   1.740 +    }
   1.741 +
   1.742 +    /// Gamma distribution with given integer shape
   1.743 +
   1.744 +    /// This function generates a gamma distribution random number.
   1.745 +    /// 
   1.746 +    ///\param k shape parameter (<tt>k>0</tt> integer)
   1.747 +    double gamma(int k) 
   1.748 +    {
   1.749 +      double s = 0;
   1.750 +      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
   1.751 +      return s;
   1.752 +    }
   1.753 +    
   1.754 +    /// Gamma distribution with given shape and scale parameter
   1.755 +
   1.756 +    /// This function generates a gamma distribution random number.
   1.757 +    /// 
   1.758 +    ///\param k shape parameter (<tt>k>0</tt>)
   1.759 +    ///\param theta scale parameter
   1.760 +    ///
   1.761 +    double gamma(double k,double theta=1.0)
   1.762 +    {
   1.763 +      double xi,nu;
   1.764 +      const double delta = k-std::floor(k);
   1.765 +      const double v0=M_E/(M_E-delta);
   1.766 +      do {
   1.767 +	double V0=1.0-real<double>();
   1.768 +	double V1=1.0-real<double>();
   1.769 +	double V2=1.0-real<double>();
   1.770 +	if(V2<=v0) 
   1.771 +	  {
   1.772 +	    xi=std::pow(V1,1.0/delta);
   1.773 +	    nu=V0*std::pow(xi,delta-1.0);
   1.774 +	  }
   1.775 +	else 
   1.776 +	  {
   1.777 +	    xi=1.0-std::log(V1);
   1.778 +	    nu=V0*std::exp(-xi);
   1.779 +	  }
   1.780 +      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
   1.781 +      return theta*(xi-gamma(int(std::floor(k))));
   1.782 +    }
   1.783 +    
   1.784 +    /// Weibull distribution
   1.785 +
   1.786 +    /// This function generates a Weibull distribution random number.
   1.787 +    /// 
   1.788 +    ///\param k shape parameter (<tt>k>0</tt>)
   1.789 +    ///\param lambda scale parameter (<tt>lambda>0</tt>)
   1.790 +    ///
   1.791 +    double weibull(double k,double lambda)
   1.792 +    {
   1.793 +      return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
   1.794 +    }  
   1.795 +      
   1.796 +    /// Pareto distribution
   1.797 +
   1.798 +    /// This function generates a Pareto distribution random number.
   1.799 +    /// 
   1.800 +    ///\param k shape parameter (<tt>k>0</tt>)
   1.801 +    ///\param x_min location parameter (<tt>x_min>0</tt>)
   1.802 +    ///
   1.803 +    double pareto(double k,double x_min)
   1.804 +    {
   1.805 +      return exponential(gamma(k,1.0/x_min));
   1.806 +    }  
   1.807 +      
   1.808 +    ///@}
   1.809 +    
   1.810 +    ///\name Two dimensional distributions
   1.811 +    ///
   1.812 +
   1.813 +    ///@{
   1.814 +    
   1.815 +    /// Uniform distribution on the full unit circle.
   1.816 +
   1.817 +    /// Uniform distribution on the full unit circle.
   1.818 +    ///
   1.819 +    dim2::Point<double> disc() 
   1.820 +    {
   1.821 +      double V1,V2;
   1.822 +      do {
   1.823 +	V1=2*real<double>()-1;
   1.824 +	V2=2*real<double>()-1;
   1.825 +	
   1.826 +      } while(V1*V1+V2*V2>=1);
   1.827 +      return dim2::Point<double>(V1,V2);
   1.828 +    }
   1.829 +    /// A kind of two dimensional Gauss distribution
   1.830 +
   1.831 +    /// This function provides a turning symmetric two-dimensional distribution.
   1.832 +    /// Both coordinates are of standard normal distribution, but they are not
   1.833 +    /// independent.
   1.834 +    ///
   1.835 +    /// \note The coordinates are the two random variables provided by
   1.836 +    /// the Box-Muller method.
   1.837 +    dim2::Point<double> gauss2()
   1.838 +    {
   1.839 +      double V1,V2,S;
   1.840 +      do {
   1.841 +	V1=2*real<double>()-1;
   1.842 +	V2=2*real<double>()-1;
   1.843 +	S=V1*V1+V2*V2;
   1.844 +      } while(S>=1);
   1.845 +      double W=std::sqrt(-2*std::log(S)/S);
   1.846 +      return dim2::Point<double>(W*V1,W*V2);
   1.847 +    }
   1.848 +    /// A kind of two dimensional exponential distribution
   1.849 +
   1.850 +    /// This function provides a turning symmetric two-dimensional distribution.
   1.851 +    /// The x-coordinate is of conditionally exponential distribution
   1.852 +    /// with the condition that x is positive and y=0. If x is negative and 
   1.853 +    /// y=0 then, -x is of exponential distribution. The same is true for the
   1.854 +    /// y-coordinate.
   1.855 +    dim2::Point<double> exponential2() 
   1.856 +    {
   1.857 +      double V1,V2,S;
   1.858 +      do {
   1.859 +	V1=2*real<double>()-1;
   1.860 +	V2=2*real<double>()-1;
   1.861 +	S=V1*V1+V2*V2;
   1.862 +      } while(S>=1);
   1.863 +      double W=-std::log(S)/S;
   1.864 +      return dim2::Point<double>(W*V1,W*V2);
   1.865 +    }
   1.866 +
   1.867 +    ///@}    
   1.868 +  };
   1.869 +
   1.870 +
   1.871 +  extern Random rnd;
   1.872 +
   1.873 +}
   1.874 +
   1.875 +#endif