... |
... |
@@ -22,99 +22,100 @@
|
22 |
22 |
*
|
23 |
23 |
* See the appropriate copyright notice below.
|
24 |
24 |
*
|
25 |
25 |
* Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
|
26 |
26 |
* All rights reserved.
|
27 |
27 |
*
|
28 |
28 |
* Redistribution and use in source and binary forms, with or without
|
29 |
29 |
* modification, are permitted provided that the following conditions
|
30 |
30 |
* are met:
|
31 |
31 |
*
|
32 |
32 |
* 1. Redistributions of source code must retain the above copyright
|
33 |
33 |
* notice, this list of conditions and the following disclaimer.
|
34 |
34 |
*
|
35 |
35 |
* 2. Redistributions in binary form must reproduce the above copyright
|
36 |
36 |
* notice, this list of conditions and the following disclaimer in the
|
37 |
37 |
* documentation and/or other materials provided with the distribution.
|
38 |
38 |
*
|
39 |
39 |
* 3. The names of its contributors may not be used to endorse or promote
|
40 |
40 |
* products derived from this software without specific prior written
|
41 |
41 |
* permission.
|
42 |
42 |
*
|
43 |
43 |
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
44 |
44 |
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
45 |
45 |
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
|
46 |
46 |
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
|
47 |
47 |
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
|
48 |
48 |
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
|
49 |
49 |
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
|
50 |
50 |
* SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
51 |
51 |
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
|
52 |
52 |
* STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
|
53 |
53 |
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
|
54 |
54 |
* OF THE POSSIBILITY OF SUCH DAMAGE.
|
55 |
55 |
*
|
56 |
56 |
*
|
57 |
57 |
* Any feedback is very welcome.
|
58 |
58 |
* http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
|
59 |
59 |
* email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
|
60 |
60 |
*/
|
61 |
61 |
|
62 |
62 |
#ifndef LEMON_RANDOM_H
|
63 |
63 |
#define LEMON_RANDOM_H
|
64 |
64 |
|
65 |
65 |
#include <algorithm>
|
66 |
66 |
#include <iterator>
|
67 |
67 |
#include <vector>
|
68 |
68 |
|
69 |
69 |
#include <ctime>
|
70 |
|
#include <cmath>
|
71 |
70 |
|
|
71 |
#include <lemon/math.h>
|
72 |
72 |
#include <lemon/dim2.h>
|
|
73 |
|
73 |
74 |
///\ingroup misc
|
74 |
75 |
///\file
|
75 |
76 |
///\brief Mersenne Twister random number generator
|
76 |
77 |
|
77 |
78 |
namespace lemon {
|
78 |
79 |
|
79 |
80 |
namespace _random_bits {
|
80 |
81 |
|
81 |
82 |
template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
|
82 |
83 |
struct RandomTraits {};
|
83 |
84 |
|
84 |
85 |
template <typename _Word>
|
85 |
86 |
struct RandomTraits<_Word, 32> {
|
86 |
87 |
|
87 |
88 |
typedef _Word Word;
|
88 |
89 |
static const int bits = 32;
|
89 |
90 |
|
90 |
91 |
static const int length = 624;
|
91 |
92 |
static const int shift = 397;
|
92 |
93 |
|
93 |
94 |
static const Word mul = 0x6c078965u;
|
94 |
95 |
static const Word arrayInit = 0x012BD6AAu;
|
95 |
96 |
static const Word arrayMul1 = 0x0019660Du;
|
96 |
97 |
static const Word arrayMul2 = 0x5D588B65u;
|
97 |
98 |
|
98 |
99 |
static const Word mask = 0x9908B0DFu;
|
99 |
100 |
static const Word loMask = (1u << 31) - 1;
|
100 |
101 |
static const Word hiMask = ~loMask;
|
101 |
102 |
|
102 |
103 |
|
103 |
104 |
static Word tempering(Word rnd) {
|
104 |
105 |
rnd ^= (rnd >> 11);
|
105 |
106 |
rnd ^= (rnd << 7) & 0x9D2C5680u;
|
106 |
107 |
rnd ^= (rnd << 15) & 0xEFC60000u;
|
107 |
108 |
rnd ^= (rnd >> 18);
|
108 |
109 |
return rnd;
|
109 |
110 |
}
|
110 |
111 |
|
111 |
112 |
};
|
112 |
113 |
|
113 |
114 |
template <typename _Word>
|
114 |
115 |
struct RandomTraits<_Word, 64> {
|
115 |
116 |
|
116 |
117 |
typedef _Word Word;
|
117 |
118 |
static const int bits = 64;
|
118 |
119 |
|
119 |
120 |
static const int length = 312;
|
120 |
121 |
static const int shift = 156;
|
... |
... |
@@ -714,97 +715,97 @@
|
714 |
715 |
V1=2*real<double>()-1;
|
715 |
716 |
V2=2*real<double>()-1;
|
716 |
717 |
S=V1*V1+V2*V2;
|
717 |
718 |
} while(S>=1);
|
718 |
719 |
return std::sqrt(-2*std::log(S)/S)*V1;
|
719 |
720 |
}
|
720 |
721 |
/// Gauss distribution with given mean and standard deviation
|
721 |
722 |
|
722 |
723 |
/// Gauss distribution with given mean and standard deviation.
|
723 |
724 |
/// \sa gauss()
|
724 |
725 |
double gauss(double mean,double std_dev)
|
725 |
726 |
{
|
726 |
727 |
return gauss()*std_dev+mean;
|
727 |
728 |
}
|
728 |
729 |
|
729 |
730 |
/// Exponential distribution with given mean
|
730 |
731 |
|
731 |
732 |
/// This function generates an exponential distribution random number
|
732 |
733 |
/// with mean <tt>1/lambda</tt>.
|
733 |
734 |
///
|
734 |
735 |
double exponential(double lambda=1.0)
|
735 |
736 |
{
|
736 |
737 |
return -std::log(1.0-real<double>())/lambda;
|
737 |
738 |
}
|
738 |
739 |
|
739 |
740 |
/// Gamma distribution with given integer shape
|
740 |
741 |
|
741 |
742 |
/// This function generates a gamma distribution random number.
|
742 |
743 |
///
|
743 |
744 |
///\param k shape parameter (<tt>k>0</tt> integer)
|
744 |
745 |
double gamma(int k)
|
745 |
746 |
{
|
746 |
747 |
double s = 0;
|
747 |
748 |
for(int i=0;i<k;i++) s-=std::log(1.0-real<double>());
|
748 |
749 |
return s;
|
749 |
750 |
}
|
750 |
751 |
|
751 |
752 |
/// Gamma distribution with given shape and scale parameter
|
752 |
753 |
|
753 |
754 |
/// This function generates a gamma distribution random number.
|
754 |
755 |
///
|
755 |
756 |
///\param k shape parameter (<tt>k>0</tt>)
|
756 |
757 |
///\param theta scale parameter
|
757 |
758 |
///
|
758 |
759 |
double gamma(double k,double theta=1.0)
|
759 |
760 |
{
|
760 |
761 |
double xi,nu;
|
761 |
762 |
const double delta = k-std::floor(k);
|
762 |
|
const double v0=M_E/(M_E-delta);
|
|
763 |
const double v0=E/(E-delta);
|
763 |
764 |
do {
|
764 |
765 |
double V0=1.0-real<double>();
|
765 |
766 |
double V1=1.0-real<double>();
|
766 |
767 |
double V2=1.0-real<double>();
|
767 |
768 |
if(V2<=v0)
|
768 |
769 |
{
|
769 |
770 |
xi=std::pow(V1,1.0/delta);
|
770 |
771 |
nu=V0*std::pow(xi,delta-1.0);
|
771 |
772 |
}
|
772 |
773 |
else
|
773 |
774 |
{
|
774 |
775 |
xi=1.0-std::log(V1);
|
775 |
776 |
nu=V0*std::exp(-xi);
|
776 |
777 |
}
|
777 |
778 |
} while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
|
778 |
779 |
return theta*(xi-gamma(int(std::floor(k))));
|
779 |
780 |
}
|
780 |
781 |
|
781 |
782 |
/// Weibull distribution
|
782 |
783 |
|
783 |
784 |
/// This function generates a Weibull distribution random number.
|
784 |
785 |
///
|
785 |
786 |
///\param k shape parameter (<tt>k>0</tt>)
|
786 |
787 |
///\param lambda scale parameter (<tt>lambda>0</tt>)
|
787 |
788 |
///
|
788 |
789 |
double weibull(double k,double lambda)
|
789 |
790 |
{
|
790 |
791 |
return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
|
791 |
792 |
}
|
792 |
793 |
|
793 |
794 |
/// Pareto distribution
|
794 |
795 |
|
795 |
796 |
/// This function generates a Pareto distribution random number.
|
796 |
797 |
///
|
797 |
798 |
///\param k shape parameter (<tt>k>0</tt>)
|
798 |
799 |
///\param x_min location parameter (<tt>x_min>0</tt>)
|
799 |
800 |
///
|
800 |
801 |
double pareto(double k,double x_min)
|
801 |
802 |
{
|
802 |
803 |
return exponential(gamma(k,1.0/x_min));
|
803 |
804 |
}
|
804 |
805 |
|
805 |
806 |
///@}
|
806 |
807 |
|
807 |
808 |
///\name Two dimensional distributions
|
808 |
809 |
///
|
809 |
810 |
|
810 |
811 |
///@{
|