Port random.h & Co. from svn -r3422 + some cleanups
authorAlpar Juttner <alpar@cs.elte.hu>
Fri, 21 Dec 2007 00:07:03 +0000
changeset 1099e499ca560b
parent 9 6ae00527c178
child 11 ea5945b2da9c
Port random.h & Co. from svn -r3422 + some cleanups
- gauss(double std_dev) has been remove for clarity
lemon/Makefile.am
lemon/random.cc
lemon/random.h
test/Makefile.am
test/random_test.cc
     1.1 --- a/lemon/Makefile.am	Thu Dec 20 23:46:50 2007 +0000
     1.2 +++ b/lemon/Makefile.am	Fri Dec 21 00:07:03 2007 +0000
     1.3 @@ -7,13 +7,16 @@
     1.4  lib_LTLIBRARIES += lemon/libemon.la
     1.5  
     1.6  lemon_libemon_la_SOURCES = \
     1.7 -        lemon/base.cc
     1.8 +        lemon/base.cc \
     1.9 +        lemon/random.cc
    1.10 +
    1.11  
    1.12  lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS)
    1.13  lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS)
    1.14  
    1.15  lemon_HEADERS += \
    1.16          lemon/dim2.h \
    1.17 +        lemon/random.h \
    1.18  	lemon/list_graph.h \
    1.19          lemon/tolerance.h
    1.20  
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/random.cc	Fri Dec 21 00:07:03 2007 +0000
     2.3 @@ -0,0 +1,29 @@
     2.4 +/* -*- C++ -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library
     2.7 + *
     2.8 + * Copyright (C) 2003-2007
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +///\file
    2.23 +///\brief Instantiation of the Random class.
    2.24 +
    2.25 +#include <lemon/random.h>
    2.26 +
    2.27 +namespace lemon {
    2.28 +  /// \brief Global random number generator instance
    2.29 +  ///
    2.30 +  /// A global mersenne twister random number generator instance
    2.31 +  Random rnd;
    2.32 +}
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/lemon/random.h	Fri Dec 21 00:07:03 2007 +0000
     3.3 @@ -0,0 +1,850 @@
     3.4 +/* -*- C++ -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library
     3.7 + *
     3.8 + * Copyright (C) 2003-2007
     3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    3.11 + *
    3.12 + * Permission to use, modify and distribute this software is granted
    3.13 + * provided that this copyright notice appears in all copies. For
    3.14 + * precise terms see the accompanying LICENSE file.
    3.15 + *
    3.16 + * This software is provided "AS IS" with no warranty of any kind,
    3.17 + * express or implied, and with no claim as to its suitability for any
    3.18 + * purpose.
    3.19 + *
    3.20 + */
    3.21 +
    3.22 +/*
    3.23 + * This file contains the reimplemented version of the Mersenne Twister
    3.24 + * Generator of Matsumoto and Nishimura.
    3.25 + *
    3.26 + * See the appropriate copyright notice below.
    3.27 + * 
    3.28 + * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
    3.29 + * All rights reserved.                          
    3.30 + *
    3.31 + * Redistribution and use in source and binary forms, with or without
    3.32 + * modification, are permitted provided that the following conditions
    3.33 + * are met:
    3.34 + *
    3.35 + * 1. Redistributions of source code must retain the above copyright
    3.36 + *    notice, this list of conditions and the following disclaimer.
    3.37 + *
    3.38 + * 2. Redistributions in binary form must reproduce the above copyright
    3.39 + *    notice, this list of conditions and the following disclaimer in the
    3.40 + *    documentation and/or other materials provided with the distribution.
    3.41 + *
    3.42 + * 3. The names of its contributors may not be used to endorse or promote 
    3.43 + *    products derived from this software without specific prior written 
    3.44 + *    permission.
    3.45 + *
    3.46 + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
    3.47 + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
    3.48 + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
    3.49 + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
    3.50 + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
    3.51 + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    3.52 + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
    3.53 + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
    3.54 + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
    3.55 + * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    3.56 + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
    3.57 + * OF THE POSSIBILITY OF SUCH DAMAGE.
    3.58 + *
    3.59 + *
    3.60 + * Any feedback is very welcome.
    3.61 + * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
    3.62 + * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
    3.63 + */
    3.64 +
    3.65 +#ifndef LEMON_RANDOM_H
    3.66 +#define LEMON_RANDOM_H
    3.67 +
    3.68 +#include <algorithm>
    3.69 +#include <iterator>
    3.70 +#include <vector>
    3.71 +
    3.72 +#include <ctime>
    3.73 +#include <cmath>
    3.74 +
    3.75 +#include <lemon/dim2.h>
    3.76 +///\ingroup misc
    3.77 +///\file
    3.78 +///\brief Mersenne Twister random number generator
    3.79 +///
    3.80 +///\author Balazs Dezso
    3.81 +
    3.82 +namespace lemon {
    3.83 +
    3.84 +  namespace _random_bits {
    3.85 +    
    3.86 +    template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
    3.87 +    struct RandomTraits {};
    3.88 +
    3.89 +    template <typename _Word>
    3.90 +    struct RandomTraits<_Word, 32> {
    3.91 +
    3.92 +      typedef _Word Word;
    3.93 +      static const int bits = 32;
    3.94 +
    3.95 +      static const int length = 624;
    3.96 +      static const int shift = 397;
    3.97 +      
    3.98 +      static const Word mul = 0x6c078965u;
    3.99 +      static const Word arrayInit = 0x012BD6AAu;
   3.100 +      static const Word arrayMul1 = 0x0019660Du;
   3.101 +      static const Word arrayMul2 = 0x5D588B65u;
   3.102 +
   3.103 +      static const Word mask = 0x9908B0DFu;
   3.104 +      static const Word loMask = (1u << 31) - 1;
   3.105 +      static const Word hiMask = ~loMask;
   3.106 +
   3.107 +
   3.108 +      static Word tempering(Word rnd) {
   3.109 +        rnd ^= (rnd >> 11);
   3.110 +        rnd ^= (rnd << 7) & 0x9D2C5680u;
   3.111 +        rnd ^= (rnd << 15) & 0xEFC60000u;
   3.112 +        rnd ^= (rnd >> 18);
   3.113 +        return rnd;
   3.114 +      }
   3.115 +
   3.116 +    };
   3.117 +
   3.118 +    template <typename _Word>
   3.119 +    struct RandomTraits<_Word, 64> {
   3.120 +
   3.121 +      typedef _Word Word;
   3.122 +      static const int bits = 64;
   3.123 +
   3.124 +      static const int length = 312;
   3.125 +      static const int shift = 156;
   3.126 +
   3.127 +      static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du);
   3.128 +      static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu);
   3.129 +      static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u);
   3.130 +      static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu);
   3.131 +
   3.132 +      static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u);
   3.133 +      static const Word loMask = (Word(1u) << 31) - 1;
   3.134 +      static const Word hiMask = ~loMask;
   3.135 +
   3.136 +      static Word tempering(Word rnd) {
   3.137 +        rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u));
   3.138 +        rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u));
   3.139 +        rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u));
   3.140 +        rnd ^= (rnd >> 43);
   3.141 +        return rnd;
   3.142 +      }
   3.143 +
   3.144 +    };
   3.145 +
   3.146 +    template <typename _Word>
   3.147 +    class RandomCore {
   3.148 +    public:
   3.149 +
   3.150 +      typedef _Word Word;
   3.151 +
   3.152 +    private:
   3.153 +
   3.154 +      static const int bits = RandomTraits<Word>::bits;
   3.155 +
   3.156 +      static const int length = RandomTraits<Word>::length;
   3.157 +      static const int shift = RandomTraits<Word>::shift;
   3.158 +
   3.159 +    public:
   3.160 +
   3.161 +      void initState() {
   3.162 +        static const Word seedArray[4] = {
   3.163 +          0x12345u, 0x23456u, 0x34567u, 0x45678u
   3.164 +        };
   3.165 +    
   3.166 +        initState(seedArray, seedArray + 4);
   3.167 +      }
   3.168 +
   3.169 +      void initState(Word seed) {
   3.170 +
   3.171 +        static const Word mul = RandomTraits<Word>::mul;
   3.172 +
   3.173 +        current = state; 
   3.174 +
   3.175 +        Word *curr = state + length - 1;
   3.176 +        curr[0] = seed; --curr;
   3.177 +        for (int i = 1; i < length; ++i) {
   3.178 +          curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
   3.179 +          --curr;
   3.180 +        }
   3.181 +      }
   3.182 +
   3.183 +      template <typename Iterator>
   3.184 +      void initState(Iterator begin, Iterator end) {
   3.185 +
   3.186 +        static const Word init = RandomTraits<Word>::arrayInit;
   3.187 +        static const Word mul1 = RandomTraits<Word>::arrayMul1;
   3.188 +        static const Word mul2 = RandomTraits<Word>::arrayMul2;
   3.189 +
   3.190 +
   3.191 +        Word *curr = state + length - 1; --curr;
   3.192 +        Iterator it = begin; int cnt = 0;
   3.193 +        int num;
   3.194 +
   3.195 +        initState(init);
   3.196 +
   3.197 +        num = length > end - begin ? length : end - begin;
   3.198 +        while (num--) {
   3.199 +          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 
   3.200 +            + *it + cnt;
   3.201 +          ++it; ++cnt;
   3.202 +          if (it == end) {
   3.203 +            it = begin; cnt = 0;
   3.204 +          }
   3.205 +          if (curr == state) {
   3.206 +            curr = state + length - 1; curr[0] = state[0];
   3.207 +          }
   3.208 +          --curr;
   3.209 +        }
   3.210 +
   3.211 +        num = length - 1; cnt = length - (curr - state) - 1;
   3.212 +        while (num--) {
   3.213 +          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
   3.214 +            - cnt;
   3.215 +          --curr; ++cnt;
   3.216 +          if (curr == state) {
   3.217 +            curr = state + length - 1; curr[0] = state[0]; --curr;
   3.218 +            cnt = 1;
   3.219 +          }
   3.220 +        }
   3.221 +        
   3.222 +        state[length - 1] = Word(1) << (bits - 1);
   3.223 +      }
   3.224 +      
   3.225 +      void copyState(const RandomCore& other) {
   3.226 +        std::copy(other.state, other.state + length, state);
   3.227 +        current = state + (other.current - other.state);
   3.228 +      }
   3.229 +
   3.230 +      Word operator()() {
   3.231 +        if (current == state) fillState();
   3.232 +        --current;
   3.233 +        Word rnd = *current;
   3.234 +        return RandomTraits<Word>::tempering(rnd);
   3.235 +      }
   3.236 +
   3.237 +    private:
   3.238 +
   3.239 +  
   3.240 +      void fillState() {
   3.241 +        static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
   3.242 +        static const Word loMask = RandomTraits<Word>::loMask;
   3.243 +        static const Word hiMask = RandomTraits<Word>::hiMask;
   3.244 +
   3.245 +        current = state + length; 
   3.246 +
   3.247 +        register Word *curr = state + length - 1;
   3.248 +        register long num;
   3.249 +      
   3.250 +        num = length - shift;
   3.251 +        while (num--) {
   3.252 +          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   3.253 +            curr[- shift] ^ mask[curr[-1] & 1ul];
   3.254 +          --curr;
   3.255 +        }
   3.256 +        num = shift - 1;
   3.257 +        while (num--) {
   3.258 +          curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
   3.259 +            curr[length - shift] ^ mask[curr[-1] & 1ul];
   3.260 +          --curr;
   3.261 +        }
   3.262 +        curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
   3.263 +          curr[length - shift] ^ mask[curr[length - 1] & 1ul];
   3.264 +
   3.265 +      }
   3.266 +
   3.267 +  
   3.268 +      Word *current;
   3.269 +      Word state[length];
   3.270 +      
   3.271 +    };
   3.272 +
   3.273 +
   3.274 +    template <typename Result, 
   3.275 +              int shift = (std::numeric_limits<Result>::digits + 1) / 2>
   3.276 +    struct Masker {
   3.277 +      static Result mask(const Result& result) {
   3.278 +        return Masker<Result, (shift + 1) / 2>::
   3.279 +          mask(static_cast<Result>(result | (result >> shift)));
   3.280 +      }
   3.281 +    };
   3.282 +    
   3.283 +    template <typename Result>
   3.284 +    struct Masker<Result, 1> {
   3.285 +      static Result mask(const Result& result) {
   3.286 +        return static_cast<Result>(result | (result >> 1));
   3.287 +      }
   3.288 +    };
   3.289 +
   3.290 +    template <typename Result, typename Word, 
   3.291 +              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   3.292 +              bool last = rest <= std::numeric_limits<Word>::digits>
   3.293 +    struct IntConversion {
   3.294 +      static const int bits = std::numeric_limits<Word>::digits;
   3.295 +    
   3.296 +      static Result convert(RandomCore<Word>& rnd) {
   3.297 +        return static_cast<Result>(rnd() >> (bits - rest)) << shift;
   3.298 +      }
   3.299 +      
   3.300 +    }; 
   3.301 +
   3.302 +    template <typename Result, typename Word, int rest, int shift> 
   3.303 +    struct IntConversion<Result, Word, rest, shift, false> {
   3.304 +      static const int bits = std::numeric_limits<Word>::digits;
   3.305 +
   3.306 +      static Result convert(RandomCore<Word>& rnd) {
   3.307 +        return (static_cast<Result>(rnd()) << shift) | 
   3.308 +          IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
   3.309 +      }
   3.310 +    };
   3.311 +
   3.312 +
   3.313 +    template <typename Result, typename Word,
   3.314 +              bool one_word = (std::numeric_limits<Word>::digits < 
   3.315 +			       std::numeric_limits<Result>::digits) >
   3.316 +    struct Mapping {
   3.317 +      static Result map(RandomCore<Word>& rnd, const Result& bound) {
   3.318 +        Word max = Word(bound - 1);
   3.319 +        Result mask = Masker<Result>::mask(bound - 1);
   3.320 +        Result num;
   3.321 +        do {
   3.322 +          num = IntConversion<Result, Word>::convert(rnd) & mask; 
   3.323 +        } while (num > max);
   3.324 +        return num;
   3.325 +      }
   3.326 +    };
   3.327 +
   3.328 +    template <typename Result, typename Word>
   3.329 +    struct Mapping<Result, Word, false> {
   3.330 +      static Result map(RandomCore<Word>& rnd, const Result& bound) {
   3.331 +        Word max = Word(bound - 1);
   3.332 +        Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
   3.333 +          ::mask(max);
   3.334 +        Word num;
   3.335 +        do {
   3.336 +          num = rnd() & mask;
   3.337 +        } while (num > max);
   3.338 +        return num;
   3.339 +      }
   3.340 +    };
   3.341 +
   3.342 +    template <typename Result, int exp, bool pos = (exp >= 0)>
   3.343 +    struct ShiftMultiplier {
   3.344 +      static const Result multiplier() {
   3.345 +        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   3.346 +        res *= res;
   3.347 +        if ((exp & 1) == 1) res *= static_cast<Result>(2.0);
   3.348 +        return res; 
   3.349 +      }
   3.350 +    };
   3.351 +
   3.352 +    template <typename Result, int exp>
   3.353 +    struct ShiftMultiplier<Result, exp, false> {
   3.354 +      static const Result multiplier() {
   3.355 +        Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
   3.356 +        res *= res;
   3.357 +        if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
   3.358 +        return res; 
   3.359 +      }
   3.360 +    };
   3.361 +
   3.362 +    template <typename Result>
   3.363 +    struct ShiftMultiplier<Result, 0, true> {
   3.364 +      static const Result multiplier() {
   3.365 +        return static_cast<Result>(1.0); 
   3.366 +      }
   3.367 +    };
   3.368 +
   3.369 +    template <typename Result>
   3.370 +    struct ShiftMultiplier<Result, -20, true> {
   3.371 +      static const Result multiplier() {
   3.372 +        return static_cast<Result>(1.0/1048576.0); 
   3.373 +      }
   3.374 +    };
   3.375 +    
   3.376 +    template <typename Result>
   3.377 +    struct ShiftMultiplier<Result, -32, true> {
   3.378 +      static const Result multiplier() {
   3.379 +        return static_cast<Result>(1.0/424967296.0); 
   3.380 +      }
   3.381 +    };
   3.382 +
   3.383 +    template <typename Result>
   3.384 +    struct ShiftMultiplier<Result, -53, true> {
   3.385 +      static const Result multiplier() {
   3.386 +        return static_cast<Result>(1.0/9007199254740992.0); 
   3.387 +      }
   3.388 +    };
   3.389 +
   3.390 +    template <typename Result>
   3.391 +    struct ShiftMultiplier<Result, -64, true> {
   3.392 +      static const Result multiplier() {
   3.393 +        return static_cast<Result>(1.0/18446744073709551616.0); 
   3.394 +      }
   3.395 +    };
   3.396 +
   3.397 +    template <typename Result, int exp>
   3.398 +    struct Shifting {
   3.399 +      static Result shift(const Result& result) {
   3.400 +        return result * ShiftMultiplier<Result, exp>::multiplier();
   3.401 +      }
   3.402 +    };
   3.403 +
   3.404 +    template <typename Result, typename Word,
   3.405 +              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
   3.406 +              bool last = rest <= std::numeric_limits<Word>::digits>
   3.407 +    struct RealConversion{ 
   3.408 +      static const int bits = std::numeric_limits<Word>::digits;
   3.409 +
   3.410 +      static Result convert(RandomCore<Word>& rnd) {
   3.411 +        return Shifting<Result, - shift - rest>::
   3.412 +          shift(static_cast<Result>(rnd() >> (bits - rest)));
   3.413 +      }
   3.414 +    };
   3.415 +
   3.416 +    template <typename Result, typename Word, int rest, int shift>
   3.417 +    struct RealConversion<Result, Word, rest, shift, false> { 
   3.418 +      static const int bits = std::numeric_limits<Word>::digits;
   3.419 +
   3.420 +      static Result convert(RandomCore<Word>& rnd) {
   3.421 +        return Shifting<Result, - shift - bits>::
   3.422 +          shift(static_cast<Result>(rnd())) +
   3.423 +          RealConversion<Result, Word, rest-bits, shift + bits>::
   3.424 +          convert(rnd);
   3.425 +      }
   3.426 +    };
   3.427 +
   3.428 +    template <typename Result, typename Word>
   3.429 +    struct Initializer {
   3.430 +
   3.431 +      template <typename Iterator>
   3.432 +      static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
   3.433 +        std::vector<Word> ws;
   3.434 +        for (Iterator it = begin; it != end; ++it) {
   3.435 +          ws.push_back(Word(*it));
   3.436 +        }
   3.437 +        rnd.initState(ws.begin(), ws.end());
   3.438 +      }
   3.439 +
   3.440 +      static void init(RandomCore<Word>& rnd, Result seed) {
   3.441 +        rnd.initState(seed);
   3.442 +      }
   3.443 +    };
   3.444 +
   3.445 +    template <typename Word>
   3.446 +    struct BoolConversion {
   3.447 +      static bool convert(RandomCore<Word>& rnd) {
   3.448 +        return (rnd() & 1) == 1;
   3.449 +      }
   3.450 +    };
   3.451 +
   3.452 +    template <typename Word>
   3.453 +    struct BoolProducer {
   3.454 +      Word buffer;
   3.455 +      int num;
   3.456 +      
   3.457 +      BoolProducer() : num(0) {}
   3.458 +
   3.459 +      bool convert(RandomCore<Word>& rnd) {
   3.460 +        if (num == 0) {
   3.461 +          buffer = rnd();
   3.462 +          num = RandomTraits<Word>::bits;
   3.463 +        }
   3.464 +        bool r = (buffer & 1);
   3.465 +        buffer >>= 1;
   3.466 +        --num;
   3.467 +        return r;
   3.468 +      }
   3.469 +    };
   3.470 +
   3.471 +  }
   3.472 +
   3.473 +  /// \ingroup misc
   3.474 +  ///
   3.475 +  /// \brief Mersenne Twister random number generator
   3.476 +  ///
   3.477 +  /// The Mersenne Twister is a twisted generalized feedback
   3.478 +  /// shift-register generator of Matsumoto and Nishimura. The period
   3.479 +  /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
   3.480 +  /// equi-distributed in 623 dimensions for 32-bit numbers. The time
   3.481 +  /// performance of this generator is comparable to the commonly used
   3.482 +  /// generators.
   3.483 +  ///
   3.484 +  /// This implementation is specialized for both 32-bit and 64-bit
   3.485 +  /// architectures. The generators differ sligthly in the
   3.486 +  /// initialization and generation phase so they produce two
   3.487 +  /// completly different sequences.
   3.488 +  ///
   3.489 +  /// The generator gives back random numbers of serveral types. To
   3.490 +  /// get a random number from a range of a floating point type you
   3.491 +  /// can use one form of the \c operator() or the \c real() member
   3.492 +  /// function. If you want to get random number from the {0, 1, ...,
   3.493 +  /// n-1} integer range use the \c operator[] or the \c integer()
   3.494 +  /// method. And to get random number from the whole range of an
   3.495 +  /// integer type you can use the argumentless \c integer() or \c
   3.496 +  /// uinteger() functions. After all you can get random bool with
   3.497 +  /// equal chance of true and false or given probability of true
   3.498 +  /// result with the \c boolean() member functions.
   3.499 +  ///
   3.500 +  ///\code
   3.501 +  /// // The commented code is identical to the other
   3.502 +  /// double a = rnd();                     // [0.0, 1.0)
   3.503 +  /// // double a = rnd.real();             // [0.0, 1.0)
   3.504 +  /// double b = rnd(100.0);                // [0.0, 100.0)
   3.505 +  /// // double b = rnd.real(100.0);        // [0.0, 100.0)
   3.506 +  /// double c = rnd(1.0, 2.0);             // [1.0, 2.0)
   3.507 +  /// // double c = rnd.real(1.0, 2.0);     // [1.0, 2.0)
   3.508 +  /// int d = rnd[100000];                  // 0..99999
   3.509 +  /// // int d = rnd.integer(100000);       // 0..99999
   3.510 +  /// int e = rnd[6] + 1;                   // 1..6
   3.511 +  /// // int e = rnd.integer(1, 1 + 6);     // 1..6
   3.512 +  /// int b = rnd.uinteger<int>();          // 0 .. 2^31 - 1
   3.513 +  /// int c = rnd.integer<int>();           // - 2^31 .. 2^31 - 1
   3.514 +  /// bool g = rnd.boolean();               // P(g = true) = 0.5
   3.515 +  /// bool h = rnd.boolean(0.8);            // P(h = true) = 0.8
   3.516 +  ///\endcode
   3.517 +  ///
   3.518 +  /// The lemon provides a global instance of the random number
   3.519 +  /// generator which name is \ref lemon::rnd "rnd". Usually it is a
   3.520 +  /// good programming convenience to use this global generator to get
   3.521 +  /// random numbers.
   3.522 +  ///
   3.523 +  /// \author Balazs Dezso
   3.524 +  class Random {
   3.525 +  private:
   3.526 +
   3.527 +    // architecture word
   3.528 +    typedef unsigned long Word;
   3.529 +    
   3.530 +    _random_bits::RandomCore<Word> core;
   3.531 +    _random_bits::BoolProducer<Word> bool_producer;
   3.532 +    
   3.533 +
   3.534 +  public:
   3.535 +
   3.536 +    /// \brief Constructor
   3.537 +    ///
   3.538 +    /// Constructor with constant seeding.
   3.539 +    Random() { core.initState(); }
   3.540 +
   3.541 +    /// \brief Constructor
   3.542 +    ///
   3.543 +    /// Constructor with seed. The current number type will be converted
   3.544 +    /// to the architecture word type.
   3.545 +    template <typename Number>
   3.546 +    Random(Number seed) { 
   3.547 +      _random_bits::Initializer<Number, Word>::init(core, seed);
   3.548 +    }
   3.549 +
   3.550 +    /// \brief Constructor
   3.551 +    ///
   3.552 +    /// Constructor with array seeding. The given range should contain
   3.553 +    /// any number type and the numbers will be converted to the
   3.554 +    /// architecture word type.
   3.555 +    template <typename Iterator>
   3.556 +    Random(Iterator begin, Iterator end) { 
   3.557 +      typedef typename std::iterator_traits<Iterator>::value_type Number;
   3.558 +      _random_bits::Initializer<Number, Word>::init(core, begin, end);
   3.559 +    }
   3.560 +
   3.561 +    /// \brief Copy constructor
   3.562 +    ///
   3.563 +    /// Copy constructor. The generated sequence will be identical to
   3.564 +    /// the other sequence. It can be used to save the current state
   3.565 +    /// of the generator and later use it to generate the same
   3.566 +    /// sequence.
   3.567 +    Random(const Random& other) {
   3.568 +      core.copyState(other.core);
   3.569 +    }
   3.570 +
   3.571 +    /// \brief Assign operator
   3.572 +    ///
   3.573 +    /// Assign operator. The generated sequence will be identical to
   3.574 +    /// the other sequence. It can be used to save the current state
   3.575 +    /// of the generator and later use it to generate the same
   3.576 +    /// sequence.
   3.577 +    Random& operator=(const Random& other) {
   3.578 +      if (&other != this) {
   3.579 +        core.copyState(other.core);
   3.580 +      }
   3.581 +      return *this;
   3.582 +    }
   3.583 +
   3.584 +    /// \brief Returns a random real number from the range [0, 1)
   3.585 +    ///
   3.586 +    /// It returns a random real number from the range [0, 1). The
   3.587 +    /// default Number type is double.
   3.588 +    template <typename Number>
   3.589 +    Number real() {
   3.590 +      return _random_bits::RealConversion<Number, Word>::convert(core);
   3.591 +    }
   3.592 +
   3.593 +    double real() {
   3.594 +      return real<double>();
   3.595 +    }
   3.596 +
   3.597 +    /// \brief Returns a random real number the range [0, b)
   3.598 +    ///
   3.599 +    /// It returns a random real number from the range [0, b).
   3.600 +    template <typename Number>
   3.601 +    Number real(Number b) { 
   3.602 +      return real<Number>() * b; 
   3.603 +    }
   3.604 +
   3.605 +    /// \brief Returns a random real number from the range [a, b)
   3.606 +    ///
   3.607 +    /// It returns a random real number from the range [a, b).
   3.608 +    template <typename Number>
   3.609 +    Number real(Number a, Number b) { 
   3.610 +      return real<Number>() * (b - a) + a; 
   3.611 +    }
   3.612 +
   3.613 +    /// \brief Returns a random real number from the range [0, 1)
   3.614 +    ///
   3.615 +    /// It returns a random double from the range [0, 1).
   3.616 +    double operator()() {
   3.617 +      return real<double>();
   3.618 +    }
   3.619 +
   3.620 +    /// \brief Returns a random real number from the range [0, b)
   3.621 +    ///
   3.622 +    /// It returns a random real number from the range [0, b).
   3.623 +    template <typename Number>
   3.624 +    Number operator()(Number b) { 
   3.625 +      return real<Number>() * b; 
   3.626 +    }
   3.627 +
   3.628 +    /// \brief Returns a random real number from the range [a, b)
   3.629 +    ///
   3.630 +    /// It returns a random real number from the range [a, b).
   3.631 +    template <typename Number>
   3.632 +    Number operator()(Number a, Number b) { 
   3.633 +      return real<Number>() * (b - a) + a; 
   3.634 +    }
   3.635 +
   3.636 +    /// \brief Returns a random integer from a range
   3.637 +    ///
   3.638 +    /// It returns a random integer from the range {0, 1, ..., b - 1}.
   3.639 +    template <typename Number>
   3.640 +    Number integer(Number b) {
   3.641 +      return _random_bits::Mapping<Number, Word>::map(core, b);
   3.642 +    }
   3.643 +
   3.644 +    /// \brief Returns a random integer from a range
   3.645 +    ///
   3.646 +    /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
   3.647 +    template <typename Number>
   3.648 +    Number integer(Number a, Number b) {
   3.649 +      return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
   3.650 +    }
   3.651 +
   3.652 +    /// \brief Returns a random integer from a range
   3.653 +    ///
   3.654 +    /// It returns a random integer from the range {0, 1, ..., b - 1}.
   3.655 +    template <typename Number>
   3.656 +    Number operator[](Number b) {
   3.657 +      return _random_bits::Mapping<Number, Word>::map(core, b);
   3.658 +    }
   3.659 +
   3.660 +    /// \brief Returns a random non-negative integer
   3.661 +    ///
   3.662 +    /// It returns a random non-negative integer uniformly from the
   3.663 +    /// whole range of the current \c Number type.  The default result
   3.664 +    /// type of this function is unsigned int.
   3.665 +    template <typename Number>
   3.666 +    Number uinteger() {
   3.667 +      return _random_bits::IntConversion<Number, Word>::convert(core);
   3.668 +    }
   3.669 +
   3.670 +    unsigned int uinteger() {
   3.671 +      return uinteger<unsigned int>();
   3.672 +    }
   3.673 +
   3.674 +    /// \brief Returns a random integer
   3.675 +    ///
   3.676 +    /// It returns a random integer uniformly from the whole range of
   3.677 +    /// the current \c Number type. The default result type of this
   3.678 +    /// function is int.
   3.679 +    template <typename Number>
   3.680 +    Number integer() {
   3.681 +      static const int nb = std::numeric_limits<Number>::digits + 
   3.682 +        (std::numeric_limits<Number>::is_signed ? 1 : 0);
   3.683 +      return _random_bits::IntConversion<Number, Word, nb>::convert(core);
   3.684 +    }
   3.685 +
   3.686 +    int integer() {
   3.687 +      return integer<int>();
   3.688 +    }
   3.689 +    
   3.690 +    /// \brief Returns a random bool
   3.691 +    ///
   3.692 +    /// It returns a random bool. The generator holds a buffer for
   3.693 +    /// random bits. Every time when it become empty the generator makes
   3.694 +    /// a new random word and fill the buffer up.
   3.695 +    bool boolean() {
   3.696 +      return bool_producer.convert(core);
   3.697 +    }
   3.698 +
   3.699 +    ///\name Nonuniform distributions
   3.700 +    ///
   3.701 +    
   3.702 +    ///@{
   3.703 +    
   3.704 +    /// \brief Returns a random bool
   3.705 +    ///
   3.706 +    /// It returns a random bool with given probability of true result
   3.707 +    bool boolean(double p) {
   3.708 +      return operator()() < p;
   3.709 +    }
   3.710 +
   3.711 +    /// Standard Gauss distribution
   3.712 +
   3.713 +    /// Standard Gauss distribution.
   3.714 +    /// \note The Cartesian form of the Box-Muller
   3.715 +    /// transformation is used to generate a random normal distribution.
   3.716 +    /// \todo Consider using the "ziggurat" method instead.
   3.717 +    double gauss() 
   3.718 +    {
   3.719 +      double V1,V2,S;
   3.720 +      do {
   3.721 +	V1=2*real<double>()-1;
   3.722 +	V2=2*real<double>()-1;
   3.723 +	S=V1*V1+V2*V2;
   3.724 +      } while(S>=1);
   3.725 +      return std::sqrt(-2*std::log(S)/S)*V1;
   3.726 +    }
   3.727 +    /// Gauss distribution with given mean and standard deviation
   3.728 +
   3.729 +    /// \sa gauss()
   3.730 +    ///
   3.731 +    double gauss(double mean,double std_dev)
   3.732 +    {
   3.733 +      return gauss()*std_dev+mean;
   3.734 +    }
   3.735 +
   3.736 +    /// Exponential distribution with given mean
   3.737 +
   3.738 +    /// This function generates an exponential distribution random number
   3.739 +    /// with mean <tt>1/lambda</tt>.
   3.740 +    ///
   3.741 +    double exponential(double lambda=1.0)
   3.742 +    {
   3.743 +      return -std::log(real<double>())/lambda;
   3.744 +    }
   3.745 +
   3.746 +    /// Gamma distribution with given integer shape
   3.747 +
   3.748 +    /// This function generates a gamma distribution random number.
   3.749 +    /// 
   3.750 +    ///\param k shape parameter (<tt>k>0</tt> integer)
   3.751 +    double gamma(int k) 
   3.752 +    {
   3.753 +      double s = 0;
   3.754 +      for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
   3.755 +      return s;
   3.756 +    }
   3.757 +    
   3.758 +    /// Gamma distribution with given shape and scale parameter
   3.759 +
   3.760 +    /// This function generates a gamma distribution random number.
   3.761 +    /// 
   3.762 +    ///\param k shape parameter (<tt>k>0</tt>)
   3.763 +    ///\param theta scale parameter
   3.764 +    ///
   3.765 +    double gamma(double k,double theta=1.0)
   3.766 +    {
   3.767 +      double xi,nu;
   3.768 +      const double delta = k-std::floor(k);
   3.769 +      const double v0=M_E/(M_E-delta);
   3.770 +      do {
   3.771 +	double V0=1.0-real<double>();
   3.772 +	double V1=1.0-real<double>();
   3.773 +	double V2=1.0-real<double>();
   3.774 +	if(V2<=v0) 
   3.775 +	  {
   3.776 +	    xi=std::pow(V1,1.0/delta);
   3.777 +	    nu=V0*std::pow(xi,delta-1.0);
   3.778 +	  }
   3.779 +	else 
   3.780 +	  {
   3.781 +	    xi=1.0-std::log(V1);
   3.782 +	    nu=V0*std::exp(-xi);
   3.783 +	  }
   3.784 +      } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
   3.785 +      return theta*(xi-gamma(int(std::floor(k))));
   3.786 +    }
   3.787 +    
   3.788 +      
   3.789 +    ///@}
   3.790 +    
   3.791 +    ///\name Two dimensional distributions
   3.792 +    ///
   3.793 +
   3.794 +    ///@{
   3.795 +    
   3.796 +    /// Uniform distribution on the full unit circle.
   3.797 +    dim2::Point<double> disc() 
   3.798 +    {
   3.799 +      double V1,V2;
   3.800 +      do {
   3.801 +	V1=2*real<double>()-1;
   3.802 +	V2=2*real<double>()-1;
   3.803 +	
   3.804 +      } while(V1*V1+V2*V2>=1);
   3.805 +      return dim2::Point<double>(V1,V2);
   3.806 +    }
   3.807 +    /// A kind of two dimensional Gauss distribution
   3.808 +
   3.809 +    /// This function provides a turning symmetric two-dimensional distribution.
   3.810 +    /// Both coordinates are of standard normal distribution, but they are not
   3.811 +    /// independent.
   3.812 +    ///
   3.813 +    /// \note The coordinates are the two random variables provided by
   3.814 +    /// the Box-Muller method.
   3.815 +    dim2::Point<double> gauss2()
   3.816 +    {
   3.817 +      double V1,V2,S;
   3.818 +      do {
   3.819 +	V1=2*real<double>()-1;
   3.820 +	V2=2*real<double>()-1;
   3.821 +	S=V1*V1+V2*V2;
   3.822 +      } while(S>=1);
   3.823 +      double W=std::sqrt(-2*std::log(S)/S);
   3.824 +      return dim2::Point<double>(W*V1,W*V2);
   3.825 +    }
   3.826 +    /// A kind of two dimensional exponential distribution
   3.827 +
   3.828 +    /// This function provides a turning symmetric two-dimensional distribution.
   3.829 +    /// The x-coordinate is of conditionally exponential distribution
   3.830 +    /// with the condition that x is positive and y=0. If x is negative and 
   3.831 +    /// y=0 then, -x is of exponential distribution. The same is true for the
   3.832 +    /// y-coordinate.
   3.833 +    dim2::Point<double> exponential2() 
   3.834 +    {
   3.835 +      double V1,V2,S;
   3.836 +      do {
   3.837 +	V1=2*real<double>()-1;
   3.838 +	V2=2*real<double>()-1;
   3.839 +	S=V1*V1+V2*V2;
   3.840 +      } while(S>=1);
   3.841 +      double W=-std::log(S)/S;
   3.842 +      return dim2::Point<double>(W*V1,W*V2);
   3.843 +    }
   3.844 +
   3.845 +    ///@}    
   3.846 +  };
   3.847 +
   3.848 +
   3.849 +  extern Random rnd;
   3.850 +
   3.851 +}
   3.852 +
   3.853 +#endif
     4.1 --- a/test/Makefile.am	Thu Dec 20 23:46:50 2007 +0000
     4.2 +++ b/test/Makefile.am	Fri Dec 21 00:07:03 2007 +0000
     4.3 @@ -3,15 +3,17 @@
     4.4  
     4.5  noinst_HEADERS += \
     4.6          test/test_tools.h
     4.7 - 
     4.8 +
     4.9  check_PROGRAMS += \
    4.10          test/dim_test \
    4.11 +        test/random_test \
    4.12          test/test_tools_fail \
    4.13          test/test_tools_pass
    4.14 - 
    4.15 +
    4.16  TESTS += $(check_PROGRAMS)
    4.17  XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
    4.18  
    4.19  test_dim_test_SOURCES = test/dim_test.cc
    4.20 +test_random_test_SOURCES = test/random_test.cc
    4.21  test_test_tools_fail_SOURCES = test/test_tools_fail.cc
    4.22  test_test_tools_pass_SOURCES = test/test_tools_pass.cc
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/test/random_test.cc	Fri Dec 21 00:07:03 2007 +0000
     5.3 @@ -0,0 +1,36 @@
     5.4 +/* -*- C++ -*-
     5.5 + *
     5.6 + * This file is a part of LEMON, a generic C++ optimization library
     5.7 + *
     5.8 + * Copyright (C) 2003-2007
     5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    5.11 + *
    5.12 + * Permission to use, modify and distribute this software is granted
    5.13 + * provided that this copyright notice appears in all copies. For
    5.14 + * precise terms see the accompanying LICENSE file.
    5.15 + *
    5.16 + * This software is provided "AS IS" with no warranty of any kind,
    5.17 + * express or implied, and with no claim as to its suitability for any
    5.18 + * purpose.
    5.19 + *
    5.20 + */
    5.21 +
    5.22 +#include <lemon/random.h>
    5.23 +#include "test_tools.h"
    5.24 +
    5.25 +///\file \brief Test cases for random.h
    5.26 +///
    5.27 +///\todo To be extended
    5.28 +///
    5.29 +
    5.30 +int main()
    5.31 +{
    5.32 +  double a=lemon::rnd();
    5.33 +  check(a<1.0&&a>0.0,"This should be in [0,1)");
    5.34 +  a=lemon::rnd.gauss();
    5.35 +  a=lemon::rnd.gamma(3.45,0);
    5.36 +  a=lemon::rnd.gamma(4);
    5.37 +  //Does gamma work with integer k?
    5.38 +  a=lemon::rnd.gamma(4.0,0);
    5.39 +}