1 | MODULE mt19937ar |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE mt19937ar *** |
---|
4 | !! Mersen Twister random generator |
---|
5 | !!====================================================================== |
---|
6 | !! |
---|
7 | !! ** Purpose : |
---|
8 | !! ** Method : |
---|
9 | !! References : This a a Fortran90 translation of mt19937ar.c: |
---|
10 | !! A C-program for MT19937, with initialization improved 2002/1/26. |
---|
11 | !! Coded by Takuji Nishimura and Makoto Matsumoto. |
---|
12 | !! |
---|
13 | !! Before using, initialize the state by using init_genrand(seed) |
---|
14 | !! or init_by_array(init_key, key_length). |
---|
15 | !! |
---|
16 | !! Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, |
---|
17 | !! All rights reserved. |
---|
18 | !! |
---|
19 | !! Redistribution and use in source and binary forms, with or without |
---|
20 | !! modification, are permitted provided that the following conditions |
---|
21 | !! are met: |
---|
22 | !! |
---|
23 | !! 1. Redistributions of source code must retain the above copyright |
---|
24 | !! notice, this list of conditions and the following disclaimer. |
---|
25 | !! |
---|
26 | !! 2. Redistributions in binary form must reproduce the above copyright |
---|
27 | !! notice, this list of conditions and the following disclaimer in the |
---|
28 | !! documentation and/or other materials provided with the distribution. |
---|
29 | !! |
---|
30 | !! 3. The names of its contributors may not be used to endorse or promote |
---|
31 | !! products derived from this software without specific prior written |
---|
32 | !! permission. |
---|
33 | !! |
---|
34 | !! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
---|
35 | !! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
---|
36 | !! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
---|
37 | !! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR |
---|
38 | !! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
---|
39 | !! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
---|
40 | !! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
---|
41 | !! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
---|
42 | !! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
---|
43 | !! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
---|
44 | !! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
---|
45 | !! |
---|
46 | !! |
---|
47 | !! Any feedback is very welcome. |
---|
48 | !! http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html |
---|
49 | !! email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) |
---|
50 | !! |
---|
51 | !! Until 2001/4/6, MT had been distributed under GNU Public License, |
---|
52 | !! but after 2001/4/6, we decided to let MT be used for any purpose, |
---|
53 | !! including commercial use. 2002-versions mt19937ar.c, mt19937ar-cok.c |
---|
54 | !! are considered to be usable freely. |
---|
55 | !! |
---|
56 | !! History: |
---|
57 | !! ! 97 - 02 (Makoto Matsumoto/Takuji Nishimura) original c version |
---|
58 | !! ! 10-03 (A. Vidard) F90 translation and NEMOTAM adaptation |
---|
59 | !-----------------------------------------------------------------------! |
---|
60 | USE par_kind, ONLY : wp |
---|
61 | PRIVATE |
---|
62 | INTEGER, PARAMETER :: jpn = 624 ! period parameters |
---|
63 | INTEGER, PARAMETER :: jpm = 397 |
---|
64 | |
---|
65 | INTEGER*8, PARAMETER :: jpa = Z'9908b0df' ! constant vector A |
---|
66 | INTEGER*8, PARAMETER :: jpmsb = Z'80000000' ! most significant w-r bit mask |
---|
67 | INTEGER*8, PARAMETER :: jplsb = Z'7fffffff' ! least significant r bit mask |
---|
68 | ! |
---|
69 | INTEGER*8, PARAMETER :: jpall = Z'ffffffff' ! all bits mask |
---|
70 | |
---|
71 | INTEGER*8, PARAMETER :: jptmp1= Z'9d2c5680' ! Tampering paramters |
---|
72 | INTEGER*8, PARAMETER :: jptmp2= Z'efc60000' |
---|
73 | |
---|
74 | INTEGER :: mti |
---|
75 | INTEGER*8, DIMENSION(0:jpn-1) :: mt |
---|
76 | LOGICAL :: l_mtinit = .FALSE. |
---|
77 | |
---|
78 | PUBLIC & |
---|
79 | & l_mtinit, & |
---|
80 | & mtrand_int31, & |
---|
81 | & mtrand_int32, & |
---|
82 | & mtrand_real1, & |
---|
83 | & mtrand_real2, & |
---|
84 | & mtrand_real3, & |
---|
85 | & mtrand_real53,& |
---|
86 | & init_by_array,& |
---|
87 | & init_mtrand |
---|
88 | CONTAINS |
---|
89 | SUBROUTINE init_mtrand(ks) |
---|
90 | !----------------------------------------------------------------------- |
---|
91 | !! initialize mt(0:N-1) with a seed |
---|
92 | !----------------------------------------------------------------------- |
---|
93 | INTEGER :: ks |
---|
94 | ! |
---|
95 | mt(0) = IAND(INT(ks,kind=8),jpall) |
---|
96 | DO mti = 1, jpn-1 |
---|
97 | mt(mti) = 1812433253 * IEOR(mt(mti-1),ISHFT(mt(mti-1),-30)) + mti |
---|
98 | ! See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. |
---|
99 | ! In the previous versions, MSBs of the seed affect |
---|
100 | ! only MSBs of the array mt[]. |
---|
101 | ! modified by Makoto Matsumoto |
---|
102 | mt(mti) = IAND(mt(mti),jpall) |
---|
103 | ! for >32 bit machines |
---|
104 | END DO |
---|
105 | l_mtinit = .TRUE. |
---|
106 | ! |
---|
107 | END SUBROUTINE init_mtrand |
---|
108 | SUBROUTINE init_by_array(kinit_key,key_length) |
---|
109 | !----------------------------------------------------------------------- |
---|
110 | ! initialize by an array with array-length |
---|
111 | ! init_key is the array for initializing keys |
---|
112 | ! key_length is its length |
---|
113 | !----------------------------------------------------------------------- |
---|
114 | INTEGER :: kinit_key(0:*) |
---|
115 | INTEGER :: key_length |
---|
116 | INTEGER :: ji,jj,jk |
---|
117 | ! |
---|
118 | CALL init_mtrand(19650218) |
---|
119 | ji = 1 |
---|
120 | jj = 0 |
---|
121 | DO jk = MAX(jpn,key_length),1,-1 |
---|
122 | mt(ji) = IEOR(mt(ji),IEOR(mt(ji-1),ISHFT(mt(ji-1),-30)) * 1664525) & |
---|
123 | & + kinit_key(jj) + jj |
---|
124 | mt(ji) = IAND(mt(ji),jpall) |
---|
125 | ji = ji + 1 |
---|
126 | jj = jj + 1 |
---|
127 | IF (ji .GE. jpn) THEN |
---|
128 | mt(0) = mt(jpn-1) |
---|
129 | ji = 1 |
---|
130 | ENDIF |
---|
131 | IF (jj .GE. key_length) THEN |
---|
132 | jj = 0 |
---|
133 | ENDIF |
---|
134 | END DO |
---|
135 | DO jk = jpn-1,1,-1 |
---|
136 | mt(ji) = IEOR(mt(ji),IEOR(mt(ji-1),ISHFT(mt(ji-1),-30)) * 1566083941) - ji ! non linear |
---|
137 | mt(ji) = IAND(mt(ji),jpall) ! for WORDSIZE > 32 machines |
---|
138 | ji = ji + 1 |
---|
139 | IF (ji .GE. jpn) THEN |
---|
140 | mt(0) = mt(jpn-1) |
---|
141 | ji = 1 |
---|
142 | ENDIF |
---|
143 | END DO |
---|
144 | mt(0) = jpmsb ! MSB is 1; assuring non-zero initial array |
---|
145 | ! |
---|
146 | END SUBROUTINE init_by_array |
---|
147 | FUNCTION mtrand_int32() |
---|
148 | !----------------------------------------------------------------------- |
---|
149 | ! generates a random number on [0,0xffffffff]-interval |
---|
150 | !----------------------------------------------------------------------- |
---|
151 | INTEGER :: mtrand_int32 |
---|
152 | INTEGER :: jk |
---|
153 | INTEGER*8 :: iy |
---|
154 | INTEGER*8, DIMENSION(0:1) :: mag01 |
---|
155 | ! |
---|
156 | mag01(0) = 0 |
---|
157 | mag01(1) = jpa ! mag01[x] = x * vector A for x=0,1 |
---|
158 | ! |
---|
159 | IF(.NOT. l_mtinit)THEN ! if init_mtrand() has not been called |
---|
160 | CALL init_mtrand(5489) ! a default initial seed is used |
---|
161 | ENDIF |
---|
162 | ! |
---|
163 | IF (mti .GE. jpn) THEN |
---|
164 | DO jk = 0,jpn-jpm-1 |
---|
165 | iy = IOR(IAND(mt(jk),jpmsb),IAND(mt(jk+1),jplsb)) |
---|
166 | mt(jk) = IEOR(IEOR(mt(jk+jpm),ISHFT(iy,-1)),mag01(IAND(iy,INT(1,kind=8)))) |
---|
167 | END DO |
---|
168 | DO jk = jpn-jpm,jpn-1-1 |
---|
169 | iy = IOR(IAND(mt(jk),jpmsb),IAND(mt(jk+1),jplsb)) |
---|
170 | mt(jk) = IEOR(IEOR(mt(jk+(jpm-jpn)),ISHFT(iy,-1)),mag01(IAND(iy,INT(1,kind=8)))) |
---|
171 | END DO |
---|
172 | iy = IOR(IAND(mt(jpn-1),jpmsb),IAND(mt(0),jplsb)) |
---|
173 | mt(jk) = IEOR(IEOR(mt(jpm-1),ISHFT(iy,-1)),mag01(IAND(iy,INT(1,kind=8)))) |
---|
174 | mti = 0 |
---|
175 | ENDIF |
---|
176 | ! |
---|
177 | iy = mt(mti) |
---|
178 | mti = mti + 1 |
---|
179 | ! |
---|
180 | iy = IEOR(iy,ISHFT(iy,-11)) |
---|
181 | iy = IEOR(iy,IAND(ISHFT(iy,7),jptmp1)) |
---|
182 | iy = IEOR(iy,IAND(ISHFT(iy,15),jptmp2)) |
---|
183 | iy = IEOR(iy,ISHFT(iy,-18)) |
---|
184 | ! |
---|
185 | mtrand_int32 = iy |
---|
186 | END FUNCTION mtrand_int32 |
---|
187 | FUNCTION mtrand_int31() |
---|
188 | !----------------------------------------------------------------------- |
---|
189 | ! generates a random number on [0,0x7fffffff]-interval |
---|
190 | !----------------------------------------------------------------------- |
---|
191 | INTEGER :: mtrand_int31 |
---|
192 | mtrand_int31 = INT(ISHFT(mtrand_int32(),-1)) |
---|
193 | END FUNCTION mtrand_int31 |
---|
194 | FUNCTION mtrand_real1() |
---|
195 | !----------------------------------------------------------------------- |
---|
196 | ! generates a random number on [0,1]-real-interval |
---|
197 | !----------------------------------------------------------------------- |
---|
198 | REAL(wp) :: mtrand_real1, zr |
---|
199 | zr = REAL(mtrand_int32(),wp) |
---|
200 | IF (zr.LT.0._wp) zr = zr + 2._wp**32 |
---|
201 | mtrand_real1 = zr / 4294967295._wp |
---|
202 | END FUNCTION mtrand_real1 |
---|
203 | FUNCTION mtrand_real2() |
---|
204 | !----------------------------------------------------------------------- |
---|
205 | ! generates a random number on [0,1)-real-interval |
---|
206 | !----------------------------------------------------------------------- |
---|
207 | REAL(wp) :: mtrand_real2,zr |
---|
208 | zr = REAL(mtrand_int32(),wp) |
---|
209 | IF (zr.LT.0._wp) zr = zr + 2._wp**32 |
---|
210 | mtrand_real2=zr / 4294967296._wp |
---|
211 | END FUNCTION mtrand_real2 |
---|
212 | FUNCTION mtrand_real3() |
---|
213 | !----------------------------------------------------------------------- |
---|
214 | ! generates a random number on (0,1)-real-interval |
---|
215 | !----------------------------------------------------------------------- |
---|
216 | REAL(wp) :: mtrand_real3, zr |
---|
217 | zr = REAL(mtrand_int32(),wp) |
---|
218 | IF (zr.LT.0._wp) zr = zr + 2._wp**32 |
---|
219 | mtrand_real3 = ( zr + 0.5_wp ) / 4294967296._wp |
---|
220 | END FUNCTION mtrand_real3 |
---|
221 | FUNCTION mtrand_res53() |
---|
222 | !----------------------------------------------------------------------- |
---|
223 | ! generates a random number on [0,1) with 53-bit resolution |
---|
224 | !----------------------------------------------------------------------- |
---|
225 | REAL(wp) :: mtrand_res53 |
---|
226 | REAL(wp) :: za,zb |
---|
227 | za = ISHFT(mtrand_int32(),-5) |
---|
228 | zb = ISHFT(mtrand_int32(),-6) |
---|
229 | IF (za.LT.0._wp) za = za + 2._wp**32 |
---|
230 | IF (zb.LT.0._wp) zb = zb + 2._wp**32 |
---|
231 | mtrand_res53 = ( za * 67108864._wp + zb ) / 9007199254740992._wp |
---|
232 | END FUNCTION mtrand_res53 |
---|
233 | END MODULE mt19937ar |
---|