Changeset 209:765619b7cbb2 in lemon for lemon/random.h
 Timestamp:
 07/13/08 20:51:02 (16 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/random.h
r178 r209 1 /* * C++*2 * 3 * This file is a part of LEMON, a generic C++ optimization library 1 /* * mode: C++; indenttabsmode: nil; * 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 5 * Copyright (C) 20032008 … … 22 22 * 23 23 * See the appropriate copyright notice below. 24 * 24 * 25 25 * Copyright (C) 1997  2002, Makoto Matsumoto and Takuji Nishimura, 26 * All rights reserved. 26 * All rights reserved. 27 27 * 28 28 * Redistribution and use in source and binary forms, with or without … … 37 37 * documentation and/or other materials provided with the distribution. 38 38 * 39 * 3. The names of its contributors may not be used to endorse or promote 40 * products derived from this software without specific prior written 39 * 3. The names of its contributors may not be used to endorse or promote 40 * products derived from this software without specific prior written 41 41 * permission. 42 42 * … … 88 88 89 89 namespace _random_bits { 90 90 91 91 template <typename _Word, int _bits = std::numeric_limits<_Word>::digits> 92 92 struct RandomTraits {}; … … 100 100 static const int length = 624; 101 101 static const int shift = 397; 102 102 103 103 static const Word mul = 0x6c078965u; 104 104 static const Word arrayInit = 0x012BD6AAu; … … 168 168 0x12345u, 0x23456u, 0x34567u, 0x45678u 169 169 }; 170 170 171 171 initState(seedArray, seedArray + 4); 172 172 } … … 176 176 static const Word mul = RandomTraits<Word>::mul; 177 177 178 current = state; 178 current = state; 179 179 180 180 Word *curr = state + length  1; … … 202 202 num = length > end  begin ? length : end  begin; 203 203 while (num) { 204 curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits  2))) * mul1)) 204 curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits  2))) * mul1)) 205 205 + *it + cnt; 206 206 ++it; ++cnt; … … 224 224 } 225 225 } 226 226 227 227 state[length  1] = Word(1) << (bits  1); 228 228 } 229 229 230 230 void copyState(const RandomCore& other) { 231 231 std::copy(other.state, other.state + length, state); … … 242 242 private: 243 243 244 244 245 245 void fillState() { 246 246 static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask }; … … 248 248 static const Word hiMask = RandomTraits<Word>::hiMask; 249 249 250 current = state + length; 250 current = state + length; 251 251 252 252 register Word *curr = state + length  1; 253 253 register long num; 254 254 255 255 num = length  shift; 256 256 while (num) { … … 270 270 } 271 271 272 272 273 273 Word *current; 274 274 Word state[length]; 275 276 }; 277 278 279 template <typename Result, 275 276 }; 277 278 279 template <typename Result, 280 280 int shift = (std::numeric_limits<Result>::digits + 1) / 2> 281 281 struct Masker { … … 285 285 } 286 286 }; 287 287 288 288 template <typename Result> 289 289 struct Masker<Result, 1> { … … 293 293 }; 294 294 295 template <typename Result, typename Word, 296 int rest = std::numeric_limits<Result>::digits, int shift = 0, 295 template <typename Result, typename Word, 296 int rest = std::numeric_limits<Result>::digits, int shift = 0, 297 297 bool last = rest <= std::numeric_limits<Word>::digits> 298 298 struct IntConversion { 299 299 static const int bits = std::numeric_limits<Word>::digits; 300 300 301 301 static Result convert(RandomCore<Word>& rnd) { 302 302 return static_cast<Result>(rnd() >> (bits  rest)) << shift; 303 303 } 304 305 }; 306 307 template <typename Result, typename Word, int rest, int shift> 304 305 }; 306 307 template <typename Result, typename Word, int rest, int shift> 308 308 struct IntConversion<Result, Word, rest, shift, false> { 309 309 static const int bits = std::numeric_limits<Word>::digits; 310 310 311 311 static Result convert(RandomCore<Word>& rnd) { 312 return (static_cast<Result>(rnd()) << shift)  312 return (static_cast<Result>(rnd()) << shift)  313 313 IntConversion<Result, Word, rest  bits, shift + bits>::convert(rnd); 314 314 } … … 317 317 318 318 template <typename Result, typename Word, 319 bool one_word = (std::numeric_limits<Word>::digits < 320 319 bool one_word = (std::numeric_limits<Word>::digits < 320 std::numeric_limits<Result>::digits) > 321 321 struct Mapping { 322 322 static Result map(RandomCore<Word>& rnd, const Result& bound) { … … 325 325 Result num; 326 326 do { 327 num = IntConversion<Result, Word>::convert(rnd) & mask; 327 num = IntConversion<Result, Word>::convert(rnd) & mask; 328 328 } while (num > max); 329 329 return num; … … 351 351 res *= res; 352 352 if ((exp & 1) == 1) res *= static_cast<Result>(2.0); 353 return res; 353 return res; 354 354 } 355 355 }; … … 361 361 res *= res; 362 362 if ((exp & 1) == 1) res *= static_cast<Result>(0.5); 363 return res; 363 return res; 364 364 } 365 365 }; … … 368 368 struct ShiftMultiplier<Result, 0, true> { 369 369 static const Result multiplier() { 370 return static_cast<Result>(1.0); 370 return static_cast<Result>(1.0); 371 371 } 372 372 }; … … 375 375 struct ShiftMultiplier<Result, 20, true> { 376 376 static const Result multiplier() { 377 return static_cast<Result>(1.0/1048576.0); 378 } 379 }; 380 377 return static_cast<Result>(1.0/1048576.0); 378 } 379 }; 380 381 381 template <typename Result> 382 382 struct ShiftMultiplier<Result, 32, true> { 383 383 static const Result multiplier() { 384 return static_cast<Result>(1.0/424967296.0); 384 return static_cast<Result>(1.0/424967296.0); 385 385 } 386 386 }; … … 389 389 struct ShiftMultiplier<Result, 53, true> { 390 390 static const Result multiplier() { 391 return static_cast<Result>(1.0/9007199254740992.0); 391 return static_cast<Result>(1.0/9007199254740992.0); 392 392 } 393 393 }; … … 396 396 struct ShiftMultiplier<Result, 64, true> { 397 397 static const Result multiplier() { 398 return static_cast<Result>(1.0/18446744073709551616.0); 398 return static_cast<Result>(1.0/18446744073709551616.0); 399 399 } 400 400 }; … … 408 408 409 409 template <typename Result, typename Word, 410 int rest = std::numeric_limits<Result>::digits, int shift = 0, 410 int rest = std::numeric_limits<Result>::digits, int shift = 0, 411 411 bool last = rest <= std::numeric_limits<Word>::digits> 412 struct RealConversion{ 412 struct RealConversion{ 413 413 static const int bits = std::numeric_limits<Word>::digits; 414 414 … … 420 420 421 421 template <typename Result, typename Word, int rest, int shift> 422 struct RealConversion<Result, Word, rest, shift, false> { 422 struct RealConversion<Result, Word, rest, shift, false> { 423 423 static const int bits = std::numeric_limits<Word>::digits; 424 424 … … 459 459 Word buffer; 460 460 int num; 461 461 462 462 BoolProducer() : num(0) {} 463 463 … … 530 530 // Architecture word 531 531 typedef unsigned long Word; 532 532 533 533 _random_bits::RandomCore<Word> core; 534 534 _random_bits::BoolProducer<Word> bool_producer; 535 535 536 536 537 537 public: … … 555 555 /// to the architecture word type. 556 556 template <typename Number> 557 Random(Number seed) { 557 Random(Number seed) { 558 558 _random_bits::Initializer<Number, Word>::init(core, seed); 559 559 } … … 565 565 /// architecture word type. 566 566 template <typename Iterator> 567 Random(Iterator begin, Iterator end) { 567 Random(Iterator begin, Iterator end) { 568 568 typedef typename std::iterator_traits<Iterator>::value_type Number; 569 569 _random_bits::Initializer<Number, Word>::init(core, begin, end); … … 598 598 /// converted to the architecture word type. 599 599 template <typename Number> 600 void seed(Number seed) { 600 void seed(Number seed) { 601 601 _random_bits::Initializer<Number, Word>::init(core, seed); 602 602 } … … 608 608 /// architecture word type. 609 609 template <typename Iterator> 610 void seed(Iterator begin, Iterator end) { 610 void seed(Iterator begin, Iterator end) { 611 611 typedef typename std::iterator_traits<Iterator>::value_type Number; 612 612 _random_bits::Initializer<Number, Word>::init(core, begin, end); … … 626 626 return false; 627 627 } 628 628 629 629 /// \brief Seeding from file 630 630 /// … … 641 641 /// \return True when the seeding successes. 642 642 #ifndef WIN32 643 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 643 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 644 644 #else 645 bool seedFromFile(const std::string& file = "", int offset = 0) 645 bool seedFromFile(const std::string& file = "", int offset = 0) 646 646 #endif 647 647 { … … 661 661 /// random sequence. 662 662 /// \return Currently always true. 663 bool seedFromTime() { 663 bool seedFromTime() { 664 664 #ifndef WIN32 665 665 timeval tv; … … 697 697 /// It returns a random real number from the range [0, b). 698 698 template <typename Number> 699 Number real(Number b) { 700 return real<Number>() * b; 699 Number real(Number b) { 700 return real<Number>() * b; 701 701 } 702 702 … … 705 705 /// It returns a random real number from the range [a, b). 706 706 template <typename Number> 707 Number real(Number a, Number b) { 708 return real<Number>() * (b  a) + a; 707 Number real(Number a, Number b) { 708 return real<Number>() * (b  a) + a; 709 709 } 710 710 … … 726 726 /// It returns a random real number from the range [0, b). 727 727 template <typename Number> 728 Number operator()(Number b) { 729 return real<Number>() * b; 728 Number operator()(Number b) { 729 return real<Number>() * b; 730 730 } 731 731 … … 734 734 /// It returns a random real number from the range [a, b). 735 735 template <typename Number> 736 Number operator()(Number a, Number b) { 737 return real<Number>() * (b  a) + a; 736 Number operator()(Number a, Number b) { 737 return real<Number>() * (b  a) + a; 738 738 } 739 739 … … 785 785 template <typename Number> 786 786 Number integer() { 787 static const int nb = std::numeric_limits<Number>::digits + 787 static const int nb = std::numeric_limits<Number>::digits + 788 788 (std::numeric_limits<Number>::is_signed ? 1 : 0); 789 789 return _random_bits::IntConversion<Number, Word, nb>::convert(core); … … 793 793 return integer<int>(); 794 794 } 795 795 796 796 /// \brief Returns a random bool 797 797 /// … … 807 807 ///\name Nonuniform distributions 808 808 /// 809 809 810 810 ///@{ 811 811 812 812 /// \brief Returns a random bool 813 813 /// … … 823 823 /// transformation is used to generate a random normal distribution. 824 824 /// \todo Consider using the "ziggurat" method instead. 825 double gauss() 825 double gauss() 826 826 { 827 827 double V1,V2,S; 828 828 do { 829 830 831 829 V1=2*real<double>()1; 830 V2=2*real<double>()1; 831 S=V1*V1+V2*V2; 832 832 } while(S>=1); 833 833 return std::sqrt(2*std::log(S)/S)*V1; … … 855 855 856 856 /// This function generates a gamma distribution random number. 857 /// 857 /// 858 858 ///\param k shape parameter (<tt>k>0</tt> integer) 859 double gamma(int k) 859 double gamma(int k) 860 860 { 861 861 double s = 0; … … 863 863 return s; 864 864 } 865 865 866 866 /// Gamma distribution with given shape and scale parameter 867 867 868 868 /// This function generates a gamma distribution random number. 869 /// 869 /// 870 870 ///\param k shape parameter (<tt>k>0</tt>) 871 871 ///\param theta scale parameter … … 877 877 const double v0=E/(Edelta); 878 878 do { 879 880 881 882 if(V2<=v0) 883 884 885 886 887 else 888 889 890 891 879 double V0=1.0real<double>(); 880 double V1=1.0real<double>(); 881 double V2=1.0real<double>(); 882 if(V2<=v0) 883 { 884 xi=std::pow(V1,1.0/delta); 885 nu=V0*std::pow(xi,delta1.0); 886 } 887 else 888 { 889 xi=1.0std::log(V1); 890 nu=V0*std::exp(xi); 891 } 892 892 } while(nu>std::pow(xi,delta1.0)*std::exp(xi)); 893 893 return theta*(xi+gamma(int(std::floor(k)))); 894 894 } 895 895 896 896 /// Weibull distribution 897 897 898 898 /// This function generates a Weibull distribution random number. 899 /// 899 /// 900 900 ///\param k shape parameter (<tt>k>0</tt>) 901 901 ///\param lambda scale parameter (<tt>lambda>0</tt>) … … 904 904 { 905 905 return lambda*pow(std::log(1.0real<double>()),1.0/k); 906 } 907 906 } 907 908 908 /// Pareto distribution 909 909 910 910 /// This function generates a Pareto distribution random number. 911 /// 911 /// 912 912 ///\param k shape parameter (<tt>k>0</tt>) 913 913 ///\param x_min location parameter (<tt>x_min>0</tt>) … … 916 916 { 917 917 return exponential(gamma(k,1.0/x_min))+x_min; 918 } 919 918 } 919 920 920 /// Poisson distribution 921 921 922 922 /// This function generates a Poisson distribution random number with 923 923 /// parameter \c lambda. 924 /// 924 /// 925 925 /// The probability mass function of this distribusion is 926 926 /// \f[ \frac{e^{\lambda}\lambda^k}{k!} \f] … … 928 928 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the 929 929 /// return value. 930 930 931 931 int poisson(double lambda) 932 932 { … … 935 935 double p = 1.0; 936 936 do { 937 938 937 k++; 938 p*=real<double>(); 939 939 } while (p>=l); 940 940 return k1; 941 } 942 941 } 942 943 943 ///@} 944 944 945 945 ///\name Two dimensional distributions 946 946 /// 947 947 948 948 ///@{ 949 949 950 950 /// Uniform distribution on the full unit circle 951 951 952 952 /// Uniform distribution on the full unit circle. 953 953 /// 954 dim2::Point<double> disc() 954 dim2::Point<double> disc() 955 955 { 956 956 double V1,V2; 957 957 do { 958 959 960 958 V1=2*real<double>()1; 959 V2=2*real<double>()1; 960 961 961 } while(V1*V1+V2*V2>=1); 962 962 return dim2::Point<double>(V1,V2); … … 974 974 double V1,V2,S; 975 975 do { 976 977 978 976 V1=2*real<double>()1; 977 V2=2*real<double>()1; 978 S=V1*V1+V2*V2; 979 979 } while(S>=1); 980 980 double W=std::sqrt(2*std::log(S)/S); … … 985 985 /// This function provides a turning symmetric twodimensional distribution. 986 986 /// The xcoordinate is of conditionally exponential distribution 987 /// with the condition that x is positive and y=0. If x is negative and 987 /// with the condition that x is positive and y=0. If x is negative and 988 988 /// y=0 then, x is of exponential distribution. The same is true for the 989 989 /// ycoordinate. 990 dim2::Point<double> exponential2() 990 dim2::Point<double> exponential2() 991 991 { 992 992 double V1,V2,S; 993 993 do { 994 995 996 994 V1=2*real<double>()1; 995 V2=2*real<double>()1; 996 S=V1*V1+V2*V2; 997 997 } while(S>=1); 998 998 double W=std::log(S)/S; … … 1000 1000 } 1001 1001 1002 ///@} 1002 ///@} 1003 1003 }; 1004 1004
Note: See TracChangeset
for help on using the changeset viewer.