-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrandom_numbers.cpp
178 lines (142 loc) · 3.51 KB
/
random_numbers.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
/*Random number generators. ranf0 is to be used any time we need a uniform random
number, excluding as a test number for a Gaussian generator. ranf1,2,3 are to be
used with gaussian1,2,3, respectively.*/
double ranf0()
/* Uniform random number generator x(n+1)= a*x(n) mod c
with a = pow(7,5) and c = pow(2,31)-1.
Copyright (c) Tao Pang 1997. */
{
const int ia=16807,ic=2147483647,iq=127773,ir=2836;
int il,ih,it;
double rc;
extern int iseed0;
ih = iseed0/iq;
il = iseed0%iq;
it = ia*il-ir*ih;
if (it > 0)
{
iseed0 = it;
}
else
{
iseed0 = ic+it;
}
rc = ic;
return iseed0/rc;
}//end of ranf0()
double ranf2()
/* Uniform random number generator x(n+1)= a*x(n) mod c
with a = pow(7,5) and c = pow(2,31)-1.
Copyright (c) Tao Pang 1997. */
{
const int ia=16807,ic=2147483647,iq=127773,ir=2836;
int il,ih,it;
double rc;
extern int iseed2;
ih = iseed2/iq;
il = iseed2%iq;
it = ia*il-ir*ih;
if (it > 0)
{
iseed2 = it;
}
else
{
iseed2 = ic+it;
}
rc = ic;
return iseed2/rc;
}//end of ranf2()
double ranf6()
/* Uniform random number generator x(n+1)= a*x(n) mod c
* with a = pow(7,5) and c = pow(2,31)-1.
* Copyright (c) Tao Pang 1997. */
{
const int ia=16807,ic=2147483647,iq=127773,ir=2836;
int il,ih,it;
double rc;
extern int iseed6;
ih = iseed6/iq;
il = iseed6%iq;
it = ia*il-ir*ih;
if (it > 0)
{
iseed6 = it;
}
else
{
iseed6 = ic+it;
}
rc = ic;
return iseed6/rc;
}//end of ranf6()
/*Takes uniformly distributed random numbers and gives random number according to
Gaussian distribution by the rejection method.*/
//rotational brownian motion
double gaussian2()
{
double y, py, ptest;
y = gauss_range2*(ranf2() - 0.5);
ptest = gauss_prefact2*ranf2();//gauss_prefact is the maximum value of P_y(y), the probability density function
py = gauss_prefact2*exp(gauss_exp_const2*y*y);
if(ptest < py)
return(y);
else
return gaussian2();
}
///////////////////////////////////////////////////////////////////
double ranf5()
/* Uniform random number generator x(n+1)= a*x(n) mod c
* with a = pow(7,5) and c = pow(2,31)-1.
* Copyright (c) Tao Pang 1997. */
{
const int ia=16807,ic=2147483647,iq=127773,ir=2836;
int il,ih,it;
double rc;
extern int iseed5;
ih = iseed5/iq;
il = iseed5%iq;
it = ia*il-ir*ih;
if (it > 0)
{
iseed5 = it;
}
else
{
iseed5 = ic+it;
}
rc = ic;
return iseed5/rc;
}//end of ranf5()
//gaussian random number generator for heavy polymer mono rotation
double gaussian6()
{
double y, py, ptest;
y = gauss_range6*(ranf6() - 0.5);
ptest = gauss_prefact6*ranf6();//gauss_prefact is the maximum value of P_y(y), the probability density function
py = gauss_prefact6*exp(gauss_exp_const6*y*y);
if(ptest < py)
return(y);
else
return gaussian6();
}
/*****************/
double gaussian_std()
{
double y, py, ptest;
y = gauss_range_std*(ranf_std() - 0.5);
ptest = ranf_std();//gauss_prefact is the maximum value of P_y(y), the probability density function
py = exp(gauss_exp_const_std*y*y);
if(ptest < py)
return(y);
else
return gaussian_std();
}
double gaussian_inverf(unsigned long long step)
{
double r;
do{
r = RNUM0.get_double();
}while(r < 1.e-16);
return (SQRT_TWO* erf_inv( 2.*(r-0.5) ));
}