# HG changeset patch # User Alpar Juttner # Date 1199355209 -3600 # Node ID 45f8b617339ede421b3776fde736505a48be1488 # Parent 40d6f625e549fabdc92788f7cdb20aceea5f8b78# Parent 9244aa9ec12e88ef84bb6a217f83637092d899e7 Merge diff -r 40d6f625e549 -r 45f8b617339e .hgignore --- a/.hgignore Thu Jan 03 11:12:27 2008 +0100 +++ b/.hgignore Thu Jan 03 11:13:29 2008 +0100 @@ -5,6 +5,9 @@ *~ *.o .#.* +*.log +*.lo +*.tar.* Makefile.in aclocal.m4 config.h.in @@ -19,11 +22,13 @@ lemon/libemon.la lemon/stamp-h2 doc/Doxyfile -lemon/.dirstamp -lemon/.libs/* +.dirstamp +.libs/* +.deps/* syntax: regexp -html/.* -autom4te.cache/.* -build-aux/.* -objs.*/.* +^doc/html/.* +^autom4te.cache/.* +^build-aux/.* +^objs.*/.* +^test/[a-z_]*$ \ No newline at end of file diff -r 40d6f625e549 -r 45f8b617339e lemon/Makefile.am --- a/lemon/Makefile.am Thu Jan 03 11:12:27 2008 +0100 +++ b/lemon/Makefile.am Thu Jan 03 11:13:29 2008 +0100 @@ -7,13 +7,16 @@ lib_LTLIBRARIES += lemon/libemon.la lemon_libemon_la_SOURCES = \ - lemon/base.cc + lemon/base.cc \ + lemon/random.cc + lemon_libemon_la_CXXFLAGS = $(GLPK_CFLAGS) $(CPLEX_CFLAGS) $(SOPLEX_CXXFLAGS) lemon_libemon_la_LDFLAGS = $(GLPK_LIBS) $(CPLEX_LIBS) $(SOPLEX_LIBS) lemon_HEADERS += \ lemon/dim2.h \ + lemon/random.h \ lemon/list_graph.h \ lemon/tolerance.h diff -r 40d6f625e549 -r 45f8b617339e lemon/bits/invalid.h --- a/lemon/bits/invalid.h Thu Jan 03 11:12:27 2008 +0100 +++ b/lemon/bits/invalid.h Thu Jan 03 11:13:29 2008 +0100 @@ -24,7 +24,7 @@ namespace lemon { - /// \brief Dummy type to make it easier to make invalid iterators. + /// \brief Dummy type to make it easier to create invalid iterators. /// /// See \ref INVALID for the usage. struct Invalid { @@ -34,8 +34,8 @@ bool operator< (Invalid) { return false; } }; - /// Invalid iterators. - + /// \brief Invalid iterators. + /// /// \ref Invalid is a global type that converts to each iterator /// in such a way that the value of the target iterator will be invalid. diff -r 40d6f625e549 -r 45f8b617339e lemon/dim2.h --- a/lemon/dim2.h Thu Jan 03 11:12:27 2008 +0100 +++ b/lemon/dim2.h Thu Jan 03 11:13:29 2008 +0100 @@ -34,9 +34,6 @@ /// can be used to determine /// the rectangular bounding box of a set of /// \ref lemon::dim2::Point "dim2::Point"'s. -/// -///\author Attila Bernath - namespace lemon { @@ -62,9 +59,9 @@ typedef T Value; - ///First co-ordinate + ///First coordinate T x; - ///Second co-ordinate + ///Second coordinate T y; ///Default constructor @@ -75,8 +72,8 @@ ///The dimension of the vector. - ///This class give back always 2. - /// + ///The dimension of the vector. + ///This function always returns 2. int size() const { return 2; } ///Subscripting operator @@ -138,7 +135,7 @@ return b+=u; } - ///Return the neg of the vectors + ///Return the negative of the vector Point operator-() const { Point b=*this; b.x=-b.x; b.y=-b.y; @@ -175,9 +172,9 @@ }; - ///Return an Point + ///Return a Point - ///Return an Point + ///Return a Point. ///\relates Point template inline Point makePoint(const T& x, const T& y) { @@ -186,7 +183,7 @@ ///Return a vector multiplied by a scalar - ///Return a vector multiplied by a scalar + ///Return a vector multiplied by a scalar. ///\relates Point template Point operator*(const T &u,const Point &x) { return x*u; @@ -194,7 +191,7 @@ ///Read a plainvector from a stream - ///Read a plainvector from a stream + ///Read a plainvector from a stream. ///\relates Point /// template @@ -222,7 +219,7 @@ ///Write a plainvector to a stream - ///Write a plainvector to a stream + ///Write a plainvector to a stream. ///\relates Point /// template @@ -234,7 +231,7 @@ ///Rotate by 90 degrees - ///Returns its parameter rotated by 90 degrees in positive direction. + ///Returns the parameter rotated by 90 degrees in positive direction. ///\relates Point /// template @@ -245,7 +242,7 @@ ///Rotate by 180 degrees - ///Returns its parameter rotated by 180 degrees. + ///Returns the parameter rotated by 180 degrees. ///\relates Point /// template @@ -256,7 +253,7 @@ ///Rotate by 270 degrees - ///Returns its parameter rotated by 90 degrees in negative direction. + ///Returns the parameter rotated by 90 degrees in negative direction. ///\relates Point /// template @@ -271,7 +268,6 @@ /// A class to calculate or store the bounding box of plainvectors. /// - ///\author Attila Bernath template class BoundingBox { Point bottom_left, top_right; @@ -286,9 +282,11 @@ ///Construct an instance from two points - ///Construct an instance from two points - ///\warning The coordinates of the bottom-left corner must be no more - ///than those of the top-right one + ///Construct an instance from two points. + ///\param a The bottom left corner. + ///\param b The top right corner. + ///\warning The coordinates of the bottom left corner must be no more + ///than those of the top right one. BoundingBox(Point a,Point b) { bottom_left=a; @@ -298,9 +296,13 @@ ///Construct an instance from four numbers - ///Construct an instance from four numbers - ///\warning The coordinates of the bottom-left corner must be no more - ///than those of the top-right one + ///Construct an instance from four numbers. + ///\param l The left side of the box. + ///\param b The bottom of the box. + ///\param r The right side of the box. + ///\param t The top of the box. + ///\warning The left side must be no more than the right side and + ///bottom must be no more than the top. BoundingBox(T l,T b,T r,T t) { bottom_left=Point(l,b); @@ -308,7 +310,12 @@ _empty = false; } - ///Were any points added? + ///Return \c true if the bounding box is empty. + + ///Return \c true if the bounding box is empty (i.e. return \c false + ///if at least one point was added to the box or the coordinates of + ///the box were set). + ///The coordinates of an empty bounding box are not defined. bool empty() const { return _empty; } @@ -329,7 +336,7 @@ ///Set the bottom left corner ///Set the bottom left corner. - ///It should only bee used for non-empty box. + ///It should only be used for non-empty box. void bottomLeft(Point p) { bottom_left = p; } @@ -345,7 +352,7 @@ ///Set the top right corner ///Set the top right corner. - ///It should only bee used for non-empty box. + ///It should only be used for non-empty box. void topRight(Point p) { top_right = p; } @@ -361,7 +368,7 @@ ///Set the bottom right corner ///Set the bottom right corner. - ///It should only bee used for non-empty box. + ///It should only be used for non-empty box. void bottomRight(Point p) { top_right.x = p.x; bottom_left.y = p.y; @@ -378,7 +385,7 @@ ///Set the top left corner ///Set the top left corner. - ///It should only bee used for non-empty box. + ///It should only be used for non-empty box. void topLeft(Point p) { top_right.y = p.y; bottom_left.x = p.x; @@ -395,7 +402,7 @@ ///Set the bottom of the box ///Set the bottom of the box. - ///It should only bee used for non-empty box. + ///It should only be used for non-empty box. void bottom(T t) { bottom_left.y = t; } @@ -411,7 +418,7 @@ ///Set the top of the box ///Set the top of the box. - ///It should only bee used for non-empty box. + ///It should only be used for non-empty box. void top(T t) { top_right.y = t; } @@ -427,7 +434,7 @@ ///Set the left side of the box ///Set the left side of the box. - ///It should only bee used for non-empty box + ///It should only be used for non-empty box. void left(T t) { bottom_left.x = t; } @@ -443,7 +450,7 @@ ///Set the right side of the box ///Set the right side of the box. - ///It should only bee used for non-empty box + ///It should only be used for non-empty box. void right(T t) { top_right.x = t; } @@ -465,7 +472,7 @@ } ///Checks whether a point is inside a bounding box - bool inside(const Point& u){ + bool inside(const Point& u) const { if (_empty) return false; else{ @@ -475,6 +482,9 @@ } ///Increments a bounding box with a point + + ///Increments a bounding box with a point. + /// BoundingBox& add(const Point& u){ if (_empty){ bottom_left=top_right=u; @@ -489,7 +499,10 @@ return *this; } - ///Increments a bounding to contain another bounding box + ///Increments a bounding box to contain another bounding box + + ///Increments a bounding box to contain another bounding box. + /// BoundingBox& add(const BoundingBox &u){ if ( !u.empty() ){ this->add(u.bottomLeft()); @@ -499,24 +512,31 @@ } ///Intersection of two bounding boxes - BoundingBox operator &(const BoundingBox& u){ + + ///Intersection of two bounding boxes. + /// + BoundingBox operator&(const BoundingBox& u) const { BoundingBox b; - b.bottom_left.x=std::max(this->bottom_left.x,u.bottom_left.x); - b.bottom_left.y=std::max(this->bottom_left.y,u.bottom_left.y); - b.top_right.x=std::min(this->top_right.x,u.top_right.x); - b.top_right.y=std::min(this->top_right.y,u.top_right.y); - b._empty = this->_empty || u._empty || - b.bottom_left.x>top_right.x && b.bottom_left.y>top_right.y; + if (this->_empty || u._empty) { + b._empty = true; + } else { + b.bottom_left.x = std::max(this->bottom_left.x,u.bottom_left.x); + b.bottom_left.y = std::max(this->bottom_left.y,u.bottom_left.y); + b.top_right.x = std::min(this->top_right.x,u.top_right.x); + b.top_right.y = std::min(this->top_right.y,u.top_right.y); + b._empty = b.bottom_left.x > b.top_right.x || + b.bottom_left.y > b.top_right.y; + } return b; } };//class Boundingbox - ///Map of x-coordinates of a dim2::Point<>-map + ///Map of x-coordinates of a \ref Point "Point"-map ///\ingroup maps - ///Map of x-coordinates of a dim2::Point<>-map + ///Map of x-coordinates of a \ref Point "Point"-map. /// template class XMap @@ -570,7 +590,7 @@ ///Returns a \ref ConstXMap class - ///This function just returns an \ref ConstXMap class. + ///This function just returns a \ref ConstXMap class. /// ///\ingroup maps ///\relates ConstXMap @@ -580,10 +600,10 @@ return ConstXMap(m); } - ///Map of y-coordinates of a dim2::Point<>-map + ///Map of y-coordinates of a \ref Point "Point"-map ///\ingroup maps - ///Map of y-coordinates of a dim2::Point<>-map + ///Map of y-coordinates of a \ref Point "Point"-map. /// template class YMap @@ -599,9 +619,9 @@ void set(Key k,Value v) {_map.set(k,typename M::Value(_map[k].x,v));} }; - ///Returns an \ref YMap class + ///Returns a \ref YMap class - ///This function just returns an \ref YMap class. + ///This function just returns a \ref YMap class. /// ///\ingroup maps ///\relates YMap @@ -637,7 +657,7 @@ ///Returns a \ref ConstYMap class - ///This function just returns an \ref ConstYMap class. + ///This function just returns a \ref ConstYMap class. /// ///\ingroup maps ///\relates ConstYMap @@ -649,10 +669,10 @@ ///\brief Map of the \ref Point::normSquare() "normSquare()" - ///of an \ref Point "Point"-map + ///of a \ref Point "Point"-map /// ///Map of the \ref Point::normSquare() "normSquare()" - ///of an \ref Point "Point"-map + ///of a \ref Point "Point"-map. ///\ingroup maps /// template @@ -670,7 +690,7 @@ ///Returns a \ref NormSquareMap class - ///This function just returns an \ref NormSquareMap class. + ///This function just returns a \ref NormSquareMap class. /// ///\ingroup maps ///\relates NormSquareMap diff -r 40d6f625e549 -r 45f8b617339e lemon/random.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lemon/random.cc Thu Jan 03 11:13:29 2008 +0100 @@ -0,0 +1,29 @@ +/* -*- C++ -*- + * + * This file is a part of LEMON, a generic C++ optimization library + * + * Copyright (C) 2003-2007 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +///\file +///\brief Instantiation of the Random class. + +#include + +namespace lemon { + /// \brief Global random number generator instance + /// + /// A global Mersenne Twister random number generator instance. + Random rnd; +} diff -r 40d6f625e549 -r 45f8b617339e lemon/random.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lemon/random.h Thu Jan 03 11:13:29 2008 +0100 @@ -0,0 +1,872 @@ +/* -*- C++ -*- + * + * This file is a part of LEMON, a generic C++ optimization library + * + * Copyright (C) 2003-2007 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +/* + * This file contains the reimplemented version of the Mersenne Twister + * Generator of Matsumoto and Nishimura. + * + * See the appropriate copyright notice below. + * + * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * 3. The names of its contributors may not be used to endorse or promote + * products derived from this software without specific prior written + * permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS + * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, + * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED + * OF THE POSSIBILITY OF SUCH DAMAGE. + * + * + * Any feedback is very welcome. + * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) + */ + +#ifndef LEMON_RANDOM_H +#define LEMON_RANDOM_H + +#include +#include +#include + +#include +#include + +#include +///\ingroup misc +///\file +///\brief Mersenne Twister random number generator + +namespace lemon { + + namespace _random_bits { + + template ::digits> + struct RandomTraits {}; + + template + struct RandomTraits<_Word, 32> { + + typedef _Word Word; + static const int bits = 32; + + static const int length = 624; + static const int shift = 397; + + static const Word mul = 0x6c078965u; + static const Word arrayInit = 0x012BD6AAu; + static const Word arrayMul1 = 0x0019660Du; + static const Word arrayMul2 = 0x5D588B65u; + + static const Word mask = 0x9908B0DFu; + static const Word loMask = (1u << 31) - 1; + static const Word hiMask = ~loMask; + + + static Word tempering(Word rnd) { + rnd ^= (rnd >> 11); + rnd ^= (rnd << 7) & 0x9D2C5680u; + rnd ^= (rnd << 15) & 0xEFC60000u; + rnd ^= (rnd >> 18); + return rnd; + } + + }; + + template + struct RandomTraits<_Word, 64> { + + typedef _Word Word; + static const int bits = 64; + + static const int length = 312; + static const int shift = 156; + + static const Word mul = Word(0x5851F42Du) << 32 | Word(0x4C957F2Du); + static const Word arrayInit = Word(0x00000000u) << 32 |Word(0x012BD6AAu); + static const Word arrayMul1 = Word(0x369DEA0Fu) << 32 |Word(0x31A53F85u); + static const Word arrayMul2 = Word(0x27BB2EE6u) << 32 |Word(0x87B0B0FDu); + + static const Word mask = Word(0xB5026F5Au) << 32 | Word(0xA96619E9u); + static const Word loMask = (Word(1u) << 31) - 1; + static const Word hiMask = ~loMask; + + static Word tempering(Word rnd) { + rnd ^= (rnd >> 29) & (Word(0x55555555u) << 32 | Word(0x55555555u)); + rnd ^= (rnd << 17) & (Word(0x71D67FFFu) << 32 | Word(0xEDA60000u)); + rnd ^= (rnd << 37) & (Word(0xFFF7EEE0u) << 32 | Word(0x00000000u)); + rnd ^= (rnd >> 43); + return rnd; + } + + }; + + template + class RandomCore { + public: + + typedef _Word Word; + + private: + + static const int bits = RandomTraits::bits; + + static const int length = RandomTraits::length; + static const int shift = RandomTraits::shift; + + public: + + void initState() { + static const Word seedArray[4] = { + 0x12345u, 0x23456u, 0x34567u, 0x45678u + }; + + initState(seedArray, seedArray + 4); + } + + void initState(Word seed) { + + static const Word mul = RandomTraits::mul; + + current = state; + + Word *curr = state + length - 1; + curr[0] = seed; --curr; + for (int i = 1; i < length; ++i) { + curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i); + --curr; + } + } + + template + void initState(Iterator begin, Iterator end) { + + static const Word init = RandomTraits::arrayInit; + static const Word mul1 = RandomTraits::arrayMul1; + static const Word mul2 = RandomTraits::arrayMul2; + + + Word *curr = state + length - 1; --curr; + Iterator it = begin; int cnt = 0; + int num; + + initState(init); + + num = length > end - begin ? length : end - begin; + while (num--) { + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) + + *it + cnt; + ++it; ++cnt; + if (it == end) { + it = begin; cnt = 0; + } + if (curr == state) { + curr = state + length - 1; curr[0] = state[0]; + } + --curr; + } + + num = length - 1; cnt = length - (curr - state) - 1; + while (num--) { + curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul2)) + - cnt; + --curr; ++cnt; + if (curr == state) { + curr = state + length - 1; curr[0] = state[0]; --curr; + cnt = 1; + } + } + + state[length - 1] = Word(1) << (bits - 1); + } + + void copyState(const RandomCore& other) { + std::copy(other.state, other.state + length, state); + current = state + (other.current - other.state); + } + + Word operator()() { + if (current == state) fillState(); + --current; + Word rnd = *current; + return RandomTraits::tempering(rnd); + } + + private: + + + void fillState() { + static const Word mask[2] = { 0x0ul, RandomTraits::mask }; + static const Word loMask = RandomTraits::loMask; + static const Word hiMask = RandomTraits::hiMask; + + current = state + length; + + register Word *curr = state + length - 1; + register long num; + + num = length - shift; + while (num--) { + curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ + curr[- shift] ^ mask[curr[-1] & 1ul]; + --curr; + } + num = shift - 1; + while (num--) { + curr[0] = (((curr[0] & hiMask) | (curr[-1] & loMask)) >> 1) ^ + curr[length - shift] ^ mask[curr[-1] & 1ul]; + --curr; + } + curr[0] = (((curr[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ + curr[length - shift] ^ mask[curr[length - 1] & 1ul]; + + } + + + Word *current; + Word state[length]; + + }; + + + template ::digits + 1) / 2> + struct Masker { + static Result mask(const Result& result) { + return Masker:: + mask(static_cast(result | (result >> shift))); + } + }; + + template + struct Masker { + static Result mask(const Result& result) { + return static_cast(result | (result >> 1)); + } + }; + + template ::digits, int shift = 0, + bool last = rest <= std::numeric_limits::digits> + struct IntConversion { + static const int bits = std::numeric_limits::digits; + + static Result convert(RandomCore& rnd) { + return static_cast(rnd() >> (bits - rest)) << shift; + } + + }; + + template + struct IntConversion { + static const int bits = std::numeric_limits::digits; + + static Result convert(RandomCore& rnd) { + return (static_cast(rnd()) << shift) | + IntConversion::convert(rnd); + } + }; + + + template ::digits < + std::numeric_limits::digits) > + struct Mapping { + static Result map(RandomCore& rnd, const Result& bound) { + Word max = Word(bound - 1); + Result mask = Masker::mask(bound - 1); + Result num; + do { + num = IntConversion::convert(rnd) & mask; + } while (num > max); + return num; + } + }; + + template + struct Mapping { + static Result map(RandomCore& rnd, const Result& bound) { + Word max = Word(bound - 1); + Word mask = Masker::digits + 1) / 2> + ::mask(max); + Word num; + do { + num = rnd() & mask; + } while (num > max); + return num; + } + }; + + template = 0)> + struct ShiftMultiplier { + static const Result multiplier() { + Result res = ShiftMultiplier::multiplier(); + res *= res; + if ((exp & 1) == 1) res *= static_cast(2.0); + return res; + } + }; + + template + struct ShiftMultiplier { + static const Result multiplier() { + Result res = ShiftMultiplier::multiplier(); + res *= res; + if ((exp & 1) == 1) res *= static_cast(0.5); + return res; + } + }; + + template + struct ShiftMultiplier { + static const Result multiplier() { + return static_cast(1.0); + } + }; + + template + struct ShiftMultiplier { + static const Result multiplier() { + return static_cast(1.0/1048576.0); + } + }; + + template + struct ShiftMultiplier { + static const Result multiplier() { + return static_cast(1.0/424967296.0); + } + }; + + template + struct ShiftMultiplier { + static const Result multiplier() { + return static_cast(1.0/9007199254740992.0); + } + }; + + template + struct ShiftMultiplier { + static const Result multiplier() { + return static_cast(1.0/18446744073709551616.0); + } + }; + + template + struct Shifting { + static Result shift(const Result& result) { + return result * ShiftMultiplier::multiplier(); + } + }; + + template ::digits, int shift = 0, + bool last = rest <= std::numeric_limits::digits> + struct RealConversion{ + static const int bits = std::numeric_limits::digits; + + static Result convert(RandomCore& rnd) { + return Shifting:: + shift(static_cast(rnd() >> (bits - rest))); + } + }; + + template + struct RealConversion { + static const int bits = std::numeric_limits::digits; + + static Result convert(RandomCore& rnd) { + return Shifting:: + shift(static_cast(rnd())) + + RealConversion:: + convert(rnd); + } + }; + + template + struct Initializer { + + template + static void init(RandomCore& rnd, Iterator begin, Iterator end) { + std::vector ws; + for (Iterator it = begin; it != end; ++it) { + ws.push_back(Word(*it)); + } + rnd.initState(ws.begin(), ws.end()); + } + + static void init(RandomCore& rnd, Result seed) { + rnd.initState(seed); + } + }; + + template + struct BoolConversion { + static bool convert(RandomCore& rnd) { + return (rnd() & 1) == 1; + } + }; + + template + struct BoolProducer { + Word buffer; + int num; + + BoolProducer() : num(0) {} + + bool convert(RandomCore& rnd) { + if (num == 0) { + buffer = rnd(); + num = RandomTraits::bits; + } + bool r = (buffer & 1); + buffer >>= 1; + --num; + return r; + } + }; + + } + + /// \ingroup misc + /// + /// \brief Mersenne Twister random number generator + /// + /// The Mersenne Twister is a twisted generalized feedback + /// shift-register generator of Matsumoto and Nishimura. The period + /// of this generator is \f$ 2^{19937} - 1 \f$ and it is + /// equi-distributed in 623 dimensions for 32-bit numbers. The time + /// performance of this generator is comparable to the commonly used + /// generators. + /// + /// This implementation is specialized for both 32-bit and 64-bit + /// architectures. The generators differ sligthly in the + /// initialization and generation phase so they produce two + /// completly different sequences. + /// + /// The generator gives back random numbers of serveral types. To + /// get a random number from a range of a floating point type you + /// can use one form of the \c operator() or the \c real() member + /// function. If you want to get random number from the {0, 1, ..., + /// n-1} integer range use the \c operator[] or the \c integer() + /// method. And to get random number from the whole range of an + /// integer type you can use the argumentless \c integer() or \c + /// uinteger() functions. After all you can get random bool with + /// equal chance of true and false or given probability of true + /// result with the \c boolean() member functions. + /// + ///\code + /// // The commented code is identical to the other + /// double a = rnd(); // [0.0, 1.0) + /// // double a = rnd.real(); // [0.0, 1.0) + /// double b = rnd(100.0); // [0.0, 100.0) + /// // double b = rnd.real(100.0); // [0.0, 100.0) + /// double c = rnd(1.0, 2.0); // [1.0, 2.0) + /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) + /// int d = rnd[100000]; // 0..99999 + /// // int d = rnd.integer(100000); // 0..99999 + /// int e = rnd[6] + 1; // 1..6 + /// // int e = rnd.integer(1, 1 + 6); // 1..6 + /// int b = rnd.uinteger(); // 0 .. 2^31 - 1 + /// int c = rnd.integer(); // - 2^31 .. 2^31 - 1 + /// bool g = rnd.boolean(); // P(g = true) = 0.5 + /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 + ///\endcode + /// + /// The lemon provides a global instance of the random number + /// generator which name is \ref lemon::rnd "rnd". Usually it is a + /// good programming convenience to use this global generator to get + /// random numbers. + class Random { + private: + + // Architecture word + typedef unsigned long Word; + + _random_bits::RandomCore core; + _random_bits::BoolProducer bool_producer; + + + public: + + /// \brief Constructor + /// + /// Constructor with constant seeding. + Random() { core.initState(); } + + /// \brief Constructor + /// + /// Constructor with seed. The current number type will be converted + /// to the architecture word type. + template + Random(Number seed) { + _random_bits::Initializer::init(core, seed); + } + + /// \brief Constructor + /// + /// Constructor with array seeding. The given range should contain + /// any number type and the numbers will be converted to the + /// architecture word type. + template + Random(Iterator begin, Iterator end) { + typedef typename std::iterator_traits::value_type Number; + _random_bits::Initializer::init(core, begin, end); + } + + /// \brief Copy constructor + /// + /// Copy constructor. The generated sequence will be identical to + /// the other sequence. It can be used to save the current state + /// of the generator and later use it to generate the same + /// sequence. + Random(const Random& other) { + core.copyState(other.core); + } + + /// \brief Assign operator + /// + /// Assign operator. The generated sequence will be identical to + /// the other sequence. It can be used to save the current state + /// of the generator and later use it to generate the same + /// sequence. + Random& operator=(const Random& other) { + if (&other != this) { + core.copyState(other.core); + } + return *this; + } + + /// \brief Returns a random real number from the range [0, 1) + /// + /// It returns a random real number from the range [0, 1). The + /// default Number type is double. + template + Number real() { + return _random_bits::RealConversion::convert(core); + } + + double real() { + return real(); + } + + /// \brief Returns a random real number the range [0, b) + /// + /// It returns a random real number from the range [0, b). + template + Number real(Number b) { + return real() * b; + } + + /// \brief Returns a random real number from the range [a, b) + /// + /// It returns a random real number from the range [a, b). + template + Number real(Number a, Number b) { + return real() * (b - a) + a; + } + + /// \brief Returns a random real number from the range [0, 1) + /// + /// It returns a random double from the range [0, 1). + double operator()() { + return real(); + } + + /// \brief Returns a random real number from the range [0, b) + /// + /// It returns a random real number from the range [0, b). + template + Number operator()(Number b) { + return real() * b; + } + + /// \brief Returns a random real number from the range [a, b) + /// + /// It returns a random real number from the range [a, b). + template + Number operator()(Number a, Number b) { + return real() * (b - a) + a; + } + + /// \brief Returns a random integer from a range + /// + /// It returns a random integer from the range {0, 1, ..., b - 1}. + template + Number integer(Number b) { + return _random_bits::Mapping::map(core, b); + } + + /// \brief Returns a random integer from a range + /// + /// It returns a random integer from the range {a, a + 1, ..., b - 1}. + template + Number integer(Number a, Number b) { + return _random_bits::Mapping::map(core, b - a) + a; + } + + /// \brief Returns a random integer from a range + /// + /// It returns a random integer from the range {0, 1, ..., b - 1}. + template + Number operator[](Number b) { + return _random_bits::Mapping::map(core, b); + } + + /// \brief Returns a random non-negative integer + /// + /// It returns a random non-negative integer uniformly from the + /// whole range of the current \c Number type. The default result + /// type of this function is unsigned int. + template + Number uinteger() { + return _random_bits::IntConversion::convert(core); + } + + unsigned int uinteger() { + return uinteger(); + } + + /// \brief Returns a random integer + /// + /// It returns a random integer uniformly from the whole range of + /// the current \c Number type. The default result type of this + /// function is int. + template + Number integer() { + static const int nb = std::numeric_limits::digits + + (std::numeric_limits::is_signed ? 1 : 0); + return _random_bits::IntConversion::convert(core); + } + + int integer() { + return integer(); + } + + /// \brief Returns a random bool + /// + /// It returns a random bool. The generator holds a buffer for + /// random bits. Every time when it become empty the generator makes + /// a new random word and fill the buffer up. + bool boolean() { + return bool_producer.convert(core); + } + + ///\name Nonuniform distributions + /// + + ///@{ + + /// \brief Returns a random bool + /// + /// It returns a random bool with given probability of true result + bool boolean(double p) { + return operator()() < p; + } + + /// Standard Gauss distribution + + /// Standard Gauss distribution. + /// \note The Cartesian form of the Box-Muller + /// transformation is used to generate a random normal distribution. + /// \todo Consider using the "ziggurat" method instead. + double gauss() + { + double V1,V2,S; + do { + V1=2*real()-1; + V2=2*real()-1; + S=V1*V1+V2*V2; + } while(S>=1); + return std::sqrt(-2*std::log(S)/S)*V1; + } + /// Gauss distribution with given mean and standard deviation + + /// Gauss distribution with given mean and standard deviation + /// \sa gauss() + double gauss(double mean,double std_dev) + { + return gauss()*std_dev+mean; + } + + /// Exponential distribution with given mean + + /// This function generates an exponential distribution random number + /// with mean 1/lambda. + /// + double exponential(double lambda=1.0) + { + return -std::log(1.0-real())/lambda; + } + + /// Gamma distribution with given integer shape + + /// This function generates a gamma distribution random number. + /// + ///\param k shape parameter (k>0 integer) + double gamma(int k) + { + double s = 0; + for(int i=0;i()); + return s; + } + + /// Gamma distribution with given shape and scale parameter + + /// This function generates a gamma distribution random number. + /// + ///\param k shape parameter (k>0) + ///\param theta scale parameter + /// + double gamma(double k,double theta=1.0) + { + double xi,nu; + const double delta = k-std::floor(k); + const double v0=M_E/(M_E-delta); + do { + double V0=1.0-real(); + double V1=1.0-real(); + double V2=1.0-real(); + if(V2<=v0) + { + xi=std::pow(V1,1.0/delta); + nu=V0*std::pow(xi,delta-1.0); + } + else + { + xi=1.0-std::log(V1); + nu=V0*std::exp(-xi); + } + } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); + return theta*(xi-gamma(int(std::floor(k)))); + } + + /// Weibull distribution + + /// This function generates a Weibull distribution random number. + /// + ///\param k shape parameter (k>0) + ///\param lambda scale parameter (lambda>0) + /// + double weibull(double k,double lambda) + { + return lambda*pow(-std::log(1.0-real()),1.0/k); + } + + /// Pareto distribution + + /// This function generates a Pareto distribution random number. + /// + ///\param k shape parameter (k>0) + ///\param x_min location parameter (x_min>0) + /// + double pareto(double k,double x_min) + { + return exponential(gamma(k,1.0/x_min)); + } + + ///@} + + ///\name Two dimensional distributions + /// + + ///@{ + + /// Uniform distribution on the full unit circle. + + /// Uniform distribution on the full unit circle. + /// + dim2::Point disc() + { + double V1,V2; + do { + V1=2*real()-1; + V2=2*real()-1; + + } while(V1*V1+V2*V2>=1); + return dim2::Point(V1,V2); + } + /// A kind of two dimensional Gauss distribution + + /// This function provides a turning symmetric two-dimensional distribution. + /// Both coordinates are of standard normal distribution, but they are not + /// independent. + /// + /// \note The coordinates are the two random variables provided by + /// the Box-Muller method. + dim2::Point gauss2() + { + double V1,V2,S; + do { + V1=2*real()-1; + V2=2*real()-1; + S=V1*V1+V2*V2; + } while(S>=1); + double W=std::sqrt(-2*std::log(S)/S); + return dim2::Point(W*V1,W*V2); + } + /// A kind of two dimensional exponential distribution + + /// This function provides a turning symmetric two-dimensional distribution. + /// The x-coordinate is of conditionally exponential distribution + /// with the condition that x is positive and y=0. If x is negative and + /// y=0 then, -x is of exponential distribution. The same is true for the + /// y-coordinate. + dim2::Point exponential2() + { + double V1,V2,S; + do { + V1=2*real()-1; + V2=2*real()-1; + S=V1*V1+V2*V2; + } while(S>=1); + double W=-std::log(S)/S; + return dim2::Point(W*V1,W*V2); + } + + ///@} + }; + + + extern Random rnd; + +} + +#endif diff -r 40d6f625e549 -r 45f8b617339e lemon/tolerance.h --- a/lemon/tolerance.h Thu Jan 03 11:12:27 2008 +0100 +++ b/lemon/tolerance.h Thu Jan 03 11:13:29 2008 +0100 @@ -48,9 +48,13 @@ ///\sa Tolerance ///\sa Tolerance ///\sa Tolerance +#if defined __GNUC__ && !defined __STRICT_ANSI__ ///\sa Tolerance +#endif ///\sa Tolerance +#if defined __GNUC__ && !defined __STRICT_ANSI__ ///\sa Tolerance +#endif template class Tolerance @@ -130,7 +134,7 @@ ///Returns \c true if \c a is \e surely negative bool negative(Value a) const { return -_epsilon>a; } ///Returns \c true if \c a is \e surely non-zero - bool nonZero(Value a) const { return positive(a)||negative(a); }; + bool nonZero(Value a) const { return positive(a)||negative(a); } ///@} @@ -181,7 +185,7 @@ ///Returns \c true if \c a is \e surely negative bool negative(Value a) const { return -_epsilon>a; } ///Returns \c true if \c a is \e surely non-zero - bool nonZero(Value a) const { return positive(a)||negative(a); }; + bool nonZero(Value a) const { return positive(a)||negative(a); } ///@} @@ -232,7 +236,7 @@ ///Returns \c true if \c a is \e surely negative bool negative(Value a) const { return -_epsilon>a; } ///Returns \c true if \c a is \e surely non-zero - bool nonZero(Value a) const { return positive(a)||negative(a); }; + bool nonZero(Value a) const { return positive(a)||negative(a); } ///@} @@ -265,7 +269,7 @@ ///Returns \c true if \c a is \e surely negative static bool negative(Value a) { return 0>a; } ///Returns \c true if \c a is \e surely non-zero - static bool nonZero(Value a) { return a!=0; }; + static bool nonZero(Value a) { return a!=0; } ///@} @@ -298,7 +302,7 @@ ///Returns \c true if \c a is \e surely negative static bool negative(Value) { return false; } ///Returns \c true if \c a is \e surely non-zero - static bool nonZero(Value a) { return a!=0; }; + static bool nonZero(Value a) { return a!=0; } ///@} @@ -332,7 +336,7 @@ ///Returns \c true if \c a is \e surely negative static bool negative(Value a) { return 0>a; } ///Returns \c true if \c a is \e surely non-zero - static bool nonZero(Value a) { return a!=0;}; + static bool nonZero(Value a) { return a!=0;} ///@} @@ -365,7 +369,7 @@ ///Returns \c true if \c a is \e surely negative static bool negative(Value) { return false; } ///Returns \c true if \c a is \e surely non-zero - static bool nonZero(Value a) { return a!=0;}; + static bool nonZero(Value a) { return a!=0;} ///@} @@ -402,7 +406,7 @@ ///Returns \c true if \c a is \e surely negative static bool negative(Value a) { return 0>a; } ///Returns \c true if \c a is \e surely non-zero - static bool nonZero(Value a) { return a!=0;}; + static bool nonZero(Value a) { return a!=0;} ///@} @@ -437,7 +441,7 @@ ///Returns \c true if \c a is \e surely negative static bool negative(Value) { return false; } ///Returns \c true if \c a is \e surely non-zero - static bool nonZero(Value a) { return a!=0;}; + static bool nonZero(Value a) { return a!=0;} ///@} diff -r 40d6f625e549 -r 45f8b617339e test/Makefile.am --- a/test/Makefile.am Thu Jan 03 11:12:27 2008 +0100 +++ b/test/Makefile.am Thu Jan 03 11:13:29 2008 +0100 @@ -3,15 +3,17 @@ noinst_HEADERS += \ test/test_tools.h - + check_PROGRAMS += \ test/dim_test \ + test/random_test \ test/test_tools_fail \ test/test_tools_pass - + TESTS += $(check_PROGRAMS) XFAIL_TESTS += test/test_tools_fail$(EXEEXT) test_dim_test_SOURCES = test/dim_test.cc +test_random_test_SOURCES = test/random_test.cc test_test_tools_fail_SOURCES = test/test_tools_fail.cc test_test_tools_pass_SOURCES = test/test_tools_pass.cc diff -r 40d6f625e549 -r 45f8b617339e test/dim_test.cc --- a/test/dim_test.cc Thu Jan 03 11:12:27 2008 +0100 +++ b/test/dim_test.cc Thu Jan 03 11:13:29 2008 +0100 @@ -22,65 +22,66 @@ using namespace std; using namespace lemon; + int main() { - - cout << "Testing classes `dim2::Point' and `dim2::BoundingBox'." << endl; + cout << "Testing classes 'dim2::Point' and 'dim2::BoundingBox'." << endl; typedef dim2::Point Point; - - Point seged; - check(seged.size()==2, "Wrong vector addition"); + + Point p; + check(p.size()==2, "Wrong vector initialization."); Point a(1,2); Point b(3,4); + check(a[0]==1 && a[1]==2, "Wrong vector initialization."); - check(a[0]==1 && a[1]==2, "Wrong vector addition"); + p = a+b; + check(p.x==4 && p.y==6, "Wrong vector addition."); - seged = a+b; - check(seged.x==4 && seged.y==6, "Wrong vector addition"); + p = a-b; + check(p.x==-2 && p.y==-2, "Wrong vector subtraction."); - seged = a-b; - check(seged.x==-2 && seged.y==-2, "a-b"); - - check(a.normSquare()==5,"Wrong norm calculation"); - check(a*b==11, "a*b"); + check(a.normSquare()==5,"Wrong vector norm calculation."); + check(a*b==11, "Wrong vector scalar product."); int l=2; - seged = a*l; - check(seged.x==2 && seged.y==4, "a*l"); + p = a*l; + check(p.x==2 && p.y==4, "Wrong vector multiplication by a scalar."); - seged = b/l; - check(seged.x==1 && seged.y==2, "b/l"); + p = b/l; + check(p.x==1 && p.y==2, "Wrong vector division by a scalar."); typedef dim2::BoundingBox BB; - BB doboz1; - check(doboz1.empty(), "It should be empty."); - - doboz1.add(a); - check(!doboz1.empty(), "It should not be empty."); - doboz1.add(b); + BB box1; + check(box1.empty(), "It should be empty."); - check(doboz1.bottomLeft().x==1 && - doboz1.bottomLeft().y==2 && - doboz1.topRight().x==3 && - doboz1.topRight().y==4, - "added points to box"); + box1.add(a); + check(!box1.empty(), "It should not be empty."); + box1.add(b); - seged.x=2;seged.y=3; - check(doboz1.inside(seged),"It should be inside."); + check(box1.bottomLeft().x==1 && + box1.bottomLeft().y==2 && + box1.topRight().x==3 && + box1.topRight().y==4, + "Wrong addition of points to box."); - seged.x=1;seged.y=3; - check(doboz1.inside(seged),"It should be inside."); + p.x=2; p.y=3; + check(box1.inside(p), "It should be inside."); - seged.x=0;seged.y=3; - check(!doboz1.inside(seged),"It should not be inside."); + p.x=1; p.y=3; + check(box1.inside(p), "It should be inside."); - BB doboz2(seged); - check(!doboz2.empty(), + p.x=0; p.y=3; + check(!box1.inside(p), "It should not be inside."); + + BB box2(p); + check(!box2.empty(), "It should not be empty. Constructed from 1 point."); - doboz2.add(doboz1); - check(doboz2.inside(seged), + box2.add(box1); + check(box2.inside(p), "It should be inside. Incremented a box with another one."); + + return 0; } diff -r 40d6f625e549 -r 45f8b617339e test/random_test.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/random_test.cc Thu Jan 03 11:13:29 2008 +0100 @@ -0,0 +1,36 @@ +/* -*- C++ -*- + * + * This file is a part of LEMON, a generic C++ optimization library + * + * Copyright (C) 2003-2007 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include "test_tools.h" + +///\file \brief Test cases for random.h +/// +///\todo To be extended +/// + +int main() +{ + double a=lemon::rnd(); + check(a<1.0&&a>0.0,"This should be in [0,1)"); + a=lemon::rnd.gauss(); + a=lemon::rnd.gamma(3.45,0); + a=lemon::rnd.gamma(4); + //Does gamma work with integer k? + a=lemon::rnd.gamma(4.0,0); +}