forked from qutip/qutip
/
mt19937.f90
276 lines (212 loc) · 8.07 KB
/
mt19937.f90
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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
! A Fortran-program for MT19937: Real number version
! Code converted using TO_F90 by Alan Miller
! Date: 1999-11-26 Time: 17:09:23
! Latest revision - 5 February 2002
! A new seed initialization routine has been added based upon the new
! C version dated 26 January 2002.
! This version assumes that integer overflows do NOT cause crashes.
! This version is compatible with Lahey's ELF90 compiler,
! and should be compatible with most full Fortran 90 or 95 compilers.
! Notice the strange way in which umask is specified for ELF90.
! genrand() generates one pseudorandom real number (double) which is
! uniformly distributed on [0,1]-interval, for each call.
! sgenrand(seed) set initial values to the working area of 624 words.
! Before genrand(), sgenrand(seed) must be called once. (seed is any 32-bit
! integer except for 0).
! Integer generator is obtained by modifying two lines.
! Coded by Takuji Nishimura, considering the suggestions by
! Topher Cooper and Marc Rieffel in July-Aug. 1997.
! This library is free software; you can redistribute it and/or modify it
! under the terms of the GNU Library General Public License as published by
! the Free Software Foundation; either version 2 of the License, or (at your
! option) any later version. This library is distributed in the hope that
! it will be useful, but WITHOUT ANY WARRANTY; without even the implied
! warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
! See the GNU Library General Public License for more details.
! You should have received a copy of the GNU Library General Public License
! along with this library; if not, write to the Free Foundation, Inc.,
! 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
! Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
! When you use this, send an email to: matumoto@math.keio.ac.jp
! with an appropriate reference to your work.
!***********************************************************************
! Fortran translation by Hiroshi Takano. Jan. 13, 1999.
! genrand() -> double precision function grnd()
! sgenrand(seed) -> subroutine sgrnd(seed)
! integer seed
! This program uses the following standard intrinsics.
! ishft(i,n): If n > 0, shifts bits in i by n positions to left.
! If n < 0, shifts bits in i by n positions to right.
! iand (i,j): Performs logical AND on corresponding bits of i and j.
! ior (i,j): Performs inclusive OR on corresponding bits of i and j.
! ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
!***********************************************************************
MODULE mt19937
use qutraj_precision
IMPLICIT NONE
!INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60)
! Period parameters
INTEGER, PARAMETER :: n = 624, n1 = n+1, m = 397, mata = -1727483681
! constant vector a
INTEGER, PARAMETER :: umask = -2147483647 - 1
! most significant w-r bits
INTEGER, PARAMETER :: lmask = 2147483647
! least significant r bits
! Tempering parameters
INTEGER, PARAMETER :: tmaskb= -1658038656, tmaskc= -272236544
! the array for the state vector
INTEGER, SAVE :: mt(0:n-1), mti = n1
! mti==N+1 means mt[N] is not initialized
PRIVATE
PUBLIC :: dp, sgrnd, grnd, init_genrand
CONTAINS
SUBROUTINE sgrnd(seed)
! This is the original version of the seeding routine.
! It was replaced in the Japanese version in C on 26 January 2002
! It is recommended that routine init_genrand is used instead.
INTEGER, INTENT(IN) :: seed
! setting initial seeds to mt[N] using the generator Line 25 of Table 1 in
! [KNUTH 1981, The Art of Computer Programming Vol. 2 (2nd Ed.), pp102]
mt(0)= IAND(seed, -1)
DO mti=1,n-1
mt(mti) = IAND(69069 * mt(mti-1), -1)
END DO
RETURN
END SUBROUTINE sgrnd
!***********************************************************************
SUBROUTINE init_genrand(seed)
! This initialization is based upon the multiplier given on p.106 of the
! 3rd edition of Knuth, The Art of Computer Programming Vol. 2.
! This version assumes that integer overflow does NOT cause a crash.
INTEGER, INTENT(IN) :: seed
INTEGER :: latest
mt(0) = seed
latest = seed
DO mti = 1, n-1
latest = IEOR( latest, ISHFT( latest, -30 ) )
latest = latest * 1812433253 + mti
mt(mti) = latest
END DO
RETURN
END SUBROUTINE init_genrand
!***********************************************************************
FUNCTION grnd() RESULT(fn_val)
use qutraj_precision
REAL (dp) :: fn_val
INTEGER, SAVE :: mag01(0:1) = (/ 0, mata /)
! mag01(x) = x * MATA for x=0,1
INTEGER :: kk, y
! These statement functions have been replaced with separate functions
! tshftu(y) = ISHFT(y,-11)
! tshfts(y) = ISHFT(y,7)
! tshftt(y) = ISHFT(y,15)
! tshftl(y) = ISHFT(y,-18)
IF(mti >= n) THEN
! generate N words at one time
IF(mti == n+1) THEN
! if sgrnd() has not been called,
CALL sgrnd(4357)
! a default initial seed is used
END IF
DO kk = 0, n-m-1
y = IOR(IAND(mt(kk),umask), IAND(mt(kk+1),lmask))
mt(kk) = IEOR(IEOR(mt(kk+m), ISHFT(y,-1)),mag01(IAND(y,1)))
END DO
DO kk = n-m, n-2
y = IOR(IAND(mt(kk),umask), IAND(mt(kk+1),lmask))
mt(kk) = IEOR(IEOR(mt(kk+(m-n)), ISHFT(y,-1)),mag01(IAND(y,1)))
END DO
y = IOR(IAND(mt(n-1),umask), IAND(mt(0),lmask))
mt(n-1) = IEOR(IEOR(mt(m-1), ISHFT(y,-1)),mag01(IAND(y,1)))
mti = 0
END IF
y = mt(mti)
mti = mti + 1
y = IEOR(y, tshftu(y))
y = IEOR(y, IAND(tshfts(y),tmaskb))
y = IEOR(y, IAND(tshftt(y),tmaskc))
y = IEOR(y, tshftl(y))
IF(y < 0) THEN
fn_val = (DBLE(y) + 2.0D0**32) / (2.0D0**32 - 1.0D0)
ELSE
fn_val = DBLE(y) / (2.0D0**32 - 1.0D0)
END IF
RETURN
END FUNCTION grnd
FUNCTION tshftu(y) RESULT(fn_val)
INTEGER, INTENT(IN) :: y
INTEGER :: fn_val
fn_val = ISHFT(y,-11)
RETURN
END FUNCTION tshftu
FUNCTION tshfts(y) RESULT(fn_val)
INTEGER, INTENT(IN) :: y
INTEGER :: fn_val
fn_val = ISHFT(y,7)
RETURN
END FUNCTION tshfts
FUNCTION tshftt(y) RESULT(fn_val)
INTEGER, INTENT(IN) :: y
INTEGER :: fn_val
fn_val = ISHFT(y,15)
RETURN
END FUNCTION tshftt
FUNCTION tshftl(y) RESULT(fn_val)
INTEGER, INTENT(IN) :: y
INTEGER :: fn_val
fn_val = ISHFT(y,-18)
RETURN
END FUNCTION tshftl
END MODULE mt19937
! this main() outputs the first 1000 generated numbers
PROGRAM main
USE mt19937
IMPLICIT NONE
INTEGER, PARAMETER :: no = 1000
INTEGER :: count, j, k, seed
REAL (dp) :: temp, r(0:7), big, small, average, sumsq, stdev
! call sgrnd(4357)
! any nonzero integer can be used as a seed
WRITE(*, '(a)', ADVANCE='NO') ' Enter integer random number seed: '
READ(*, *) seed
CALL init_genrand(seed)
WRITE(*, *) 'Seed = ', seed
big = 0.5_dp
small = 0.5_dp
count = 0
average = 0.0_dp
sumsq = 0.0_dp
DO j=0,no-1
temp = grnd()
IF (temp > big) THEN
big = temp
ELSE IF (temp < small) THEN
small = temp
END IF
CALL update(temp, count, average, sumsq)
r(MOD(j,8)) = temp
IF(MOD(j,8) == 7) THEN
WRITE(*, '(8(f8.6, '' ''))') r(0:7)
ELSE IF(j == no-1) THEN
WRITE(*, '(8(f8.6, '' ''))') (r(k),k=0,MOD(no-1,8))
END IF
END DO
stdev = SQRT( sumsq / (count - 1) )
WRITE(*, '(a, f10.6, a, f10.6)') ' Smallest = ', small, ' Largest = ', big
WRITE(*, '(a, f9.5, a, f9.5)') ' Average = ', average, ' Std. devn. = ', stdev
WRITE(*, *) ' Std. devn. should be about 1/sqrt(12) = 0.288675'
STOP
CONTAINS
SUBROUTINE update(x, n, avge, sumsq)
REAL (dp), INTENT(IN) :: x
INTEGER, INTENT(IN OUT) :: n
REAL (dp), INTENT(IN OUT) :: avge, sumsq
REAL (dp) :: dev
n = n + 1
dev = x - avge
avge = avge + dev / n
sumsq = sumsq + dev*(x - avge)
RETURN
END SUBROUTINE update
END PROGRAM main
!***********************************************************************