1.1 --- a/.hgignore Thu Jan 03 11:12:27 2008 +0100
1.2 +++ b/.hgignore Thu Jan 03 11:13:29 2008 +0100
1.3 @@ -5,6 +5,9 @@
1.4 *~
1.5 *.o
1.6 .#.*
1.7 +*.log
1.8 +*.lo
1.9 +*.tar.*
1.10 Makefile.in
1.11 aclocal.m4
1.12 config.h.in
1.13 @@ -19,11 +22,13 @@
1.14 lemon/libemon.la
1.15 lemon/stamp-h2
1.16 doc/Doxyfile
1.17 -lemon/.dirstamp
1.18 -lemon/.libs/*
1.19 +.dirstamp
1.20 +.libs/*
1.21 +.deps/*
1.22
1.23 syntax: regexp
1.24 -html/.*
1.25 -autom4te.cache/.*
1.26 -build-aux/.*
1.27 -objs.*/.*
1.28 +^doc/html/.*
1.29 +^autom4te.cache/.*
1.30 +^build-aux/.*
1.31 +^objs.*/.*
1.32 +^test/[a-z_]*$
1.33 \ No newline at end of file
2.1 --- a/lemon/Makefile.am Thu Jan 03 11:12:27 2008 +0100
2.2 +++ b/lemon/Makefile.am Thu Jan 03 11:13:29 2008 +0100
2.3 @@ -7,13 +7,16 @@
2.4 lib_LTLIBRARIES += lemon/libemon.la
2.5
2.6 lemon_libemon_la_SOURCES = \
2.7 - lemon/base.cc
2.8 + lemon/base.cc \
2.9 + lemon/random.cc
2.10 +
2.11
2.12 lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS)
2.13 lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS)
2.14
2.15 lemon_HEADERS += \
2.16 lemon/dim2.h \
2.17 + lemon/random.h \
2.18 lemon/list_graph.h \
2.19 lemon/tolerance.h
2.20
3.1 --- a/lemon/bits/invalid.h Thu Jan 03 11:12:27 2008 +0100
3.2 +++ b/lemon/bits/invalid.h Thu Jan 03 11:13:29 2008 +0100
3.3 @@ -24,7 +24,7 @@
3.4
3.5 namespace lemon {
3.6
3.7 - /// \brief Dummy type to make it easier to make invalid iterators.
3.8 + /// \brief Dummy type to make it easier to create invalid iterators.
3.9 ///
3.10 /// See \ref INVALID for the usage.
3.11 struct Invalid {
3.12 @@ -34,8 +34,8 @@
3.13 bool operator< (Invalid) { return false; }
3.14 };
3.15
3.16 - /// Invalid iterators.
3.17 -
3.18 + /// \brief Invalid iterators.
3.19 + ///
3.20 /// \ref Invalid is a global type that converts to each iterator
3.21 /// in such a way that the value of the target iterator will be invalid.
3.22
4.1 --- a/lemon/dim2.h Thu Jan 03 11:12:27 2008 +0100
4.2 +++ b/lemon/dim2.h Thu Jan 03 11:13:29 2008 +0100
4.3 @@ -34,9 +34,6 @@
4.4 /// can be used to determine
4.5 /// the rectangular bounding box of a set of
4.6 /// \ref lemon::dim2::Point "dim2::Point"'s.
4.7 -///
4.8 -///\author Attila Bernath
4.9 -
4.10
4.11 namespace lemon {
4.12
4.13 @@ -62,9 +59,9 @@
4.14
4.15 typedef T Value;
4.16
4.17 - ///First co-ordinate
4.18 + ///First coordinate
4.19 T x;
4.20 - ///Second co-ordinate
4.21 + ///Second coordinate
4.22 T y;
4.23
4.24 ///Default constructor
4.25 @@ -75,8 +72,8 @@
4.26
4.27 ///The dimension of the vector.
4.28
4.29 - ///This class give back always 2.
4.30 - ///
4.31 + ///The dimension of the vector.
4.32 + ///This function always returns 2.
4.33 int size() const { return 2; }
4.34
4.35 ///Subscripting operator
4.36 @@ -138,7 +135,7 @@
4.37 return b+=u;
4.38 }
4.39
4.40 - ///Return the neg of the vectors
4.41 + ///Return the negative of the vector
4.42 Point<T> operator-() const {
4.43 Point<T> b=*this;
4.44 b.x=-b.x; b.y=-b.y;
4.45 @@ -175,9 +172,9 @@
4.46
4.47 };
4.48
4.49 - ///Return an Point
4.50 + ///Return a Point
4.51
4.52 - ///Return an Point
4.53 + ///Return a Point.
4.54 ///\relates Point
4.55 template <typename T>
4.56 inline Point<T> makePoint(const T& x, const T& y) {
4.57 @@ -186,7 +183,7 @@
4.58
4.59 ///Return a vector multiplied by a scalar
4.60
4.61 - ///Return a vector multiplied by a scalar
4.62 + ///Return a vector multiplied by a scalar.
4.63 ///\relates Point
4.64 template<typename T> Point<T> operator*(const T &u,const Point<T> &x) {
4.65 return x*u;
4.66 @@ -194,7 +191,7 @@
4.67
4.68 ///Read a plainvector from a stream
4.69
4.70 - ///Read a plainvector from a stream
4.71 + ///Read a plainvector from a stream.
4.72 ///\relates Point
4.73 ///
4.74 template<typename T>
4.75 @@ -222,7 +219,7 @@
4.76
4.77 ///Write a plainvector to a stream
4.78
4.79 - ///Write a plainvector to a stream
4.80 + ///Write a plainvector to a stream.
4.81 ///\relates Point
4.82 ///
4.83 template<typename T>
4.84 @@ -234,7 +231,7 @@
4.85
4.86 ///Rotate by 90 degrees
4.87
4.88 - ///Returns its parameter rotated by 90 degrees in positive direction.
4.89 + ///Returns the parameter rotated by 90 degrees in positive direction.
4.90 ///\relates Point
4.91 ///
4.92 template<typename T>
4.93 @@ -245,7 +242,7 @@
4.94
4.95 ///Rotate by 180 degrees
4.96
4.97 - ///Returns its parameter rotated by 180 degrees.
4.98 + ///Returns the parameter rotated by 180 degrees.
4.99 ///\relates Point
4.100 ///
4.101 template<typename T>
4.102 @@ -256,7 +253,7 @@
4.103
4.104 ///Rotate by 270 degrees
4.105
4.106 - ///Returns its parameter rotated by 90 degrees in negative direction.
4.107 + ///Returns the parameter rotated by 90 degrees in negative direction.
4.108 ///\relates Point
4.109 ///
4.110 template<typename T>
4.111 @@ -271,7 +268,6 @@
4.112
4.113 /// A class to calculate or store the bounding box of plainvectors.
4.114 ///
4.115 - ///\author Attila Bernath
4.116 template<typename T>
4.117 class BoundingBox {
4.118 Point<T> bottom_left, top_right;
4.119 @@ -286,9 +282,11 @@
4.120
4.121 ///Construct an instance from two points
4.122
4.123 - ///Construct an instance from two points
4.124 - ///\warning The coordinates of the bottom-left corner must be no more
4.125 - ///than those of the top-right one
4.126 + ///Construct an instance from two points.
4.127 + ///\param a The bottom left corner.
4.128 + ///\param b The top right corner.
4.129 + ///\warning The coordinates of the bottom left corner must be no more
4.130 + ///than those of the top right one.
4.131 BoundingBox(Point<T> a,Point<T> b)
4.132 {
4.133 bottom_left=a;
4.134 @@ -298,9 +296,13 @@
4.135
4.136 ///Construct an instance from four numbers
4.137
4.138 - ///Construct an instance from four numbers
4.139 - ///\warning The coordinates of the bottom-left corner must be no more
4.140 - ///than those of the top-right one
4.141 + ///Construct an instance from four numbers.
4.142 + ///\param l The left side of the box.
4.143 + ///\param b The bottom of the box.
4.144 + ///\param r The right side of the box.
4.145 + ///\param t The top of the box.
4.146 + ///\warning The left side must be no more than the right side and
4.147 + ///bottom must be no more than the top.
4.148 BoundingBox(T l,T b,T r,T t)
4.149 {
4.150 bottom_left=Point<T>(l,b);
4.151 @@ -308,7 +310,12 @@
4.152 _empty = false;
4.153 }
4.154
4.155 - ///Were any points added?
4.156 + ///Return \c true if the bounding box is empty.
4.157 +
4.158 + ///Return \c true if the bounding box is empty (i.e. return \c false
4.159 + ///if at least one point was added to the box or the coordinates of
4.160 + ///the box were set).
4.161 + ///The coordinates of an empty bounding box are not defined.
4.162 bool empty() const {
4.163 return _empty;
4.164 }
4.165 @@ -329,7 +336,7 @@
4.166 ///Set the bottom left corner
4.167
4.168 ///Set the bottom left corner.
4.169 - ///It should only bee used for non-empty box.
4.170 + ///It should only be used for non-empty box.
4.171 void bottomLeft(Point<T> p) {
4.172 bottom_left = p;
4.173 }
4.174 @@ -345,7 +352,7 @@
4.175 ///Set the top right corner
4.176
4.177 ///Set the top right corner.
4.178 - ///It should only bee used for non-empty box.
4.179 + ///It should only be used for non-empty box.
4.180 void topRight(Point<T> p) {
4.181 top_right = p;
4.182 }
4.183 @@ -361,7 +368,7 @@
4.184 ///Set the bottom right corner
4.185
4.186 ///Set the bottom right corner.
4.187 - ///It should only bee used for non-empty box.
4.188 + ///It should only be used for non-empty box.
4.189 void bottomRight(Point<T> p) {
4.190 top_right.x = p.x;
4.191 bottom_left.y = p.y;
4.192 @@ -378,7 +385,7 @@
4.193 ///Set the top left corner
4.194
4.195 ///Set the top left corner.
4.196 - ///It should only bee used for non-empty box.
4.197 + ///It should only be used for non-empty box.
4.198 void topLeft(Point<T> p) {
4.199 top_right.y = p.y;
4.200 bottom_left.x = p.x;
4.201 @@ -395,7 +402,7 @@
4.202 ///Set the bottom of the box
4.203
4.204 ///Set the bottom of the box.
4.205 - ///It should only bee used for non-empty box.
4.206 + ///It should only be used for non-empty box.
4.207 void bottom(T t) {
4.208 bottom_left.y = t;
4.209 }
4.210 @@ -411,7 +418,7 @@
4.211 ///Set the top of the box
4.212
4.213 ///Set the top of the box.
4.214 - ///It should only bee used for non-empty box.
4.215 + ///It should only be used for non-empty box.
4.216 void top(T t) {
4.217 top_right.y = t;
4.218 }
4.219 @@ -427,7 +434,7 @@
4.220 ///Set the left side of the box
4.221
4.222 ///Set the left side of the box.
4.223 - ///It should only bee used for non-empty box
4.224 + ///It should only be used for non-empty box.
4.225 void left(T t) {
4.226 bottom_left.x = t;
4.227 }
4.228 @@ -443,7 +450,7 @@
4.229 ///Set the right side of the box
4.230
4.231 ///Set the right side of the box.
4.232 - ///It should only bee used for non-empty box
4.233 + ///It should only be used for non-empty box.
4.234 void right(T t) {
4.235 top_right.x = t;
4.236 }
4.237 @@ -465,7 +472,7 @@
4.238 }
4.239
4.240 ///Checks whether a point is inside a bounding box
4.241 - bool inside(const Point<T>& u){
4.242 + bool inside(const Point<T>& u) const {
4.243 if (_empty)
4.244 return false;
4.245 else{
4.246 @@ -475,6 +482,9 @@
4.247 }
4.248
4.249 ///Increments a bounding box with a point
4.250 +
4.251 + ///Increments a bounding box with a point.
4.252 + ///
4.253 BoundingBox& add(const Point<T>& u){
4.254 if (_empty){
4.255 bottom_left=top_right=u;
4.256 @@ -489,7 +499,10 @@
4.257 return *this;
4.258 }
4.259
4.260 - ///Increments a bounding to contain another bounding box
4.261 + ///Increments a bounding box to contain another bounding box
4.262 +
4.263 + ///Increments a bounding box to contain another bounding box.
4.264 + ///
4.265 BoundingBox& add(const BoundingBox &u){
4.266 if ( !u.empty() ){
4.267 this->add(u.bottomLeft());
4.268 @@ -499,24 +512,31 @@
4.269 }
4.270
4.271 ///Intersection of two bounding boxes
4.272 - BoundingBox operator &(const BoundingBox& u){
4.273 +
4.274 + ///Intersection of two bounding boxes.
4.275 + ///
4.276 + BoundingBox operator&(const BoundingBox& u) const {
4.277 BoundingBox b;
4.278 - b.bottom_left.x=std::max(this->bottom_left.x,u.bottom_left.x);
4.279 - b.bottom_left.y=std::max(this->bottom_left.y,u.bottom_left.y);
4.280 - b.top_right.x=std::min(this->top_right.x,u.top_right.x);
4.281 - b.top_right.y=std::min(this->top_right.y,u.top_right.y);
4.282 - b._empty = this->_empty || u._empty ||
4.283 - b.bottom_left.x>top_right.x && b.bottom_left.y>top_right.y;
4.284 + if (this->_empty || u._empty) {
4.285 + b._empty = true;
4.286 + } else {
4.287 + b.bottom_left.x = std::max(this->bottom_left.x,u.bottom_left.x);
4.288 + b.bottom_left.y = std::max(this->bottom_left.y,u.bottom_left.y);
4.289 + b.top_right.x = std::min(this->top_right.x,u.top_right.x);
4.290 + b.top_right.y = std::min(this->top_right.y,u.top_right.y);
4.291 + b._empty = b.bottom_left.x > b.top_right.x ||
4.292 + b.bottom_left.y > b.top_right.y;
4.293 + }
4.294 return b;
4.295 }
4.296
4.297 };//class Boundingbox
4.298
4.299
4.300 - ///Map of x-coordinates of a dim2::Point<>-map
4.301 + ///Map of x-coordinates of a \ref Point "Point"-map
4.302
4.303 ///\ingroup maps
4.304 - ///Map of x-coordinates of a dim2::Point<>-map
4.305 + ///Map of x-coordinates of a \ref Point "Point"-map.
4.306 ///
4.307 template<class M>
4.308 class XMap
4.309 @@ -570,7 +590,7 @@
4.310
4.311 ///Returns a \ref ConstXMap class
4.312
4.313 - ///This function just returns an \ref ConstXMap class.
4.314 + ///This function just returns a \ref ConstXMap class.
4.315 ///
4.316 ///\ingroup maps
4.317 ///\relates ConstXMap
4.318 @@ -580,10 +600,10 @@
4.319 return ConstXMap<M>(m);
4.320 }
4.321
4.322 - ///Map of y-coordinates of a dim2::Point<>-map
4.323 + ///Map of y-coordinates of a \ref Point "Point"-map
4.324
4.325 ///\ingroup maps
4.326 - ///Map of y-coordinates of a dim2::Point<>-map
4.327 + ///Map of y-coordinates of a \ref Point "Point"-map.
4.328 ///
4.329 template<class M>
4.330 class YMap
4.331 @@ -599,9 +619,9 @@
4.332 void set(Key k,Value v) {_map.set(k,typename M::Value(_map[k].x,v));}
4.333 };
4.334
4.335 - ///Returns an \ref YMap class
4.336 + ///Returns a \ref YMap class
4.337
4.338 - ///This function just returns an \ref YMap class.
4.339 + ///This function just returns a \ref YMap class.
4.340 ///
4.341 ///\ingroup maps
4.342 ///\relates YMap
4.343 @@ -637,7 +657,7 @@
4.344
4.345 ///Returns a \ref ConstYMap class
4.346
4.347 - ///This function just returns an \ref ConstYMap class.
4.348 + ///This function just returns a \ref ConstYMap class.
4.349 ///
4.350 ///\ingroup maps
4.351 ///\relates ConstYMap
4.352 @@ -649,10 +669,10 @@
4.353
4.354
4.355 ///\brief Map of the \ref Point::normSquare() "normSquare()"
4.356 - ///of an \ref Point "Point"-map
4.357 + ///of a \ref Point "Point"-map
4.358 ///
4.359 ///Map of the \ref Point::normSquare() "normSquare()"
4.360 - ///of an \ref Point "Point"-map
4.361 + ///of a \ref Point "Point"-map.
4.362 ///\ingroup maps
4.363 ///
4.364 template<class M>
4.365 @@ -670,7 +690,7 @@
4.366
4.367 ///Returns a \ref NormSquareMap class
4.368
4.369 - ///This function just returns an \ref NormSquareMap class.
4.370 + ///This function just returns a \ref NormSquareMap class.
4.371 ///
4.372 ///\ingroup maps
4.373 ///\relates NormSquareMap
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/lemon/random.cc Thu Jan 03 11:13:29 2008 +0100
5.3 @@ -0,0 +1,29 @@
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 +///\file
5.23 +///\brief Instantiation of the Random class.
5.24 +
5.25 +#include <lemon/random.h>
5.26 +
5.27 +namespace lemon {
5.28 + /// \brief Global random number generator instance
5.29 + ///
5.30 + /// A global Mersenne Twister random number generator instance.
5.31 + Random rnd;
5.32 +}
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/lemon/random.h Thu Jan 03 11:13:29 2008 +0100
6.3 @@ -0,0 +1,872 @@
6.4 +/* -*- C++ -*-
6.5 + *
6.6 + * This file is a part of LEMON, a generic C++ optimization library
6.7 + *
6.8 + * Copyright (C) 2003-2007
6.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
6.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
6.11 + *
6.12 + * Permission to use, modify and distribute this software is granted
6.13 + * provided that this copyright notice appears in all copies. For
6.14 + * precise terms see the accompanying LICENSE file.
6.15 + *
6.16 + * This software is provided "AS IS" with no warranty of any kind,
6.17 + * express or implied, and with no claim as to its suitability for any
6.18 + * purpose.
6.19 + *
6.20 + */
6.21 +
6.22 +/*
6.23 + * This file contains the reimplemented version of the Mersenne Twister
6.24 + * Generator of Matsumoto and Nishimura.
6.25 + *
6.26 + * See the appropriate copyright notice below.
6.27 + *
6.28 + * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
6.29 + * All rights reserved.
6.30 + *
6.31 + * Redistribution and use in source and binary forms, with or without
6.32 + * modification, are permitted provided that the following conditions
6.33 + * are met:
6.34 + *
6.35 + * 1. Redistributions of source code must retain the above copyright
6.36 + * notice, this list of conditions and the following disclaimer.
6.37 + *
6.38 + * 2. Redistributions in binary form must reproduce the above copyright
6.39 + * notice, this list of conditions and the following disclaimer in the
6.40 + * documentation and/or other materials provided with the distribution.
6.41 + *
6.42 + * 3. The names of its contributors may not be used to endorse or promote
6.43 + * products derived from this software without specific prior written
6.44 + * permission.
6.45 + *
6.46 + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
6.47 + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
6.48 + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
6.49 + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
6.50 + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
6.51 + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
6.52 + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
6.53 + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
6.54 + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
6.55 + * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
6.56 + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
6.57 + * OF THE POSSIBILITY OF SUCH DAMAGE.
6.58 + *
6.59 + *
6.60 + * Any feedback is very welcome.
6.61 + * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
6.62 + * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
6.63 + */
6.64 +
6.65 +#ifndef LEMON_RANDOM_H
6.66 +#define LEMON_RANDOM_H
6.67 +
6.68 +#include <algorithm>
6.69 +#include <iterator>
6.70 +#include <vector>
6.71 +
6.72 +#include <ctime>
6.73 +#include <cmath>
6.74 +
6.75 +#include <lemon/dim2.h>
6.76 +///\ingroup misc
6.77 +///\file
6.78 +///\brief Mersenne Twister random number generator
6.79 +
6.80 +namespace lemon {
6.81 +
6.82 + namespace _random_bits {
6.83 +
6.84 + template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
6.85 + struct RandomTraits {};
6.86 +
6.87 + template <typename _Word>
6.88 + struct RandomTraits<_Word, 32> {
6.89 +
6.90 + typedef _Word Word;
6.91 + static const int bits = 32;
6.92 +
6.93 + static const int length = 624;
6.94 + static const int shift = 397;
6.95 +
6.96 + static const Word mul = 0x6c078965u;
6.97 + static const Word arrayInit = 0x012BD6AAu;
6.98 + static const Word arrayMul1 = 0x0019660Du;
6.99 + static const Word arrayMul2 = 0x5D588B65u;
6.100 +
6.101 + static const Word mask = 0x9908B0DFu;
6.102 + static const Word loMask = (1u << 31) - 1;
6.103 + static const Word hiMask = ~loMask;
6.104 +
6.105 +
6.106 + static Word tempering(Word rnd) {
6.107 + rnd ^= (rnd >> 11);
6.108 + rnd ^= (rnd << 7) & 0x9D2C5680u;
6.109 + rnd ^= (rnd << 15) & 0xEFC60000u;
6.110 + rnd ^= (rnd >> 18);
6.111 + return rnd;
6.112 + }
6.113 +
6.114 + };
6.115 +
6.116 + template <typename _Word>
6.117 + struct RandomTraits<_Word, 64> {
6.118 +
6.119 + typedef _Word Word;
6.120 + static const int bits = 64;
6.121 +
6.122 + static const int length = 312;
6.123 + static const int shift = 156;
6.124 +
6.125 + static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du);
6.126 + static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu);
6.127 + static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u);
6.128 + static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu);
6.129 +
6.130 + static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u);
6.131 + static const Word loMask = (Word(1u) << 31) - 1;
6.132 + static const Word hiMask = ~loMask;
6.133 +
6.134 + static Word tempering(Word rnd) {
6.135 + rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u));
6.136 + rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u));
6.137 + rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u));
6.138 + rnd ^= (rnd >> 43);
6.139 + return rnd;
6.140 + }
6.141 +
6.142 + };
6.143 +
6.144 + template <typename _Word>
6.145 + class RandomCore {
6.146 + public:
6.147 +
6.148 + typedef _Word Word;
6.149 +
6.150 + private:
6.151 +
6.152 + static const int bits = RandomTraits<Word>::bits;
6.153 +
6.154 + static const int length = RandomTraits<Word>::length;
6.155 + static const int shift = RandomTraits<Word>::shift;
6.156 +
6.157 + public:
6.158 +
6.159 + void initState() {
6.160 + static const Word seedArray[4] = {
6.161 + 0x12345u, 0x23456u, 0x34567u, 0x45678u
6.162 + };
6.163 +
6.164 + initState(seedArray, seedArray + 4);
6.165 + }
6.166 +
6.167 + void initState(Word seed) {
6.168 +
6.169 + static const Word mul = RandomTraits<Word>::mul;
6.170 +
6.171 + current = state;
6.172 +
6.173 + Word *curr = state + length - 1;
6.174 + curr[0] = seed; --curr;
6.175 + for (int i = 1; i < length; ++i) {
6.176 + curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i);
6.177 + --curr;
6.178 + }
6.179 + }
6.180 +
6.181 + template <typename Iterator>
6.182 + void initState(Iterator begin, Iterator end) {
6.183 +
6.184 + static const Word init = RandomTraits<Word>::arrayInit;
6.185 + static const Word mul1 = RandomTraits<Word>::arrayMul1;
6.186 + static const Word mul2 = RandomTraits<Word>::arrayMul2;
6.187 +
6.188 +
6.189 + Word *curr = state + length - 1; --curr;
6.190 + Iterator it = begin; int cnt = 0;
6.191 + int num;
6.192 +
6.193 + initState(init);
6.194 +
6.195 + num = length > end - begin ? length : end - begin;
6.196 + while (num--) {
6.197 + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
6.198 + + *it + cnt;
6.199 + ++it; ++cnt;
6.200 + if (it == end) {
6.201 + it = begin; cnt = 0;
6.202 + }
6.203 + if (curr == state) {
6.204 + curr = state + length - 1; curr[0] = state[0];
6.205 + }
6.206 + --curr;
6.207 + }
6.208 +
6.209 + num = length - 1; cnt = length - (curr - state) - 1;
6.210 + while (num--) {
6.211 + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2))
6.212 + - cnt;
6.213 + --curr; ++cnt;
6.214 + if (curr == state) {
6.215 + curr = state + length - 1; curr[0] = state[0]; --curr;
6.216 + cnt = 1;
6.217 + }
6.218 + }
6.219 +
6.220 + state[length - 1] = Word(1) << (bits - 1);
6.221 + }
6.222 +
6.223 + void copyState(const RandomCore& other) {
6.224 + std::copy(other.state, other.state + length, state);
6.225 + current = state + (other.current - other.state);
6.226 + }
6.227 +
6.228 + Word operator()() {
6.229 + if (current == state) fillState();
6.230 + --current;
6.231 + Word rnd = *current;
6.232 + return RandomTraits<Word>::tempering(rnd);
6.233 + }
6.234 +
6.235 + private:
6.236 +
6.237 +
6.238 + void fillState() {
6.239 + static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
6.240 + static const Word loMask = RandomTraits<Word>::loMask;
6.241 + static const Word hiMask = RandomTraits<Word>::hiMask;
6.242 +
6.243 + current = state + length;
6.244 +
6.245 + register Word *curr = state + length - 1;
6.246 + register long num;
6.247 +
6.248 + num = length - shift;
6.249 + while (num--) {
6.250 + curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
6.251 + curr[- shift] ^ mask[curr[-1] & 1ul];
6.252 + --curr;
6.253 + }
6.254 + num = shift - 1;
6.255 + while (num--) {
6.256 + curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^
6.257 + curr[length - shift] ^ mask[curr[-1] & 1ul];
6.258 + --curr;
6.259 + }
6.260 + curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^
6.261 + curr[length - shift] ^ mask[curr[length - 1] & 1ul];
6.262 +
6.263 + }
6.264 +
6.265 +
6.266 + Word *current;
6.267 + Word state[length];
6.268 +
6.269 + };
6.270 +
6.271 +
6.272 + template <typename Result,
6.273 + int shift = (std::numeric_limits<Result>::digits + 1) / 2>
6.274 + struct Masker {
6.275 + static Result mask(const Result& result) {
6.276 + return Masker<Result, (shift + 1) / 2>::
6.277 + mask(static_cast<Result>(result | (result >> shift)));
6.278 + }
6.279 + };
6.280 +
6.281 + template <typename Result>
6.282 + struct Masker<Result, 1> {
6.283 + static Result mask(const Result& result) {
6.284 + return static_cast<Result>(result | (result >> 1));
6.285 + }
6.286 + };
6.287 +
6.288 + template <typename Result, typename Word,
6.289 + int rest = std::numeric_limits<Result>::digits, int shift = 0,
6.290 + bool last = rest <= std::numeric_limits<Word>::digits>
6.291 + struct IntConversion {
6.292 + static const int bits = std::numeric_limits<Word>::digits;
6.293 +
6.294 + static Result convert(RandomCore<Word>& rnd) {
6.295 + return static_cast<Result>(rnd() >> (bits - rest)) << shift;
6.296 + }
6.297 +
6.298 + };
6.299 +
6.300 + template <typename Result, typename Word, int rest, int shift>
6.301 + struct IntConversion<Result, Word, rest, shift, false> {
6.302 + static const int bits = std::numeric_limits<Word>::digits;
6.303 +
6.304 + static Result convert(RandomCore<Word>& rnd) {
6.305 + return (static_cast<Result>(rnd()) << shift) |
6.306 + IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
6.307 + }
6.308 + };
6.309 +
6.310 +
6.311 + template <typename Result, typename Word,
6.312 + bool one_word = (std::numeric_limits<Word>::digits <
6.313 + std::numeric_limits<Result>::digits) >
6.314 + struct Mapping {
6.315 + static Result map(RandomCore<Word>& rnd, const Result& bound) {
6.316 + Word max = Word(bound - 1);
6.317 + Result mask = Masker<Result>::mask(bound - 1);
6.318 + Result num;
6.319 + do {
6.320 + num = IntConversion<Result, Word>::convert(rnd) & mask;
6.321 + } while (num > max);
6.322 + return num;
6.323 + }
6.324 + };
6.325 +
6.326 + template <typename Result, typename Word>
6.327 + struct Mapping<Result, Word, false> {
6.328 + static Result map(RandomCore<Word>& rnd, const Result& bound) {
6.329 + Word max = Word(bound - 1);
6.330 + Word mask = Masker<Word, (std::numeric_limits<Result>::digits + 1) / 2>
6.331 + ::mask(max);
6.332 + Word num;
6.333 + do {
6.334 + num = rnd() & mask;
6.335 + } while (num > max);
6.336 + return num;
6.337 + }
6.338 + };
6.339 +
6.340 + template <typename Result, int exp, bool pos = (exp >= 0)>
6.341 + struct ShiftMultiplier {
6.342 + static const Result multiplier() {
6.343 + Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
6.344 + res *= res;
6.345 + if ((exp & 1) == 1) res *= static_cast<Result>(2.0);
6.346 + return res;
6.347 + }
6.348 + };
6.349 +
6.350 + template <typename Result, int exp>
6.351 + struct ShiftMultiplier<Result, exp, false> {
6.352 + static const Result multiplier() {
6.353 + Result res = ShiftMultiplier<Result, exp / 2>::multiplier();
6.354 + res *= res;
6.355 + if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
6.356 + return res;
6.357 + }
6.358 + };
6.359 +
6.360 + template <typename Result>
6.361 + struct ShiftMultiplier<Result, 0, true> {
6.362 + static const Result multiplier() {
6.363 + return static_cast<Result>(1.0);
6.364 + }
6.365 + };
6.366 +
6.367 + template <typename Result>
6.368 + struct ShiftMultiplier<Result, -20, true> {
6.369 + static const Result multiplier() {
6.370 + return static_cast<Result>(1.0/1048576.0);
6.371 + }
6.372 + };
6.373 +
6.374 + template <typename Result>
6.375 + struct ShiftMultiplier<Result, -32, true> {
6.376 + static const Result multiplier() {
6.377 + return static_cast<Result>(1.0/424967296.0);
6.378 + }
6.379 + };
6.380 +
6.381 + template <typename Result>
6.382 + struct ShiftMultiplier<Result, -53, true> {
6.383 + static const Result multiplier() {
6.384 + return static_cast<Result>(1.0/9007199254740992.0);
6.385 + }
6.386 + };
6.387 +
6.388 + template <typename Result>
6.389 + struct ShiftMultiplier<Result, -64, true> {
6.390 + static const Result multiplier() {
6.391 + return static_cast<Result>(1.0/18446744073709551616.0);
6.392 + }
6.393 + };
6.394 +
6.395 + template <typename Result, int exp>
6.396 + struct Shifting {
6.397 + static Result shift(const Result& result) {
6.398 + return result * ShiftMultiplier<Result, exp>::multiplier();
6.399 + }
6.400 + };
6.401 +
6.402 + template <typename Result, typename Word,
6.403 + int rest = std::numeric_limits<Result>::digits, int shift = 0,
6.404 + bool last = rest <= std::numeric_limits<Word>::digits>
6.405 + struct RealConversion{
6.406 + static const int bits = std::numeric_limits<Word>::digits;
6.407 +
6.408 + static Result convert(RandomCore<Word>& rnd) {
6.409 + return Shifting<Result, - shift - rest>::
6.410 + shift(static_cast<Result>(rnd() >> (bits - rest)));
6.411 + }
6.412 + };
6.413 +
6.414 + template <typename Result, typename Word, int rest, int shift>
6.415 + struct RealConversion<Result, Word, rest, shift, false> {
6.416 + static const int bits = std::numeric_limits<Word>::digits;
6.417 +
6.418 + static Result convert(RandomCore<Word>& rnd) {
6.419 + return Shifting<Result, - shift - bits>::
6.420 + shift(static_cast<Result>(rnd())) +
6.421 + RealConversion<Result, Word, rest-bits, shift + bits>::
6.422 + convert(rnd);
6.423 + }
6.424 + };
6.425 +
6.426 + template <typename Result, typename Word>
6.427 + struct Initializer {
6.428 +
6.429 + template <typename Iterator>
6.430 + static void init(RandomCore<Word>& rnd, Iterator begin, Iterator end) {
6.431 + std::vector<Word> ws;
6.432 + for (Iterator it = begin; it != end; ++it) {
6.433 + ws.push_back(Word(*it));
6.434 + }
6.435 + rnd.initState(ws.begin(), ws.end());
6.436 + }
6.437 +
6.438 + static void init(RandomCore<Word>& rnd, Result seed) {
6.439 + rnd.initState(seed);
6.440 + }
6.441 + };
6.442 +
6.443 + template <typename Word>
6.444 + struct BoolConversion {
6.445 + static bool convert(RandomCore<Word>& rnd) {
6.446 + return (rnd() & 1) == 1;
6.447 + }
6.448 + };
6.449 +
6.450 + template <typename Word>
6.451 + struct BoolProducer {
6.452 + Word buffer;
6.453 + int num;
6.454 +
6.455 + BoolProducer() : num(0) {}
6.456 +
6.457 + bool convert(RandomCore<Word>& rnd) {
6.458 + if (num == 0) {
6.459 + buffer = rnd();
6.460 + num = RandomTraits<Word>::bits;
6.461 + }
6.462 + bool r = (buffer & 1);
6.463 + buffer >>= 1;
6.464 + --num;
6.465 + return r;
6.466 + }
6.467 + };
6.468 +
6.469 + }
6.470 +
6.471 + /// \ingroup misc
6.472 + ///
6.473 + /// \brief Mersenne Twister random number generator
6.474 + ///
6.475 + /// The Mersenne Twister is a twisted generalized feedback
6.476 + /// shift-register generator of Matsumoto and Nishimura. The period
6.477 + /// of this generator is \f$ 2^{19937} - 1 \f$ and it is
6.478 + /// equi-distributed in 623 dimensions for 32-bit numbers. The time
6.479 + /// performance of this generator is comparable to the commonly used
6.480 + /// generators.
6.481 + ///
6.482 + /// This implementation is specialized for both 32-bit and 64-bit
6.483 + /// architectures. The generators differ sligthly in the
6.484 + /// initialization and generation phase so they produce two
6.485 + /// completly different sequences.
6.486 + ///
6.487 + /// The generator gives back random numbers of serveral types. To
6.488 + /// get a random number from a range of a floating point type you
6.489 + /// can use one form of the \c operator() or the \c real() member
6.490 + /// function. If you want to get random number from the {0, 1, ...,
6.491 + /// n-1} integer range use the \c operator[] or the \c integer()
6.492 + /// method. And to get random number from the whole range of an
6.493 + /// integer type you can use the argumentless \c integer() or \c
6.494 + /// uinteger() functions. After all you can get random bool with
6.495 + /// equal chance of true and false or given probability of true
6.496 + /// result with the \c boolean() member functions.
6.497 + ///
6.498 + ///\code
6.499 + /// // The commented code is identical to the other
6.500 + /// double a = rnd(); // [0.0, 1.0)
6.501 + /// // double a = rnd.real(); // [0.0, 1.0)
6.502 + /// double b = rnd(100.0); // [0.0, 100.0)
6.503 + /// // double b = rnd.real(100.0); // [0.0, 100.0)
6.504 + /// double c = rnd(1.0, 2.0); // [1.0, 2.0)
6.505 + /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0)
6.506 + /// int d = rnd[100000]; // 0..99999
6.507 + /// // int d = rnd.integer(100000); // 0..99999
6.508 + /// int e = rnd[6] + 1; // 1..6
6.509 + /// // int e = rnd.integer(1, 1 + 6); // 1..6
6.510 + /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1
6.511 + /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1
6.512 + /// bool g = rnd.boolean(); // P(g = true) = 0.5
6.513 + /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8
6.514 + ///\endcode
6.515 + ///
6.516 + /// The lemon provides a global instance of the random number
6.517 + /// generator which name is \ref lemon::rnd "rnd". Usually it is a
6.518 + /// good programming convenience to use this global generator to get
6.519 + /// random numbers.
6.520 + class Random {
6.521 + private:
6.522 +
6.523 + // Architecture word
6.524 + typedef unsigned long Word;
6.525 +
6.526 + _random_bits::RandomCore<Word> core;
6.527 + _random_bits::BoolProducer<Word> bool_producer;
6.528 +
6.529 +
6.530 + public:
6.531 +
6.532 + /// \brief Constructor
6.533 + ///
6.534 + /// Constructor with constant seeding.
6.535 + Random() { core.initState(); }
6.536 +
6.537 + /// \brief Constructor
6.538 + ///
6.539 + /// Constructor with seed. The current number type will be converted
6.540 + /// to the architecture word type.
6.541 + template <typename Number>
6.542 + Random(Number seed) {
6.543 + _random_bits::Initializer<Number, Word>::init(core, seed);
6.544 + }
6.545 +
6.546 + /// \brief Constructor
6.547 + ///
6.548 + /// Constructor with array seeding. The given range should contain
6.549 + /// any number type and the numbers will be converted to the
6.550 + /// architecture word type.
6.551 + template <typename Iterator>
6.552 + Random(Iterator begin, Iterator end) {
6.553 + typedef typename std::iterator_traits<Iterator>::value_type Number;
6.554 + _random_bits::Initializer<Number, Word>::init(core, begin, end);
6.555 + }
6.556 +
6.557 + /// \brief Copy constructor
6.558 + ///
6.559 + /// Copy constructor. The generated sequence will be identical to
6.560 + /// the other sequence. It can be used to save the current state
6.561 + /// of the generator and later use it to generate the same
6.562 + /// sequence.
6.563 + Random(const Random& other) {
6.564 + core.copyState(other.core);
6.565 + }
6.566 +
6.567 + /// \brief Assign operator
6.568 + ///
6.569 + /// Assign operator. The generated sequence will be identical to
6.570 + /// the other sequence. It can be used to save the current state
6.571 + /// of the generator and later use it to generate the same
6.572 + /// sequence.
6.573 + Random& operator=(const Random& other) {
6.574 + if (&other != this) {
6.575 + core.copyState(other.core);
6.576 + }
6.577 + return *this;
6.578 + }
6.579 +
6.580 + /// \brief Returns a random real number from the range [0, 1)
6.581 + ///
6.582 + /// It returns a random real number from the range [0, 1). The
6.583 + /// default Number type is double.
6.584 + template <typename Number>
6.585 + Number real() {
6.586 + return _random_bits::RealConversion<Number, Word>::convert(core);
6.587 + }
6.588 +
6.589 + double real() {
6.590 + return real<double>();
6.591 + }
6.592 +
6.593 + /// \brief Returns a random real number the range [0, b)
6.594 + ///
6.595 + /// It returns a random real number from the range [0, b).
6.596 + template <typename Number>
6.597 + Number real(Number b) {
6.598 + return real<Number>() * b;
6.599 + }
6.600 +
6.601 + /// \brief Returns a random real number from the range [a, b)
6.602 + ///
6.603 + /// It returns a random real number from the range [a, b).
6.604 + template <typename Number>
6.605 + Number real(Number a, Number b) {
6.606 + return real<Number>() * (b - a) + a;
6.607 + }
6.608 +
6.609 + /// \brief Returns a random real number from the range [0, 1)
6.610 + ///
6.611 + /// It returns a random double from the range [0, 1).
6.612 + double operator()() {
6.613 + return real<double>();
6.614 + }
6.615 +
6.616 + /// \brief Returns a random real number from the range [0, b)
6.617 + ///
6.618 + /// It returns a random real number from the range [0, b).
6.619 + template <typename Number>
6.620 + Number operator()(Number b) {
6.621 + return real<Number>() * b;
6.622 + }
6.623 +
6.624 + /// \brief Returns a random real number from the range [a, b)
6.625 + ///
6.626 + /// It returns a random real number from the range [a, b).
6.627 + template <typename Number>
6.628 + Number operator()(Number a, Number b) {
6.629 + return real<Number>() * (b - a) + a;
6.630 + }
6.631 +
6.632 + /// \brief Returns a random integer from a range
6.633 + ///
6.634 + /// It returns a random integer from the range {0, 1, ..., b - 1}.
6.635 + template <typename Number>
6.636 + Number integer(Number b) {
6.637 + return _random_bits::Mapping<Number, Word>::map(core, b);
6.638 + }
6.639 +
6.640 + /// \brief Returns a random integer from a range
6.641 + ///
6.642 + /// It returns a random integer from the range {a, a + 1, ..., b - 1}.
6.643 + template <typename Number>
6.644 + Number integer(Number a, Number b) {
6.645 + return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;
6.646 + }
6.647 +
6.648 + /// \brief Returns a random integer from a range
6.649 + ///
6.650 + /// It returns a random integer from the range {0, 1, ..., b - 1}.
6.651 + template <typename Number>
6.652 + Number operator[](Number b) {
6.653 + return _random_bits::Mapping<Number, Word>::map(core, b);
6.654 + }
6.655 +
6.656 + /// \brief Returns a random non-negative integer
6.657 + ///
6.658 + /// It returns a random non-negative integer uniformly from the
6.659 + /// whole range of the current \c Number type. The default result
6.660 + /// type of this function is unsigned int.
6.661 + template <typename Number>
6.662 + Number uinteger() {
6.663 + return _random_bits::IntConversion<Number, Word>::convert(core);
6.664 + }
6.665 +
6.666 + unsigned int uinteger() {
6.667 + return uinteger<unsigned int>();
6.668 + }
6.669 +
6.670 + /// \brief Returns a random integer
6.671 + ///
6.672 + /// It returns a random integer uniformly from the whole range of
6.673 + /// the current \c Number type. The default result type of this
6.674 + /// function is int.
6.675 + template <typename Number>
6.676 + Number integer() {
6.677 + static const int nb = std::numeric_limits<Number>::digits +
6.678 + (std::numeric_limits<Number>::is_signed ? 1 : 0);
6.679 + return _random_bits::IntConversion<Number, Word, nb>::convert(core);
6.680 + }
6.681 +
6.682 + int integer() {
6.683 + return integer<int>();
6.684 + }
6.685 +
6.686 + /// \brief Returns a random bool
6.687 + ///
6.688 + /// It returns a random bool. The generator holds a buffer for
6.689 + /// random bits. Every time when it become empty the generator makes
6.690 + /// a new random word and fill the buffer up.
6.691 + bool boolean() {
6.692 + return bool_producer.convert(core);
6.693 + }
6.694 +
6.695 + ///\name Nonuniform distributions
6.696 + ///
6.697 +
6.698 + ///@{
6.699 +
6.700 + /// \brief Returns a random bool
6.701 + ///
6.702 + /// It returns a random bool with given probability of true result
6.703 + bool boolean(double p) {
6.704 + return operator()() < p;
6.705 + }
6.706 +
6.707 + /// Standard Gauss distribution
6.708 +
6.709 + /// Standard Gauss distribution.
6.710 + /// \note The Cartesian form of the Box-Muller
6.711 + /// transformation is used to generate a random normal distribution.
6.712 + /// \todo Consider using the "ziggurat" method instead.
6.713 + double gauss()
6.714 + {
6.715 + double V1,V2,S;
6.716 + do {
6.717 + V1=2*real<double>()-1;
6.718 + V2=2*real<double>()-1;
6.719 + S=V1*V1+V2*V2;
6.720 + } while(S>=1);
6.721 + return std::sqrt(-2*std::log(S)/S)*V1;
6.722 + }
6.723 + /// Gauss distribution with given mean and standard deviation
6.724 +
6.725 + /// Gauss distribution with given mean and standard deviation
6.726 + /// \sa gauss()
6.727 + double gauss(double mean,double std_dev)
6.728 + {
6.729 + return gauss()*std_dev+mean;
6.730 + }
6.731 +
6.732 + /// Exponential distribution with given mean
6.733 +
6.734 + /// This function generates an exponential distribution random number
6.735 + /// with mean <tt>1/lambda</tt>.
6.736 + ///
6.737 + double exponential(double lambda=1.0)
6.738 + {
6.739 + return -std::log(1.0-real<double>())/lambda;
6.740 + }
6.741 +
6.742 + /// Gamma distribution with given integer shape
6.743 +
6.744 + /// This function generates a gamma distribution random number.
6.745 + ///
6.746 + ///\param k shape parameter (<tt>k>0</tt> integer)
6.747 + double gamma(int k)
6.748 + {
6.749 + double s = 0;
6.750 + for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
6.751 + return s;
6.752 + }
6.753 +
6.754 + /// Gamma distribution with given shape and scale parameter
6.755 +
6.756 + /// This function generates a gamma distribution random number.
6.757 + ///
6.758 + ///\param k shape parameter (<tt>k>0</tt>)
6.759 + ///\param theta scale parameter
6.760 + ///
6.761 + double gamma(double k,double theta=1.0)
6.762 + {
6.763 + double xi,nu;
6.764 + const double delta = k-std::floor(k);
6.765 + const double v0=M_E/(M_E-delta);
6.766 + do {
6.767 + double V0=1.0-real<double>();
6.768 + double V1=1.0-real<double>();
6.769 + double V2=1.0-real<double>();
6.770 + if(V2<=v0)
6.771 + {
6.772 + xi=std::pow(V1,1.0/delta);
6.773 + nu=V0*std::pow(xi,delta-1.0);
6.774 + }
6.775 + else
6.776 + {
6.777 + xi=1.0-std::log(V1);
6.778 + nu=V0*std::exp(-xi);
6.779 + }
6.780 + } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
6.781 + return theta*(xi-gamma(int(std::floor(k))));
6.782 + }
6.783 +
6.784 + /// Weibull distribution
6.785 +
6.786 + /// This function generates a Weibull distribution random number.
6.787 + ///
6.788 + ///\param k shape parameter (<tt>k>0</tt>)
6.789 + ///\param lambda scale parameter (<tt>lambda>0</tt>)
6.790 + ///
6.791 + double weibull(double k,double lambda)
6.792 + {
6.793 + return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
6.794 + }
6.795 +
6.796 + /// Pareto distribution
6.797 +
6.798 + /// This function generates a Pareto distribution random number.
6.799 + ///
6.800 + ///\param k shape parameter (<tt>k>0</tt>)
6.801 + ///\param x_min location parameter (<tt>x_min>0</tt>)
6.802 + ///
6.803 + double pareto(double k,double x_min)
6.804 + {
6.805 + return exponential(gamma(k,1.0/x_min));
6.806 + }
6.807 +
6.808 + ///@}
6.809 +
6.810 + ///\name Two dimensional distributions
6.811 + ///
6.812 +
6.813 + ///@{
6.814 +
6.815 + /// Uniform distribution on the full unit circle.
6.816 +
6.817 + /// Uniform distribution on the full unit circle.
6.818 + ///
6.819 + dim2::Point<double> disc()
6.820 + {
6.821 + double V1,V2;
6.822 + do {
6.823 + V1=2*real<double>()-1;
6.824 + V2=2*real<double>()-1;
6.825 +
6.826 + } while(V1*V1+V2*V2>=1);
6.827 + return dim2::Point<double>(V1,V2);
6.828 + }
6.829 + /// A kind of two dimensional Gauss distribution
6.830 +
6.831 + /// This function provides a turning symmetric two-dimensional distribution.
6.832 + /// Both coordinates are of standard normal distribution, but they are not
6.833 + /// independent.
6.834 + ///
6.835 + /// \note The coordinates are the two random variables provided by
6.836 + /// the Box-Muller method.
6.837 + dim2::Point<double> gauss2()
6.838 + {
6.839 + double V1,V2,S;
6.840 + do {
6.841 + V1=2*real<double>()-1;
6.842 + V2=2*real<double>()-1;
6.843 + S=V1*V1+V2*V2;
6.844 + } while(S>=1);
6.845 + double W=std::sqrt(-2*std::log(S)/S);
6.846 + return dim2::Point<double>(W*V1,W*V2);
6.847 + }
6.848 + /// A kind of two dimensional exponential distribution
6.849 +
6.850 + /// This function provides a turning symmetric two-dimensional distribution.
6.851 + /// The x-coordinate is of conditionally exponential distribution
6.852 + /// with the condition that x is positive and y=0. If x is negative and
6.853 + /// y=0 then, -x is of exponential distribution. The same is true for the
6.854 + /// y-coordinate.
6.855 + dim2::Point<double> exponential2()
6.856 + {
6.857 + double V1,V2,S;
6.858 + do {
6.859 + V1=2*real<double>()-1;
6.860 + V2=2*real<double>()-1;
6.861 + S=V1*V1+V2*V2;
6.862 + } while(S>=1);
6.863 + double W=-std::log(S)/S;
6.864 + return dim2::Point<double>(W*V1,W*V2);
6.865 + }
6.866 +
6.867 + ///@}
6.868 + };
6.869 +
6.870 +
6.871 + extern Random rnd;
6.872 +
6.873 +}
6.874 +
6.875 +#endif
7.1 --- a/lemon/tolerance.h Thu Jan 03 11:12:27 2008 +0100
7.2 +++ b/lemon/tolerance.h Thu Jan 03 11:13:29 2008 +0100
7.3 @@ -48,9 +48,13 @@
7.4 ///\sa Tolerance<double>
7.5 ///\sa Tolerance<long double>
7.6 ///\sa Tolerance<int>
7.7 +#if defined __GNUC__ && !defined __STRICT_ANSI__
7.8 ///\sa Tolerance<long long int>
7.9 +#endif
7.10 ///\sa Tolerance<unsigned int>
7.11 +#if defined __GNUC__ && !defined __STRICT_ANSI__
7.12 ///\sa Tolerance<unsigned long long int>
7.13 +#endif
7.14
7.15 template<class T>
7.16 class Tolerance
7.17 @@ -130,7 +134,7 @@
7.18 ///Returns \c true if \c a is \e surely negative
7.19 bool negative(Value a) const { return -_epsilon>a; }
7.20 ///Returns \c true if \c a is \e surely non-zero
7.21 - bool nonZero(Value a) const { return positive(a)||negative(a); };
7.22 + bool nonZero(Value a) const { return positive(a)||negative(a); }
7.23
7.24 ///@}
7.25
7.26 @@ -181,7 +185,7 @@
7.27 ///Returns \c true if \c a is \e surely negative
7.28 bool negative(Value a) const { return -_epsilon>a; }
7.29 ///Returns \c true if \c a is \e surely non-zero
7.30 - bool nonZero(Value a) const { return positive(a)||negative(a); };
7.31 + bool nonZero(Value a) const { return positive(a)||negative(a); }
7.32
7.33 ///@}
7.34
7.35 @@ -232,7 +236,7 @@
7.36 ///Returns \c true if \c a is \e surely negative
7.37 bool negative(Value a) const { return -_epsilon>a; }
7.38 ///Returns \c true if \c a is \e surely non-zero
7.39 - bool nonZero(Value a) const { return positive(a)||negative(a); };
7.40 + bool nonZero(Value a) const { return positive(a)||negative(a); }
7.41
7.42 ///@}
7.43
7.44 @@ -265,7 +269,7 @@
7.45 ///Returns \c true if \c a is \e surely negative
7.46 static bool negative(Value a) { return 0>a; }
7.47 ///Returns \c true if \c a is \e surely non-zero
7.48 - static bool nonZero(Value a) { return a!=0; };
7.49 + static bool nonZero(Value a) { return a!=0; }
7.50
7.51 ///@}
7.52
7.53 @@ -298,7 +302,7 @@
7.54 ///Returns \c true if \c a is \e surely negative
7.55 static bool negative(Value) { return false; }
7.56 ///Returns \c true if \c a is \e surely non-zero
7.57 - static bool nonZero(Value a) { return a!=0; };
7.58 + static bool nonZero(Value a) { return a!=0; }
7.59
7.60 ///@}
7.61
7.62 @@ -332,7 +336,7 @@
7.63 ///Returns \c true if \c a is \e surely negative
7.64 static bool negative(Value a) { return 0>a; }
7.65 ///Returns \c true if \c a is \e surely non-zero
7.66 - static bool nonZero(Value a) { return a!=0;};
7.67 + static bool nonZero(Value a) { return a!=0;}
7.68
7.69 ///@}
7.70
7.71 @@ -365,7 +369,7 @@
7.72 ///Returns \c true if \c a is \e surely negative
7.73 static bool negative(Value) { return false; }
7.74 ///Returns \c true if \c a is \e surely non-zero
7.75 - static bool nonZero(Value a) { return a!=0;};
7.76 + static bool nonZero(Value a) { return a!=0;}
7.77
7.78 ///@}
7.79
7.80 @@ -402,7 +406,7 @@
7.81 ///Returns \c true if \c a is \e surely negative
7.82 static bool negative(Value a) { return 0>a; }
7.83 ///Returns \c true if \c a is \e surely non-zero
7.84 - static bool nonZero(Value a) { return a!=0;};
7.85 + static bool nonZero(Value a) { return a!=0;}
7.86
7.87 ///@}
7.88
7.89 @@ -437,7 +441,7 @@
7.90 ///Returns \c true if \c a is \e surely negative
7.91 static bool negative(Value) { return false; }
7.92 ///Returns \c true if \c a is \e surely non-zero
7.93 - static bool nonZero(Value a) { return a!=0;};
7.94 + static bool nonZero(Value a) { return a!=0;}
7.95
7.96 ///@}
7.97
8.1 --- a/test/Makefile.am Thu Jan 03 11:12:27 2008 +0100
8.2 +++ b/test/Makefile.am Thu Jan 03 11:13:29 2008 +0100
8.3 @@ -3,15 +3,17 @@
8.4
8.5 noinst_HEADERS += \
8.6 test/test_tools.h
8.7 -
8.8 +
8.9 check_PROGRAMS += \
8.10 test/dim_test \
8.11 + test/random_test \
8.12 test/test_tools_fail \
8.13 test/test_tools_pass
8.14 -
8.15 +
8.16 TESTS += $(check_PROGRAMS)
8.17 XFAIL_TESTS += test/test_tools_fail$(EXEEXT)
8.18
8.19 test_dim_test_SOURCES = test/dim_test.cc
8.20 +test_random_test_SOURCES = test/random_test.cc
8.21 test_test_tools_fail_SOURCES = test/test_tools_fail.cc
8.22 test_test_tools_pass_SOURCES = test/test_tools_pass.cc
9.1 --- a/test/dim_test.cc Thu Jan 03 11:12:27 2008 +0100
9.2 +++ b/test/dim_test.cc Thu Jan 03 11:13:29 2008 +0100
9.3 @@ -22,65 +22,66 @@
9.4
9.5 using namespace std;
9.6 using namespace lemon;
9.7 +
9.8 int main()
9.9 {
9.10 -
9.11 - cout << "Testing classes `dim2::Point' and `dim2::BoundingBox'." << endl;
9.12 + cout << "Testing classes 'dim2::Point' and 'dim2::BoundingBox'." << endl;
9.13
9.14 typedef dim2::Point<int> Point;
9.15 -
9.16 - Point seged;
9.17 - check(seged.size()==2, "Wrong vector addition");
9.18 +
9.19 + Point p;
9.20 + check(p.size()==2, "Wrong vector initialization.");
9.21
9.22 Point a(1,2);
9.23 Point b(3,4);
9.24 + check(a[0]==1 && a[1]==2, "Wrong vector initialization.");
9.25
9.26 - check(a[0]==1 && a[1]==2, "Wrong vector addition");
9.27 + p = a+b;
9.28 + check(p.x==4 && p.y==6, "Wrong vector addition.");
9.29
9.30 - seged = a+b;
9.31 - check(seged.x==4 && seged.y==6, "Wrong vector addition");
9.32 + p = a-b;
9.33 + check(p.x==-2 && p.y==-2, "Wrong vector subtraction.");
9.34
9.35 - seged = a-b;
9.36 - check(seged.x==-2 && seged.y==-2, "a-b");
9.37 -
9.38 - check(a.normSquare()==5,"Wrong norm calculation");
9.39 - check(a*b==11, "a*b");
9.40 + check(a.normSquare()==5,"Wrong vector norm calculation.");
9.41 + check(a*b==11, "Wrong vector scalar product.");
9.42
9.43 int l=2;
9.44 - seged = a*l;
9.45 - check(seged.x==2 && seged.y==4, "a*l");
9.46 + p = a*l;
9.47 + check(p.x==2 && p.y==4, "Wrong vector multiplication by a scalar.");
9.48
9.49 - seged = b/l;
9.50 - check(seged.x==1 && seged.y==2, "b/l");
9.51 + p = b/l;
9.52 + check(p.x==1 && p.y==2, "Wrong vector division by a scalar.");
9.53
9.54 typedef dim2::BoundingBox<int> BB;
9.55 - BB doboz1;
9.56 - check(doboz1.empty(), "It should be empty.");
9.57 -
9.58 - doboz1.add(a);
9.59 - check(!doboz1.empty(), "It should not be empty.");
9.60 - doboz1.add(b);
9.61 + BB box1;
9.62 + check(box1.empty(), "It should be empty.");
9.63
9.64 - check(doboz1.bottomLeft().x==1 &&
9.65 - doboz1.bottomLeft().y==2 &&
9.66 - doboz1.topRight().x==3 &&
9.67 - doboz1.topRight().y==4,
9.68 - "added points to box");
9.69 + box1.add(a);
9.70 + check(!box1.empty(), "It should not be empty.");
9.71 + box1.add(b);
9.72
9.73 - seged.x=2;seged.y=3;
9.74 - check(doboz1.inside(seged),"It should be inside.");
9.75 + check(box1.bottomLeft().x==1 &&
9.76 + box1.bottomLeft().y==2 &&
9.77 + box1.topRight().x==3 &&
9.78 + box1.topRight().y==4,
9.79 + "Wrong addition of points to box.");
9.80
9.81 - seged.x=1;seged.y=3;
9.82 - check(doboz1.inside(seged),"It should be inside.");
9.83 + p.x=2; p.y=3;
9.84 + check(box1.inside(p), "It should be inside.");
9.85
9.86 - seged.x=0;seged.y=3;
9.87 - check(!doboz1.inside(seged),"It should not be inside.");
9.88 + p.x=1; p.y=3;
9.89 + check(box1.inside(p), "It should be inside.");
9.90
9.91 - BB doboz2(seged);
9.92 - check(!doboz2.empty(),
9.93 + p.x=0; p.y=3;
9.94 + check(!box1.inside(p), "It should not be inside.");
9.95 +
9.96 + BB box2(p);
9.97 + check(!box2.empty(),
9.98 "It should not be empty. Constructed from 1 point.");
9.99
9.100 - doboz2.add(doboz1);
9.101 - check(doboz2.inside(seged),
9.102 + box2.add(box1);
9.103 + check(box2.inside(p),
9.104 "It should be inside. Incremented a box with another one.");
9.105 +
9.106 + return 0;
9.107 }
10.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
10.2 +++ b/test/random_test.cc Thu Jan 03 11:13:29 2008 +0100
10.3 @@ -0,0 +1,36 @@
10.4 +/* -*- C++ -*-
10.5 + *
10.6 + * This file is a part of LEMON, a generic C++ optimization library
10.7 + *
10.8 + * Copyright (C) 2003-2007
10.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
10.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
10.11 + *
10.12 + * Permission to use, modify and distribute this software is granted
10.13 + * provided that this copyright notice appears in all copies. For
10.14 + * precise terms see the accompanying LICENSE file.
10.15 + *
10.16 + * This software is provided "AS IS" with no warranty of any kind,
10.17 + * express or implied, and with no claim as to its suitability for any
10.18 + * purpose.
10.19 + *
10.20 + */
10.21 +
10.22 +#include <lemon/random.h>
10.23 +#include "test_tools.h"
10.24 +
10.25 +///\file \brief Test cases for random.h
10.26 +///
10.27 +///\todo To be extended
10.28 +///
10.29 +
10.30 +int main()
10.31 +{
10.32 + double a=lemon::rnd();
10.33 + check(a<1.0&&a>0.0,"This should be in [0,1)");
10.34 + a=lemon::rnd.gauss();
10.35 + a=lemon::rnd.gamma(3.45,0);
10.36 + a=lemon::rnd.gamma(4);
10.37 + //Does gamma work with integer k?
10.38 + a=lemon::rnd.gamma(4.0,0);
10.39 +}