blob: be81d16878487dbc4f3626d15cd7af5616759fd8 [file] [log] [blame]
avm9996345a8a462022-06-04 12:41:03 +02001c
2c A C-program for MT19937, with initialization improved 2002/1/26.
3c Coded by Takuji Nishimura and Makoto Matsumoto.
4c
5c Before using, initialize the state by using init_genrand(seed)
6c or init_by_array(init_key, key_length).
7c
8c Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
9c All rights reserved.
10c Copyright (C) 2005, Mutsuo Saito,
11c All rights reserved.
12c
13c Redistribution and use in source and binary forms, with or without
14c modification, are permitted provided that the following conditions
15c are met:
16c
17c 1. Redistributions of source code must retain the above copyright
18c notice, this list of conditions and the following disclaimer.
19c
20c 2. Redistributions in binary form must reproduce the above copyright
21c notice, this list of conditions and the following disclaimer in the
22c documentation and/or other materials provided with the distribution.
23c
24c 3. The names of its contributors may not be used to endorse or promote
25c products derived from this software without specific prior written
26c permission.
27c
28c THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
29c "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
30c LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
31c A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
32c CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
33c EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
34c PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
35c PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
36c LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
37c NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
38c SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
39c
40c
41c Any feedback is very welcome.
42c http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
43c email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
44c
45c-----------------------------------------------------------------------
46c FORTRAN77 translation by Tsuyoshi TADA. (2005/12/19)
47c
48c ---------- initialize routines ----------
49c subroutine init_genrand(seed): initialize with a seed
50c subroutine init_by_array(init_key,key_length): initialize by an array
51c
52c ---------- generate functions ----------
53c integer function genrand_int32(): signed 32-bit integer
54c integer function genrand_int31(): unsigned 31-bit integer
55c double precision function genrand_real1(): [0,1] with 32-bit resolution
56c double precision function genrand_real2(): [0,1) with 32-bit resolution
57c double precision function genrand_real3(): (0,1) with 32-bit resolution
58c double precision function genrand_res53(): (0,1) with 53-bit resolution
59c
60c This program uses the following non-standard intrinsics.
61c ishft(i,n): If n>0, shifts bits in i by n positions to left.
62c If n<0, shifts bits in i by n positions to right.
63c iand (i,j): Performs logical AND on corresponding bits of i and j.
64c ior (i,j): Performs inclusive OR on corresponding bits of i and j.
65c ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
66c
67c-----------------------------------------------------------------------
68c initialize mt(0:N-1) with a seed
69c-----------------------------------------------------------------------
70 subroutine init_genrand(s)
71 integer s
72 integer N
73 integer DONE
74 integer ALLBIT_MASK
75 parameter (N=624)
76 parameter (DONE=123456789)
77 integer mti,initialized
78 integer mt(0:N-1)
79 common /mt_state1/ mti,initialized
80 common /mt_state2/ mt
81 common /mt_mask1/ ALLBIT_MASK
82c
83 call mt_initln
84 mt(0)=iand(s,ALLBIT_MASK)
85 do 100 mti=1,N-1
86 mt(mti)=1812433253*
87 & ieor(mt(mti-1),ishft(mt(mti-1),-30))+mti
88 mt(mti)=iand(mt(mti),ALLBIT_MASK)
89 100 continue
90 initialized=DONE
91c
92 return
93 end
94c-----------------------------------------------------------------------
95c initialize by an array with array-length
96c init_key is the array for initializing keys
97c key_length is its length
98c-----------------------------------------------------------------------
99 subroutine init_by_array(init_key,key_length)
100 integer init_key(0:*)
101 integer key_length
102 integer N
103 integer ALLBIT_MASK
104 integer TOPBIT_MASK
105 parameter (N=624)
106 integer i,j,k
107 integer mt(0:N-1)
108 common /mt_state2/ mt
109 common /mt_mask1/ ALLBIT_MASK
110 common /mt_mask2/ TOPBIT_MASK
111c
112 call init_genrand(19650218)
113 i=1
114 j=0
115 do 100 k=max(N,key_length),1,-1
116 mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1664525)
117 & +init_key(j)+j
118 mt(i)=iand(mt(i),ALLBIT_MASK)
119 i=i+1
120 j=j+1
121 if(i.ge.N)then
122 mt(0)=mt(N-1)
123 i=1
124 endif
125 if(j.ge.key_length)then
126 j=0
127 endif
128 100 continue
129 do 200 k=N-1,1,-1
130 mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1566083941)-i
131 mt(i)=iand(mt(i),ALLBIT_MASK)
132 i=i+1
133 if(i.ge.N)then
134 mt(0)=mt(N-1)
135 i=1
136 endif
137 200 continue
138 mt(0)=TOPBIT_MASK
139c
140 return
141 end
142c-----------------------------------------------------------------------
143c generates a random number on [0,0xffffffff]-interval
144c-----------------------------------------------------------------------
145 function genrand_int32()
146 integer genrand_int32
147 integer N,M
148 integer DONE
149 integer UPPER_MASK,LOWER_MASK,MATRIX_A
150 integer T1_MASK,T2_MASK
151 parameter (N=624)
152 parameter (M=397)
153 parameter (DONE=123456789)
154 integer mti,initialized
155 integer mt(0:N-1)
156 integer y,kk
157 integer mag01(0:1)
158 common /mt_state1/ mti,initialized
159 common /mt_state2/ mt
160 common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK
161 common /mt_mag01/ mag01
162c
163 if(initialized.ne.DONE)then
164 call init_genrand(21641)
165 endif
166c
167 if(mti.ge.N)then
168 do 100 kk=0,N-M-1
169 y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK))
170 mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
171 100 continue
172 do 200 kk=N-M,N-1-1
173 y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK))
174 mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
175 200 continue
176 y=ior(iand(mt(N-1),UPPER_MASK),iand(mt(0),LOWER_MASK))
177 mt(kk)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
178 mti=0
179 endif
180c
181 y=mt(mti)
182 mti=mti+1
183c
184 y=ieor(y,ishft(y,-11))
185 y=ieor(y,iand(ishft(y,7),T1_MASK))
186 y=ieor(y,iand(ishft(y,15),T2_MASK))
187 y=ieor(y,ishft(y,-18))
188c
189 genrand_int32=y
190 return
191 end
192c-----------------------------------------------------------------------
193c generates a random number on [0,0x7fffffff]-interval
194c-----------------------------------------------------------------------
195 function genrand_int31()
196 integer genrand_int31
197 integer genrand_int32
198 genrand_int31=int(ishft(genrand_int32(),-1))
199 return
200 end
201c-----------------------------------------------------------------------
202c generates a random number on [0,1]-real-interval
203c-----------------------------------------------------------------------
204 function genrand_real1()
205 double precision genrand_real1,r
206 integer genrand_int32
207 r=dble(genrand_int32())
208 if(r.lt.0.d0)r=r+2.d0**32
209 genrand_real1=r/4294967295.d0
210 return
211 end
212c-----------------------------------------------------------------------
213c generates a random number on [0,1)-real-interval
214c-----------------------------------------------------------------------
215 function genrand_real2()
216 double precision genrand_real2,r
217 integer genrand_int32
218 r=dble(genrand_int32())
219 if(r.lt.0.d0)r=r+2.d0**32
220 genrand_real2=r/4294967296.d0
221 return
222 end
223c-----------------------------------------------------------------------
224c generates a random number on (0,1)-real-interval
225c-----------------------------------------------------------------------
226 function genrand_real3()
227 double precision genrand_real3,r
228 integer genrand_int32
229 r=dble(genrand_int32())
230 if(r.lt.0.d0)r=r+2.d0**32
231 genrand_real3=(r+0.5d0)/4294967296.d0
232 return
233 end
234c-----------------------------------------------------------------------
235c generates a random number on [0,1) with 53-bit resolution
236c-----------------------------------------------------------------------
237 function genrand_res53()
238 double precision genrand_res53
239 integer genrand_int32
240 double precision a,b
241 a=dble(ishft(genrand_int32(),-5))
242 b=dble(ishft(genrand_int32(),-6))
243 if(a.lt.0.d0)a=a+2.d0**32
244 if(b.lt.0.d0)b=b+2.d0**32
245 genrand_res53=(a*67108864.d0+b)/9007199254740992.d0
246 return
247 end
248c-----------------------------------------------------------------------
249c initialize large number (over 32-bit constant number)
250c-----------------------------------------------------------------------
251 subroutine mt_initln
252 integer ALLBIT_MASK
253 integer TOPBIT_MASK
254 integer UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK
255 integer mag01(0:1)
256 common /mt_mask1/ ALLBIT_MASK
257 common /mt_mask2/ TOPBIT_MASK
258 common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK
259 common /mt_mag01/ mag01
260CC TOPBIT_MASK = Z'80000000'
261CC ALLBIT_MASK = Z'ffffffff'
262CC UPPER_MASK = Z'80000000'
263CC LOWER_MASK = Z'7fffffff'
264CC MATRIX_A = Z'9908b0df'
265CC T1_MASK = Z'9d2c5680'
266CC T2_MASK = Z'efc60000'
267 TOPBIT_MASK=1073741824
268 TOPBIT_MASK=ishft(TOPBIT_MASK,1)
269 ALLBIT_MASK=2147483647
270 ALLBIT_MASK=ior(ALLBIT_MASK,TOPBIT_MASK)
271 UPPER_MASK=TOPBIT_MASK
272 LOWER_MASK=2147483647
273 MATRIX_A=419999967
274 MATRIX_A=ior(MATRIX_A,TOPBIT_MASK)
275 T1_MASK=489444992
276 T1_MASK=ior(T1_MASK,TOPBIT_MASK)
277 T2_MASK=1875247104
278 T2_MASK=ior(T2_MASK,TOPBIT_MASK)
279 mag01(0)=0
280 mag01(1)=MATRIX_A
281 return
282 end