Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2012 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module mod_random
21 : use const_def, only: dp
22 :
23 : implicit none
24 :
25 : contains
26 :
27 0 : subroutine get_seed ( seed )
28 :
29 : ! *****************************************************************************80
30 : !
31 : !! GET_SEED returns a seed for the random number generator.
32 : !
33 : ! Discussion:
34 : !
35 : ! The seed depends on the current time, and ought to be (slightly)
36 : ! different every millisecond. Once the seed is obtained, a random
37 : ! number generator should be called a few times to further process
38 : ! the seed.
39 : !
40 : ! Licensing:
41 : !
42 : ! This code is distributed under the GNU LGPL license.
43 : !
44 : ! Modified:
45 : !
46 : ! 02 August 2004
47 : !
48 : ! Author:
49 : !
50 : ! John Burkardt
51 : !
52 : ! Parameters:
53 : !
54 : ! Output, integer ( kind = 4 ) SEED, a pseudorandom seed value.
55 : !
56 : implicit none
57 :
58 : integer ( kind = 4 ) :: seed
59 0 : real ( kind = 8 ) :: temp
60 : character ( len = 10 ) :: time
61 : character ( len = 8 ) :: today
62 : integer ( kind = 4 ) :: values(8)
63 : character ( len = 5 ) :: zone
64 :
65 0 : call date_and_time ( today, time, zone, values )
66 :
67 0 : temp = 0.0D+00
68 :
69 0 : temp = temp + real ( values(2) - 1, kind = 8 ) / 11.0D+00
70 0 : temp = temp + real ( values(3) - 1, kind = 8 ) / 30.0D+00
71 0 : temp = temp + real ( values(5), kind = 8 ) / 23.0D+00
72 0 : temp = temp + real ( values(6), kind = 8 ) / 59.0D+00
73 0 : temp = temp + real ( values(7), kind = 8 ) / 59.0D+00
74 0 : temp = temp + real ( values(8), kind = 8 ) / 999.0D+00
75 0 : temp = temp / 6.0D+00
76 :
77 0 : do while ( temp <= 0.0D+00 )
78 0 : temp = temp + 1.0D+00
79 : end do
80 :
81 0 : do while ( 1.0D+00 < temp )
82 0 : temp = temp - 1.0D+00
83 : end do
84 :
85 0 : seed = int ( real ( huge ( 1 ), kind = 8 ) * temp )
86 : !
87 : ! Never use a seed of 0 or maximum integer ( kind = 4 ).
88 : !
89 0 : if ( seed == 0 ) then
90 0 : seed = 1
91 : end if
92 :
93 0 : if ( seed == huge ( 1 ) ) then
94 0 : seed = seed - 1
95 : end if
96 :
97 0 : return
98 : end subroutine get_seed
99 :
100 :
101 0 : function i4_uniform ( a, b, seed )
102 :
103 : ! *****************************************************************************80
104 : !
105 : !! I4_UNIFORM returns a scaled pseudorandom I4.
106 : !
107 : ! Discussion:
108 : !
109 : ! An I4 is an integer ( kind = 4 ) value.
110 : !
111 : ! The pseudorandom number will be scaled to be uniformly distributed
112 : ! between A and B.
113 : !
114 : ! Licensing:
115 : !
116 : ! This code is distributed under the GNU LGPL license.
117 : !
118 : ! Modified:
119 : !
120 : ! 31 May 2007
121 : !
122 : ! Author:
123 : !
124 : ! John Burkardt
125 : !
126 : ! Reference:
127 : !
128 : ! Paul Bratley, Bennett Fox, Linus Schrage,
129 : ! A Guide to Simulation,
130 : ! Second Edition,
131 : ! Springer, 1987,
132 : ! ISBN: 0387964673,
133 : ! LC: QA76.9.C65.B73.
134 : !
135 : ! Bennett Fox,
136 : ! Algorithm 647:
137 : ! Implementation and Relative Efficiency of Quasirandom
138 : ! Sequence Generators,
139 : ! ACM Transactions on Mathematical Software,
140 : ! Volume 12, Number 4, December 1986, pages 362-376.
141 : !
142 : ! Pierre L'Ecuyer,
143 : ! Random Number Generation,
144 : ! in Handbook of Simulation,
145 : ! edited by Jerry Banks,
146 : ! Wiley, 1998,
147 : ! ISBN: 0471134031,
148 : ! LC: T57.62.H37.
149 : !
150 : ! Peter Lewis, Allen Goodman, James Miller,
151 : ! A Pseudo-Random Number Generator for the System/360,
152 : ! IBM Systems Journal,
153 : ! Volume 8, Number 2, 1969, pages 136-143.
154 : !
155 : ! Parameters:
156 : !
157 : ! Input, integer ( kind = 4 ) A, B, the limits of the interval.
158 : !
159 : ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which
160 : ! should NOT be 0. On output, SEED has been updated.
161 : !
162 : ! Output, integer ( kind = 4 ) I4_UNIFORM, a number between A and B.
163 : !
164 : implicit none
165 :
166 : integer ( kind = 4 ) :: a
167 : integer ( kind = 4 ) :: b
168 : integer ( kind = 4 ), parameter :: i4_huge = 2147483647
169 : integer ( kind = 4 ) :: i4_uniform
170 : integer ( kind = 4 ) :: k
171 0 : real ( kind = 4 ) :: r
172 : integer ( kind = 4 ) :: seed
173 : integer ( kind = 4 ) :: value
174 :
175 0 : if ( seed == 0 ) then
176 0 : write ( *, '(a)' ) ' '
177 0 : write ( *, '(a)' ) 'I4_UNIFORM - Fatal error!'
178 0 : write ( *, '(a)' ) ' Input value of SEED = 0.'
179 0 : stop
180 : end if
181 :
182 0 : k = seed / 127773
183 :
184 0 : seed = 16807 * ( seed - k * 127773 ) - k * 2836
185 :
186 0 : if ( seed < 0 ) then
187 0 : seed = seed + i4_huge
188 : end if
189 :
190 0 : r = real ( seed, kind = 4 ) * 4.656612875E-10
191 : !
192 : ! Scale R to lie between A-0.5 and B+0.5.
193 : !
194 : r = ( 1.0E+00 - r ) * ( real ( min ( a, b ), kind = 4 ) - 0.5E+00 ) &
195 0 : + r * ( real ( max ( a, b ), kind = 4 ) + 0.5E+00 )
196 : !
197 : ! Use rounding to convert R to an integer between A and B.
198 : !
199 0 : value = nint ( r, kind = 4 )
200 :
201 0 : value = max ( value, min ( a, b ) )
202 0 : value = min ( value, max ( a, b ) )
203 :
204 0 : i4_uniform = value
205 :
206 : return
207 : end function i4_uniform
208 :
209 :
210 0 : subroutine perm_uniform ( n, base, seed, p )
211 :
212 : ! *****************************************************************************80
213 : !
214 : !! PERM_UNIFORM selects a random permutation of N objects.
215 : !
216 : ! Licensing:
217 : !
218 : ! This code is distributed under the GNU LGPL license.
219 : !
220 : ! Modified:
221 : !
222 : ! 18 November 2008
223 : !
224 : ! Author:
225 : !
226 : ! John Burkardt
227 : !
228 : ! Reference:
229 : !
230 : ! Albert Nijenhuis, Herbert Wilf,
231 : ! Combinatorial Algorithms,
232 : ! Academic Press, 1978, second edition,
233 : ! ISBN 0-12-519260-6.
234 : !
235 : ! Parameters:
236 : !
237 : ! Input, integer ( kind = 4 ) N, the number of objects to be permuted.
238 : !
239 : ! Input, integer ( kind = 4 ) BASE, is 0 for a 0-based permutation and 1 for
240 : ! a 1-based permutation.
241 : !
242 : ! Input/output, integer ( kind = 4 ) SEED, a seed for the random
243 : ! number generator.
244 : !
245 : ! Output, integer ( kind = 4 ) P(N), the permutation. P(I) is the "new"
246 : ! location of the object originally at I.
247 : !
248 : implicit none
249 :
250 : integer ( kind = 4 ) :: n
251 :
252 : integer ( kind = 4 ) :: base
253 : integer ( kind = 4 ) :: i
254 : integer ( kind = 4 ) :: j
255 : integer ( kind = 4 ) :: k
256 : integer ( kind = 4 ) :: p(n)
257 : integer ( kind = 4 ) :: seed
258 :
259 0 : do i = 1, n
260 0 : p(i) = ( i - 1 ) + base
261 : end do
262 :
263 0 : do i = 1, n
264 0 : j = i4_uniform ( i, n, seed )
265 0 : k = p(i)
266 0 : p(i) = p(j)
267 0 : p(j) = k
268 : end do
269 :
270 0 : return
271 : end subroutine perm_uniform
272 :
273 :
274 :
275 25936 : function r8_uniform_01 ( seed )
276 :
277 : ! *****************************************************************************80
278 : !
279 : !! R8_UNIFORM_01 returns a unit pseudorandom R8.
280 : !
281 : ! Discussion:
282 : !
283 : ! An R8 is a real ( kind = 8 ) value.
284 : !
285 : ! For now, the input quantity SEED is an integer ( kind = 4 ) variable.
286 : !
287 : ! This routine implements the recursion
288 : !
289 : ! seed = 16807 * seed mod ( 2**31 - 1 )
290 : ! r8_uniform_01 = seed / ( 2**31 - 1 )
291 : !
292 : ! The integer ( kind = 4 ) arithmetic never requires more than 32 bits,
293 : ! including a sign bit.
294 : !
295 : ! If the initial seed is 12345, then the first three computations are
296 : !
297 : ! Input Output R8_UNIFORM_01
298 : ! SEED SEED
299 : !
300 : ! 12345 207482415 0.096616
301 : ! 207482415 1790989824 0.833995
302 : ! 1790989824 2035175616 0.947702
303 : !
304 : ! Licensing:
305 : !
306 : ! This code is distributed under the GNU LGPL license.
307 : !
308 : ! Modified:
309 : !
310 : ! 05 July 2006
311 : !
312 : ! Author:
313 : !
314 : ! John Burkardt
315 : !
316 : ! Reference:
317 : !
318 : ! Paul Bratley, Bennett Fox, Linus Schrage,
319 : ! A Guide to Simulation,
320 : ! Springer Verlag, pages 201-202, 1983.
321 : !
322 : ! Pierre L'Ecuyer,
323 : ! Random Number Generation,
324 : ! in Handbook of Simulation,
325 : ! edited by Jerry Banks,
326 : ! Wiley Interscience, page 95, 1998.
327 : !
328 : ! Bennett Fox,
329 : ! Algorithm 647:
330 : ! Implementation and Relative Efficiency of Quasirandom
331 : ! Sequence Generators,
332 : ! ACM Transactions on Mathematical Software,
333 : ! Volume 12, Number 4, pages 362-376, 1986.
334 : !
335 : ! Peter Lewis, Allen Goodman, James Miller
336 : ! A Pseudo-Random Number Generator for the System/360,
337 : ! IBM Systems Journal,
338 : ! Volume 8, pages 136-143, 1969.
339 : !
340 : ! Parameters:
341 : !
342 : ! Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should
343 : ! NOT be 0. On output, SEED has been updated.
344 : !
345 : ! Output, real ( kind = 8 ) R8_UNIFORM_01, a new pseudorandom variate,
346 : ! strictly between 0 and 1.
347 : !
348 : implicit none
349 :
350 : integer ( kind = 4 ) :: k
351 : real ( kind = 8 ) :: r8_uniform_01
352 : integer ( kind = 4 ) :: seed
353 :
354 25936 : if ( seed == 0 ) then
355 0 : write ( *, '(a)' ) ' '
356 0 : write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!'
357 0 : write ( *, '(a)' ) ' Input value of SEED = 0.'
358 0 : stop
359 : end if
360 :
361 25936 : k = seed / 127773
362 :
363 25936 : seed = 16807 * ( seed - k * 127773 ) - k * 2836
364 :
365 25936 : if ( seed < 0 ) then
366 278 : seed = seed + 2147483647
367 : end if
368 : !
369 : ! Although SEED can be represented exactly as a 32 bit integer ( kind = 4 ),
370 : ! it generally cannot be represented exactly as a 32 bit real number!
371 : !
372 25936 : r8_uniform_01 = real ( seed, kind = 8 ) * 4.656612875D-10
373 :
374 25936 : end function r8_uniform_01
375 :
376 : end module mod_random
|
