Changeset 1379:db1d342a1087 in lemon
- Timestamp:
- 10/08/15 13:48:09 (9 years ago)
- Branch:
- default
- Phase:
- public
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/random.h
r1343 r1379 297 297 template <typename Result, typename Word, 298 298 int rest = std::numeric_limits<Result>::digits, int shift = 0, 299 bool last = rest <= std::numeric_limits<Word>::digits>299 bool last = (rest <= std::numeric_limits<Word>::digits)> 300 300 struct IntConversion { 301 301 static const int bits = std::numeric_limits<Word>::digits; … … 466 466 }; 467 467 468 } 468 /// \ingroup misc 469 /// 470 /// \brief Mersenne Twister random number generator 471 /// 472 /// The Mersenne Twister is a twisted generalized feedback 473 /// shift-register generator of Matsumoto and Nishimura. The period 474 /// of this generator is \f$ 2^{19937} - 1 \f$ and it is 475 /// equi-distributed in 623 dimensions for 32-bit numbers. The time 476 /// performance of this generator is comparable to the commonly used 477 /// generators. 478 /// 479 /// This is a template version implementation both 32-bit and 480 /// 64-bit architecture optimized versions. The generators differ 481 /// sligthly in the initialization and generation phase so they 482 /// produce two completly different sequences. 483 /// 484 /// \alert Do not use this class directly, but instead one of \ref 485 /// Random, \ref Random32 or \ref Random64. 486 /// 487 /// The generator gives back random numbers of serveral types. To 488 /// get a random number from a range of a floating point type you 489 /// can use one form of the \c operator() or the \c real() member 490 /// function. If you want to get random number from the {0, 1, ..., 491 /// n-1} integer range use the \c operator[] or the \c integer() 492 /// method. And to get random number from the whole range of an 493 /// integer type you can use the argumentless \c integer() or \c 494 /// uinteger() functions. After all you can get random bool with 495 /// equal chance of true and false or given probability of true 496 /// result with the \c boolean() member functions. 497 /// 498 ///\code 499 /// // The commented code is identical to the other 500 /// double a = rnd(); // [0.0, 1.0) 501 /// // double a = rnd.real(); // [0.0, 1.0) 502 /// double b = rnd(100.0); // [0.0, 100.0) 503 /// // double b = rnd.real(100.0); // [0.0, 100.0) 504 /// double c = rnd(1.0, 2.0); // [1.0, 2.0) 505 /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) 506 /// int d = rnd[100000]; // 0..99999 507 /// // int d = rnd.integer(100000); // 0..99999 508 /// int e = rnd[6] + 1; // 1..6 509 /// // int e = rnd.integer(1, 1 + 6); // 1..6 510 /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1 511 /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1 512 /// bool g = rnd.boolean(); // P(g = true) = 0.5 513 /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 514 ///\endcode 515 /// 516 /// LEMON provides a global instance of the random number 517 /// generator which name is \ref lemon::rnd "rnd". Usually it is a 518 /// good programming convenience to use this global generator to get 519 /// random numbers. 520 /// 521 /// \sa \ref Random, \ref Random32 or \ref Random64. 522 /// 523 template<class Word> 524 class Random { 525 private: 526 527 _random_bits::RandomCore<Word> core; 528 _random_bits::BoolProducer<Word> bool_producer; 529 530 531 public: 532 533 ///\name Initialization 534 /// 535 /// @{ 536 537 /// \brief Default constructor 538 /// 539 /// Constructor with constant seeding. 540 Random() { core.initState(); } 541 542 /// \brief Constructor with seed 543 /// 544 /// Constructor with seed. The current number type will be converted 545 /// to the architecture word type. 546 template <typename Number> 547 Random(Number seed) { 548 _random_bits::Initializer<Number, Word>::init(core, seed); 549 } 550 551 /// \brief Constructor with array seeding 552 /// 553 /// Constructor with array seeding. The given range should contain 554 /// any number type and the numbers will be converted to the 555 /// architecture word type. 556 template <typename Iterator> 557 Random(Iterator begin, Iterator end) { 558 typedef typename std::iterator_traits<Iterator>::value_type Number; 559 _random_bits::Initializer<Number, Word>::init(core, begin, end); 560 } 561 562 /// \brief Copy constructor 563 /// 564 /// Copy constructor. The generated sequence will be identical to 565 /// the other sequence. It can be used to save the current state 566 /// of the generator and later use it to generate the same 567 /// sequence. 568 Random(const Random& other) { 569 core.copyState(other.core); 570 } 571 572 /// \brief Assign operator 573 /// 574 /// Assign operator. The generated sequence will be identical to 575 /// the other sequence. It can be used to save the current state 576 /// of the generator and later use it to generate the same 577 /// sequence. 578 Random& operator=(const Random& other) { 579 if (&other != this) { 580 core.copyState(other.core); 581 } 582 return *this; 583 } 584 585 /// \brief Seeding random sequence 586 /// 587 /// Seeding the random sequence. The current number type will be 588 /// converted to the architecture word type. 589 template <typename Number> 590 void seed(Number seed) { 591 _random_bits::Initializer<Number, Word>::init(core, seed); 592 } 593 594 /// \brief Seeding random sequence 595 /// 596 /// Seeding the random sequence. The given range should contain 597 /// any number type and the numbers will be converted to the 598 /// architecture word type. 599 template <typename Iterator> 600 void seed(Iterator begin, Iterator end) { 601 typedef typename std::iterator_traits<Iterator>::value_type Number; 602 _random_bits::Initializer<Number, Word>::init(core, begin, end); 603 } 604 605 /// \brief Seeding from file or from process id and time 606 /// 607 /// By default, this function calls the \c seedFromFile() member 608 /// function with the <tt>/dev/urandom</tt> file. If it does not success, 609 /// it uses the \c seedFromTime(). 610 /// \return Currently always \c true. 611 bool seed() { 612 #ifndef LEMON_WIN32 613 if (seedFromFile("/dev/urandom", 0)) return true; 614 #endif 615 if (seedFromTime()) return true; 616 return false; 617 } 618 619 /// \brief Seeding from file 620 /// 621 /// Seeding the random sequence from file. The linux kernel has two 622 /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which 623 /// could give good seed values for pseudo random generators (The 624 /// difference between two devices is that the <tt>random</tt> may 625 /// block the reading operation while the kernel can give good 626 /// source of randomness, while the <tt>urandom</tt> does not 627 /// block the input, but it could give back bytes with worse 628 /// entropy). 629 /// \param file The source file 630 /// \param offset The offset, from the file read. 631 /// \return \c true when the seeding successes. 632 #ifndef LEMON_WIN32 633 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 634 #else 635 bool seedFromFile(const std::string& file = "", int offset = 0) 636 #endif 637 { 638 std::ifstream rs(file.c_str()); 639 const int size = 4; 640 Word buf[size]; 641 if (offset != 0 && !rs.seekg(offset)) return false; 642 if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false; 643 seed(buf, buf + size); 644 return true; 645 } 646 647 /// \brief Seding from process id and time 648 /// 649 /// Seding from process id and time. This function uses the 650 /// current process id and the current time for initialize the 651 /// random sequence. 652 /// \return Currently always \c true. 653 bool seedFromTime() { 654 #ifndef LEMON_WIN32 655 timeval tv; 656 gettimeofday(&tv, 0); 657 seed(getpid() + tv.tv_sec + tv.tv_usec); 658 #else 659 seed(bits::getWinRndSeed()); 660 #endif 661 return true; 662 } 663 664 /// @} 665 666 ///\name Uniform Distributions 667 /// 668 /// @{ 669 670 /// \brief Returns a random real number from the range [0, 1) 671 /// 672 /// It returns a random real number from the range [0, 1). The 673 /// default Number type is \c double. 674 template <typename Number> 675 Number real() { 676 return _random_bits::RealConversion<Number, Word>::convert(core); 677 } 678 679 double real() { 680 return real<double>(); 681 } 682 683 /// \brief Returns a random real number from the range [0, 1) 684 /// 685 /// It returns a random double from the range [0, 1). 686 double operator()() { 687 return real<double>(); 688 } 689 690 /// \brief Returns a random real number from the range [0, b) 691 /// 692 /// It returns a random real number from the range [0, b). 693 double operator()(double b) { 694 return real<double>() * b; 695 } 696 697 /// \brief Returns a random real number from the range [a, b) 698 /// 699 /// It returns a random real number from the range [a, b). 700 double operator()(double a, double b) { 701 return real<double>() * (b - a) + a; 702 } 703 704 /// \brief Returns a random integer from a range 705 /// 706 /// It returns a random integer from the range {0, 1, ..., b - 1}. 707 template <typename Number> 708 Number integer(Number b) { 709 return _random_bits::Mapping<Number, Word>::map(core, b); 710 } 711 712 /// \brief Returns a random integer from a range 713 /// 714 /// It returns a random integer from the range {a, a + 1, ..., b - 1}. 715 template <typename Number> 716 Number integer(Number a, Number b) { 717 return _random_bits::Mapping<Number, Word>::map(core, b - a) + a; 718 } 719 720 /// \brief Returns a random integer from a range 721 /// 722 /// It returns a random integer from the range {0, 1, ..., b - 1}. 723 template <typename Number> 724 Number operator[](Number b) { 725 return _random_bits::Mapping<Number, Word>::map(core, b); 726 } 727 728 /// \brief Returns a random non-negative integer 729 /// 730 /// It returns a random non-negative integer uniformly from the 731 /// whole range of the current \c Number type. The default result 732 /// type of this function is <tt>unsigned int</tt>. 733 template <typename Number> 734 Number uinteger() { 735 return _random_bits::IntConversion<Number, Word>::convert(core); 736 } 737 738 unsigned int uinteger() { 739 return uinteger<unsigned int>(); 740 } 741 742 /// \brief Returns a random integer 743 /// 744 /// It returns a random integer uniformly from the whole range of 745 /// the current \c Number type. The default result type of this 746 /// function is \c int. 747 template <typename Number> 748 Number integer() { 749 static const int nb = std::numeric_limits<Number>::digits + 750 (std::numeric_limits<Number>::is_signed ? 1 : 0); 751 return _random_bits::IntConversion<Number, Word, nb>::convert(core); 752 } 753 754 int integer() { 755 return integer<int>(); 756 } 757 758 /// \brief Returns a random bool 759 /// 760 /// It returns a random bool. The generator holds a buffer for 761 /// random bits. Every time when it become empty the generator makes 762 /// a new random word and fill the buffer up. 763 bool boolean() { 764 return bool_producer.convert(core); 765 } 766 767 /// @} 768 769 ///\name Non-uniform Distributions 770 /// 771 ///@{ 772 773 /// \brief Returns a random bool with given probability of true result. 774 /// 775 /// It returns a random bool with given probability of true result. 776 bool boolean(double p) { 777 return operator()() < p; 778 } 779 780 /// Standard normal (Gauss) distribution 781 782 /// Standard normal (Gauss) distribution. 783 /// \note The Cartesian form of the Box-Muller 784 /// transformation is used to generate a random normal distribution. 785 double gauss() 786 { 787 double V1,V2,S; 788 do { 789 V1=2*real<double>()-1; 790 V2=2*real<double>()-1; 791 S=V1*V1+V2*V2; 792 } while(S>=1); 793 return std::sqrt(-2*std::log(S)/S)*V1; 794 } 795 /// Normal (Gauss) distribution with given mean and standard deviation 796 797 /// Normal (Gauss) distribution with given mean and standard deviation. 798 /// \sa gauss() 799 double gauss(double mean,double std_dev) 800 { 801 return gauss()*std_dev+mean; 802 } 803 804 /// Lognormal distribution 805 806 /// Lognormal distribution. The parameters are the mean and the standard 807 /// deviation of <tt>exp(X)</tt>. 808 /// 809 double lognormal(double n_mean,double n_std_dev) 810 { 811 return std::exp(gauss(n_mean,n_std_dev)); 812 } 813 /// Lognormal distribution 814 815 /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of 816 /// the mean and the standard deviation of <tt>exp(X)</tt>. 817 /// 818 double lognormal(const std::pair<double,double> ¶ms) 819 { 820 return std::exp(gauss(params.first,params.second)); 821 } 822 /// Compute the lognormal parameters from mean and standard deviation 823 824 /// This function computes the lognormal parameters from mean and 825 /// standard deviation. The return value can direcly be passed to 826 /// lognormal(). 827 std::pair<double,double> lognormalParamsFromMD(double mean, 828 double std_dev) 829 { 830 double fr=std_dev/mean; 831 fr*=fr; 832 double lg=std::log(1+fr); 833 return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg)); 834 } 835 /// Lognormal distribution with given mean and standard deviation 836 837 /// Lognormal distribution with given mean and standard deviation. 838 /// 839 double lognormalMD(double mean,double std_dev) 840 { 841 return lognormal(lognormalParamsFromMD(mean,std_dev)); 842 } 843 844 /// Exponential distribution with given mean 845 846 /// This function generates an exponential distribution random number 847 /// with mean <tt>1/lambda</tt>. 848 /// 849 double exponential(double lambda=1.0) 850 { 851 return -std::log(1.0-real<double>())/lambda; 852 } 853 854 /// Gamma distribution with given integer shape 855 856 /// This function generates a gamma distribution random number. 857 /// 858 ///\param k shape parameter (<tt>k>0</tt> integer) 859 double gamma(int k) 860 { 861 double s = 0; 862 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); 863 return s; 864 } 865 866 /// Gamma distribution with given shape and scale parameter 867 868 /// This function generates a gamma distribution random number. 869 /// 870 ///\param k shape parameter (<tt>k>0</tt>) 871 ///\param theta scale parameter 872 /// 873 double gamma(double k,double theta=1.0) 874 { 875 double xi,nu; 876 const double delta = k-std::floor(k); 877 const double v0=E/(E-delta); 878 do { 879 double V0=1.0-real<double>(); 880 double V1=1.0-real<double>(); 881 double V2=1.0-real<double>(); 882 if(V2<=v0) 883 { 884 xi=std::pow(V1,1.0/delta); 885 nu=V0*std::pow(xi,delta-1.0); 886 } 887 else 888 { 889 xi=1.0-std::log(V1); 890 nu=V0*std::exp(-xi); 891 } 892 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); 893 return theta*(xi+gamma(int(std::floor(k)))); 894 } 895 896 /// Weibull distribution 897 898 /// This function generates a Weibull distribution random number. 899 /// 900 ///\param k shape parameter (<tt>k>0</tt>) 901 ///\param lambda scale parameter (<tt>lambda>0</tt>) 902 /// 903 double weibull(double k,double lambda) 904 { 905 return lambda*pow(-std::log(1.0-real<double>()),1.0/k); 906 } 907 908 /// Pareto distribution 909 910 /// This function generates a Pareto distribution random number. 911 /// 912 ///\param k shape parameter (<tt>k>0</tt>) 913 ///\param x_min location parameter (<tt>x_min>0</tt>) 914 /// 915 double pareto(double k,double x_min) 916 { 917 return exponential(gamma(k,1.0/x_min))+x_min; 918 } 919 920 /// Poisson distribution 921 922 /// This function generates a Poisson distribution random number with 923 /// parameter \c lambda. 924 /// 925 /// The probability mass function of this distribusion is 926 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] 927 /// \note The algorithm is taken from the book of Donald E. Knuth titled 928 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the 929 /// return value. 930 931 int poisson(double lambda) 932 { 933 const double l = std::exp(-lambda); 934 int k=0; 935 double p = 1.0; 936 do { 937 k++; 938 p*=real<double>(); 939 } while (p>=l); 940 return k-1; 941 } 942 943 ///@} 944 945 ///\name Two Dimensional Distributions 946 /// 947 ///@{ 948 949 /// Uniform distribution on the full unit circle 950 951 /// Uniform distribution on the full unit circle. 952 /// 953 dim2::Point<double> disc() 954 { 955 double V1,V2; 956 do { 957 V1=2*real<double>()-1; 958 V2=2*real<double>()-1; 959 960 } while(V1*V1+V2*V2>=1); 961 return dim2::Point<double>(V1,V2); 962 } 963 /// A kind of two dimensional normal (Gauss) distribution 964 965 /// This function provides a turning symmetric two-dimensional distribution. 966 /// Both coordinates are of standard normal distribution, but they are not 967 /// independent. 968 /// 969 /// \note The coordinates are the two random variables provided by 970 /// the Box-Muller method. 971 dim2::Point<double> gauss2() 972 { 973 double V1,V2,S; 974 do { 975 V1=2*real<double>()-1; 976 V2=2*real<double>()-1; 977 S=V1*V1+V2*V2; 978 } while(S>=1); 979 double W=std::sqrt(-2*std::log(S)/S); 980 return dim2::Point<double>(W*V1,W*V2); 981 } 982 /// A kind of two dimensional exponential distribution 983 984 /// This function provides a turning symmetric two-dimensional distribution. 985 /// The x-coordinate is of conditionally exponential distribution 986 /// with the condition that x is positive and y=0. If x is negative and 987 /// y=0 then, -x is of exponential distribution. The same is true for the 988 /// y-coordinate. 989 dim2::Point<double> exponential2() 990 { 991 double V1,V2,S; 992 do { 993 V1=2*real<double>()-1; 994 V2=2*real<double>()-1; 995 S=V1*V1+V2*V2; 996 } while(S>=1); 997 double W=-std::log(S)/S; 998 return dim2::Point<double>(W*V1,W*V2); 999 } 1000 1001 ///@} 1002 }; 1003 1004 1005 }; 469 1006 470 1007 /// \ingroup misc … … 472 1009 /// \brief Mersenne Twister random number generator 473 1010 /// 474 /// Th e Mersenne Twister is a twisted generalized feedback475 /// shift-register generator of Matsumoto and Nishimura. The period476 /// of this generator is \f$ 2^{19937} - 1 \f$ and it is477 /// equi-distributed in 623 dimensions for 32-bit numbers. The time478 /// performance of this generator is comparable to the commonly used479 /// generators.1011 /// This class implements either the 32 bit or the 64 bit version of 1012 /// the Mersenne Twister random number generator algorithm 1013 /// depending the the system architecture. 1014 /// 1015 /// For the API description, see its base class \ref 1016 /// _random_bits::Random 480 1017 /// 481 /// This implementation is specialized for both 32-bit and 64-bit 482 /// architectures. The generators differ sligthly in the 483 /// initialization and generation phase so they produce two 484 /// completly different sequences. 1018 /// \sa \ref _random_bits::Random 1019 typedef _random_bits::Random<unsigned long> Random; 1020 /// \ingroup misc 485 1021 /// 486 /// The generator gives back random numbers of serveral types. To 487 /// get a random number from a range of a floating point type you 488 /// can use one form of the \c operator() or the \c real() member 489 /// function. If you want to get random number from the {0, 1, ..., 490 /// n-1} integer range use the \c operator[] or the \c integer() 491 /// method. And to get random number from the whole range of an 492 /// integer type you can use the argumentless \c integer() or \c 493 /// uinteger() functions. After all you can get random bool with 494 /// equal chance of true and false or given probability of true 495 /// result with the \c boolean() member functions. 1022 /// \brief Mersenne Twister random number generator (32 bit version) 496 1023 /// 497 ///\code 498 /// // The commented code is identical to the other 499 /// double a = rnd(); // [0.0, 1.0) 500 /// // double a = rnd.real(); // [0.0, 1.0) 501 /// double b = rnd(100.0); // [0.0, 100.0) 502 /// // double b = rnd.real(100.0); // [0.0, 100.0) 503 /// double c = rnd(1.0, 2.0); // [1.0, 2.0) 504 /// // double c = rnd.real(1.0, 2.0); // [1.0, 2.0) 505 /// int d = rnd[100000]; // 0..99999 506 /// // int d = rnd.integer(100000); // 0..99999 507 /// int e = rnd[6] + 1; // 1..6 508 /// // int e = rnd.integer(1, 1 + 6); // 1..6 509 /// int b = rnd.uinteger<int>(); // 0 .. 2^31 - 1 510 /// int c = rnd.integer<int>(); // - 2^31 .. 2^31 - 1 511 /// bool g = rnd.boolean(); // P(g = true) = 0.5 512 /// bool h = rnd.boolean(0.8); // P(h = true) = 0.8 513 ///\endcode 1024 /// This class implements the 32 bit version of the Mersenne Twister 1025 /// random number generator algorithm. It is recommended to be used 1026 /// when someone wants to make sure that the \e same pseudo random 1027 /// sequence will be generated on every platfrom. 514 1028 /// 515 /// LEMON provides a global instance of the random number 516 /// generator which name is \ref lemon::rnd "rnd". Usually it is a 517 /// good programming convenience to use this global generator to get 518 /// random numbers. 519 class Random { 520 private: 521 522 // Architecture word 523 typedef unsigned long Word; 524 525 _random_bits::RandomCore<Word> core; 526 _random_bits::BoolProducer<Word> bool_producer; 527 528 529 public: 530 531 ///\name Initialization 532 /// 533 /// @{ 534 535 /// \brief Default constructor 536 /// 537 /// Constructor with constant seeding. 538 Random() { core.initState(); } 539 540 /// \brief Constructor with seed 541 /// 542 /// Constructor with seed. The current number type will be converted 543 /// to the architecture word type. 544 template <typename Number> 545 Random(Number seed) { 546 _random_bits::Initializer<Number, Word>::init(core, seed); 547 } 548 549 /// \brief Constructor with array seeding 550 /// 551 /// Constructor with array seeding. The given range should contain 552 /// any number type and the numbers will be converted to the 553 /// architecture word type. 554 template <typename Iterator> 555 Random(Iterator begin, Iterator end) { 556 typedef typename std::iterator_traits<Iterator>::value_type Number; 557 _random_bits::Initializer<Number, Word>::init(core, begin, end); 558 } 559 560 /// \brief Copy constructor 561 /// 562 /// Copy constructor. The generated sequence will be identical to 563 /// the other sequence. It can be used to save the current state 564 /// of the generator and later use it to generate the same 565 /// sequence. 566 Random(const Random& other) { 567 core.copyState(other.core); 568 } 569 570 /// \brief Assign operator 571 /// 572 /// Assign operator. The generated sequence will be identical to 573 /// the other sequence. It can be used to save the current state 574 /// of the generator and later use it to generate the same 575 /// sequence. 576 Random& operator=(const Random& other) { 577 if (&other != this) { 578 core.copyState(other.core); 579 } 580 return *this; 581 } 582 583 /// \brief Seeding random sequence 584 /// 585 /// Seeding the random sequence. The current number type will be 586 /// converted to the architecture word type. 587 template <typename Number> 588 void seed(Number seed) { 589 _random_bits::Initializer<Number, Word>::init(core, seed); 590 } 591 592 /// \brief Seeding random sequence 593 /// 594 /// Seeding the random sequence. The given range should contain 595 /// any number type and the numbers will be converted to the 596 /// architecture word type. 597 template <typename Iterator> 598 void seed(Iterator begin, Iterator end) { 599 typedef typename std::iterator_traits<Iterator>::value_type Number; 600 _random_bits::Initializer<Number, Word>::init(core, begin, end); 601 } 602 603 /// \brief Seeding from file or from process id and time 604 /// 605 /// By default, this function calls the \c seedFromFile() member 606 /// function with the <tt>/dev/urandom</tt> file. If it does not success, 607 /// it uses the \c seedFromTime(). 608 /// \return Currently always \c true. 609 bool seed() { 610 #ifndef LEMON_WIN32 611 if (seedFromFile("/dev/urandom", 0)) return true; 1029 /// For the API description, see its base class \ref 1030 /// _random_bits::Random 1031 /// 1032 /// \sa \ref _random_bits::Random 1033 1034 typedef _random_bits::Random<unsigned int> Random32; 1035 /// \ingroup misc 1036 /// 1037 /// \brief Mersenne Twister random number generator (64 bit version) 1038 /// 1039 /// This class implements the 64 bit version of the Mersenne Twister 1040 /// random number generator algorithm. (Even though it will run 1041 /// on 32 bit architectures, too.) It is recommended to ber used when 1042 /// someone wants to make sure that the \e same pseudo random sequence 1043 /// will be generated on every platfrom. 1044 /// 1045 /// For the API description, see its base class \ref 1046 /// _random_bits::Random 1047 /// 1048 /// \sa \ref _random_bits::Random 1049 typedef _random_bits::Random<unsigned long long> Random64; 1050 1051 1052 extern Random rnd; 1053 1054 1055 } 1056 612 1057 #endif 613 if (seedFromTime()) return true;614 return false;615 }616 617 /// \brief Seeding from file618 ///619 /// Seeding the random sequence from file. The linux kernel has two620 /// devices, <tt>/dev/random</tt> and <tt>/dev/urandom</tt> which621 /// could give good seed values for pseudo random generators (The622 /// difference between two devices is that the <tt>random</tt> may623 /// block the reading operation while the kernel can give good624 /// source of randomness, while the <tt>urandom</tt> does not625 /// block the input, but it could give back bytes with worse626 /// entropy).627 /// \param file The source file628 /// \param offset The offset, from the file read.629 /// \return \c true when the seeding successes.630 #ifndef LEMON_WIN32631 bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)632 #else633 bool seedFromFile(const std::string& file = "", int offset = 0)634 #endif635 {636 std::ifstream rs(file.c_str());637 const int size = 4;638 Word buf[size];639 if (offset != 0 && !rs.seekg(offset)) return false;640 if (!rs.read(reinterpret_cast<char*>(buf), sizeof(buf))) return false;641 seed(buf, buf + size);642 return true;643 }644 645 /// \brief Seding from process id and time646 ///647 /// Seding from process id and time. This function uses the648 /// current process id and the current time for initialize the649 /// random sequence.650 /// \return Currently always \c true.651 bool seedFromTime() {652 #ifndef LEMON_WIN32653 timeval tv;654 gettimeofday(&tv, 0);655 seed(getpid() + tv.tv_sec + tv.tv_usec);656 #else657 seed(bits::getWinRndSeed());658 #endif659 return true;660 }661 662 /// @}663 664 ///\name Uniform Distributions665 ///666 /// @{667 668 /// \brief Returns a random real number from the range [0, 1)669 ///670 /// It returns a random real number from the range [0, 1). The671 /// default Number type is \c double.672 template <typename Number>673 Number real() {674 return _random_bits::RealConversion<Number, Word>::convert(core);675 }676 677 double real() {678 return real<double>();679 }680 681 /// \brief Returns a random real number from the range [0, 1)682 ///683 /// It returns a random double from the range [0, 1).684 double operator()() {685 return real<double>();686 }687 688 /// \brief Returns a random real number from the range [0, b)689 ///690 /// It returns a random real number from the range [0, b).691 double operator()(double b) {692 return real<double>() * b;693 }694 695 /// \brief Returns a random real number from the range [a, b)696 ///697 /// It returns a random real number from the range [a, b).698 double operator()(double a, double b) {699 return real<double>() * (b - a) + a;700 }701 702 /// \brief Returns a random integer from a range703 ///704 /// It returns a random integer from the range {0, 1, ..., b - 1}.705 template <typename Number>706 Number integer(Number b) {707 return _random_bits::Mapping<Number, Word>::map(core, b);708 }709 710 /// \brief Returns a random integer from a range711 ///712 /// It returns a random integer from the range {a, a + 1, ..., b - 1}.713 template <typename Number>714 Number integer(Number a, Number b) {715 return _random_bits::Mapping<Number, Word>::map(core, b - a) + a;716 }717 718 /// \brief Returns a random integer from a range719 ///720 /// It returns a random integer from the range {0, 1, ..., b - 1}.721 template <typename Number>722 Number operator[](Number b) {723 return _random_bits::Mapping<Number, Word>::map(core, b);724 }725 726 /// \brief Returns a random non-negative integer727 ///728 /// It returns a random non-negative integer uniformly from the729 /// whole range of the current \c Number type. The default result730 /// type of this function is <tt>unsigned int</tt>.731 template <typename Number>732 Number uinteger() {733 return _random_bits::IntConversion<Number, Word>::convert(core);734 }735 736 unsigned int uinteger() {737 return uinteger<unsigned int>();738 }739 740 /// \brief Returns a random integer741 ///742 /// It returns a random integer uniformly from the whole range of743 /// the current \c Number type. The default result type of this744 /// function is \c int.745 template <typename Number>746 Number integer() {747 static const int nb = std::numeric_limits<Number>::digits +748 (std::numeric_limits<Number>::is_signed ? 1 : 0);749 return _random_bits::IntConversion<Number, Word, nb>::convert(core);750 }751 752 int integer() {753 return integer<int>();754 }755 756 /// \brief Returns a random bool757 ///758 /// It returns a random bool. The generator holds a buffer for759 /// random bits. Every time when it become empty the generator makes760 /// a new random word and fill the buffer up.761 bool boolean() {762 return bool_producer.convert(core);763 }764 765 /// @}766 767 ///\name Non-uniform Distributions768 ///769 ///@{770 771 /// \brief Returns a random bool with given probability of true result.772 ///773 /// It returns a random bool with given probability of true result.774 bool boolean(double p) {775 return operator()() < p;776 }777 778 /// Standard normal (Gauss) distribution779 780 /// Standard normal (Gauss) distribution.781 /// \note The Cartesian form of the Box-Muller782 /// transformation is used to generate a random normal distribution.783 double gauss()784 {785 double V1,V2,S;786 do {787 V1=2*real<double>()-1;788 V2=2*real<double>()-1;789 S=V1*V1+V2*V2;790 } while(S>=1);791 return std::sqrt(-2*std::log(S)/S)*V1;792 }793 /// Normal (Gauss) distribution with given mean and standard deviation794 795 /// Normal (Gauss) distribution with given mean and standard deviation.796 /// \sa gauss()797 double gauss(double mean,double std_dev)798 {799 return gauss()*std_dev+mean;800 }801 802 /// Lognormal distribution803 804 /// Lognormal distribution. The parameters are the mean and the standard805 /// deviation of <tt>exp(X)</tt>.806 ///807 double lognormal(double n_mean,double n_std_dev)808 {809 return std::exp(gauss(n_mean,n_std_dev));810 }811 /// Lognormal distribution812 813 /// Lognormal distribution. The parameter is an <tt>std::pair</tt> of814 /// the mean and the standard deviation of <tt>exp(X)</tt>.815 ///816 double lognormal(const std::pair<double,double> ¶ms)817 {818 return std::exp(gauss(params.first,params.second));819 }820 /// Compute the lognormal parameters from mean and standard deviation821 822 /// This function computes the lognormal parameters from mean and823 /// standard deviation. The return value can direcly be passed to824 /// lognormal().825 std::pair<double,double> lognormalParamsFromMD(double mean,826 double std_dev)827 {828 double fr=std_dev/mean;829 fr*=fr;830 double lg=std::log(1+fr);831 return std::pair<double,double>(std::log(mean)-lg/2.0,std::sqrt(lg));832 }833 /// Lognormal distribution with given mean and standard deviation834 835 /// Lognormal distribution with given mean and standard deviation.836 ///837 double lognormalMD(double mean,double std_dev)838 {839 return lognormal(lognormalParamsFromMD(mean,std_dev));840 }841 842 /// Exponential distribution with given mean843 844 /// This function generates an exponential distribution random number845 /// with mean <tt>1/lambda</tt>.846 ///847 double exponential(double lambda=1.0)848 {849 return -std::log(1.0-real<double>())/lambda;850 }851 852 /// Gamma distribution with given integer shape853 854 /// This function generates a gamma distribution random number.855 ///856 ///\param k shape parameter (<tt>k>0</tt> integer)857 double gamma(int k)858 {859 double s = 0;860 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());861 return s;862 }863 864 /// Gamma distribution with given shape and scale parameter865 866 /// This function generates a gamma distribution random number.867 ///868 ///\param k shape parameter (<tt>k>0</tt>)869 ///\param theta scale parameter870 ///871 double gamma(double k,double theta=1.0)872 {873 double xi,nu;874 const double delta = k-std::floor(k);875 const double v0=E/(E-delta);876 do {877 double V0=1.0-real<double>();878 double V1=1.0-real<double>();879 double V2=1.0-real<double>();880 if(V2<=v0)881 {882 xi=std::pow(V1,1.0/delta);883 nu=V0*std::pow(xi,delta-1.0);884 }885 else886 {887 xi=1.0-std::log(V1);888 nu=V0*std::exp(-xi);889 }890 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));891 return theta*(xi+gamma(int(std::floor(k))));892 }893 894 /// Weibull distribution895 896 /// This function generates a Weibull distribution random number.897 ///898 ///\param k shape parameter (<tt>k>0</tt>)899 ///\param lambda scale parameter (<tt>lambda>0</tt>)900 ///901 double weibull(double k,double lambda)902 {903 return lambda*pow(-std::log(1.0-real<double>()),1.0/k);904 }905 906 /// Pareto distribution907 908 /// This function generates a Pareto distribution random number.909 ///910 ///\param k shape parameter (<tt>k>0</tt>)911 ///\param x_min location parameter (<tt>x_min>0</tt>)912 ///913 double pareto(double k,double x_min)914 {915 return exponential(gamma(k,1.0/x_min))+x_min;916 }917 918 /// Poisson distribution919 920 /// This function generates a Poisson distribution random number with921 /// parameter \c lambda.922 ///923 /// The probability mass function of this distribusion is924 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]925 /// \note The algorithm is taken from the book of Donald E. Knuth titled926 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the927 /// return value.928 929 int poisson(double lambda)930 {931 const double l = std::exp(-lambda);932 int k=0;933 double p = 1.0;934 do {935 k++;936 p*=real<double>();937 } while (p>=l);938 return k-1;939 }940 941 ///@}942 943 ///\name Two Dimensional Distributions944 ///945 ///@{946 947 /// Uniform distribution on the full unit circle948 949 /// Uniform distribution on the full unit circle.950 ///951 dim2::Point<double> disc()952 {953 double V1,V2;954 do {955 V1=2*real<double>()-1;956 V2=2*real<double>()-1;957 958 } while(V1*V1+V2*V2>=1);959 return dim2::Point<double>(V1,V2);960 }961 /// A kind of two dimensional normal (Gauss) distribution962 963 /// This function provides a turning symmetric two-dimensional distribution.964 /// Both coordinates are of standard normal distribution, but they are not965 /// independent.966 ///967 /// \note The coordinates are the two random variables provided by968 /// the Box-Muller method.969 dim2::Point<double> gauss2()970 {971 double V1,V2,S;972 do {973 V1=2*real<double>()-1;974 V2=2*real<double>()-1;975 S=V1*V1+V2*V2;976 } while(S>=1);977 double W=std::sqrt(-2*std::log(S)/S);978 return dim2::Point<double>(W*V1,W*V2);979 }980 /// A kind of two dimensional exponential distribution981 982 /// This function provides a turning symmetric two-dimensional distribution.983 /// The x-coordinate is of conditionally exponential distribution984 /// with the condition that x is positive and y=0. If x is negative and985 /// y=0 then, -x is of exponential distribution. The same is true for the986 /// y-coordinate.987 dim2::Point<double> exponential2()988 {989 double V1,V2,S;990 do {991 V1=2*real<double>()-1;992 V2=2*real<double>()-1;993 S=V1*V1+V2*V2;994 } while(S>=1);995 double W=-std::log(S)/S;996 return dim2::Point<double>(W*V1,W*V2);997 }998 999 ///@}1000 };1001 1002 1003 extern Random rnd;1004 1005 }1006 1007 #endif -
test/random_test.cc
r463 r1379 22 22 int seed_array[] = {1, 2}; 23 23 24 int rnd_seq32[] = { 25 2732, 43567, 42613, 52416, 45891, 21243, 30403, 32103, 26 62501, 33003, 12172, 5192, 32511, 50057, 43723, 7813, 27 23720, 35343, 6637, 30280, 44566, 31019, 18898, 33867, 28 5994, 1688, 11513, 59011, 48056, 25544, 39168, 25365, 29 17530, 8366, 27063, 49861, 55169, 63848, 11863, 49608 30 }; 31 int rnd_seq64[] = { 32 56382, 63883, 59577, 64750, 9644, 59886, 57647, 18152, 33 28520, 64078, 17818, 49294, 26424, 26697, 53684, 19209, 34 35404, 12121, 12837, 11827, 32156, 58333, 62553, 7907, 35 64427, 39399, 21971, 48789, 46981, 15716, 53335, 65256, 36 12999, 15308, 10906, 42162, 47587, 43006, 53921, 18716 37 }; 38 39 void seq_test() { 40 for(int i=0;i<5;i++) { 41 lemon::Random32 r32(i); 42 lemon::Random64 r64(i); 43 for(int j=0;j<8;j++) { 44 check(r32[65536]==rnd_seq32[i*8+j], "Wrong random sequence"); 45 check(r64[65536]==rnd_seq64[i*8+j], "Wrong random sequence"); 46 } 47 } 48 } 49 50 24 51 int main() 25 52 { … … 37 64 (sizeof(seed_array) / sizeof(seed_array[0]))); 38 65 66 seq_test(); 39 67 return 0; 40 68 }
Note: See TracChangeset
for help on using the changeset viewer.