[3611] | 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 |
---|