[e40e335] | 1 | ************************************************************************
|
---|
| 2 | * A C-program for MT19937: Real number version
|
---|
| 3 | * genrand() generates one pseudorandom real number (double)
|
---|
| 4 | * which is uniformly distributed on [0,1]-interval, for each
|
---|
| 5 | * call. sgenrand(seed) set initial values to the working area
|
---|
| 6 | * of 624 words. Before genrand(), sgenrand(seed) must be
|
---|
| 7 | * called once. (seed is any 32-bit integer except for 0).
|
---|
| 8 | * Integer generator is obtained by modifying two lines.
|
---|
| 9 | * Coded by Takuji Nishimura, considering the suggestions by
|
---|
| 10 | * Topher Cooper and Marc Rieffel in July-Aug. 1997.
|
---|
| 11 | *
|
---|
| 12 | * This library is free software; you can redistribute it and/or
|
---|
| 13 | * modify it under the terms of the GNU Library General Public
|
---|
| 14 | * License as published by the Free Software Foundation; either
|
---|
| 15 | * version 2 of the License, or (at your option) any later
|
---|
| 16 | * version.
|
---|
| 17 | * This library is distributed in the hope that it will be useful,
|
---|
| 18 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 19 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
---|
| 20 | * See the GNU Library General Public License for more details.
|
---|
| 21 | * You should have received a copy of the GNU Library General
|
---|
| 22 | * Public License along with this library; if not, write to the
|
---|
| 23 | * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
|
---|
| 24 | * 02111-1307 USA
|
---|
| 25 | *
|
---|
| 26 | * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
|
---|
| 27 | * When you use this, send an email to: matumoto@math.keio.ac.jp
|
---|
| 28 | * with an appropriate reference to your work.
|
---|
| 29 | *
|
---|
| 30 | ************************************************************************
|
---|
| 31 | * Fortran translation by Hiroshi Takano. Jan. 13, 1999.
|
---|
| 32 | *
|
---|
| 33 | *
|
---|
| 34 | * genrand() -> double precision function grnd()
|
---|
| 35 | * sgenrand(seed) -> subroutine sgrnd(seed)
|
---|
| 36 | * integer seed
|
---|
| 37 | *
|
---|
| 38 | * This program uses the following non-standard intrinsics.
|
---|
| 39 | * ishft(i,n): If n>0, shifts bits in i by n positions to left.
|
---|
| 40 | * If n<0, shifts bits in i by n positions to right.
|
---|
| 41 | * iand (i,j): Performs logical AND on corresponding bits of i and j.
|
---|
| 42 | * ior (i,j): Performs inclusive OR on corresponding bits of i and j.
|
---|
| 43 | * ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
|
---|
| 44 | *
|
---|
| 45 | ************************************************************************
|
---|
| 46 | * Adapted to Cray T3E, Heinrich Stamerjohanns (stamer@uni-oldenburg.de)
|
---|
| 47 | * May 03, 2000
|
---|
| 48 | *
|
---|
[32289cd] | 49 | * compile with: f90 -i32 -dp -o mt19937 mt19937.f
|
---|
[e40e335] | 50 | * Changes:
|
---|
| 51 | * Do not use UMASK as parameter due to restrictions of CF90
|
---|
| 52 | * instead use UMASK = -MAX_INT_32, UMASK = UMASK - 1
|
---|
| 53 | *
|
---|
| 54 | * compare y to 0.0 due to unknown bug in CF90
|
---|
| 55 | *
|
---|
| 56 | * Now the T3E produces the same random number as a RS6000 running AIX
|
---|
[32289cd] | 57 | *
|
---|
[e40e335] | 58 | * from other programs sgrnd() must be called with integer*4
|
---|
| 59 | *
|
---|
| 60 | * All other programs do not need to be compiled with -i32 as long as
|
---|
| 61 | * you do not use integers as parameters to call functions compiled with -i32
|
---|
| 62 | *
|
---|
| 63 | ***********************************************************************
|
---|
| 64 | subroutine sgrnd(seed)
|
---|
| 65 | *
|
---|
[32289cd] | 66 | integer seed, tmaskb, tmaskc, tshftu, tshfts, tshftt
|
---|
| 67 | integer tshftl, y
|
---|
| 68 |
|
---|
| 69 | integer N, N1, mti, mt, m, mata, lmask, mag01
|
---|
| 70 | integer kk
|
---|
[e40e335] | 71 | *
|
---|
| 72 | * Period parameters
|
---|
| 73 | parameter(N = 624)
|
---|
| 74 | *
|
---|
| 75 | dimension mt(0:N-1)
|
---|
| 76 | * the array for the state vector
|
---|
| 77 | common /myblock/mti,mt
|
---|
| 78 | save /myblock/
|
---|
| 79 | *
|
---|
| 80 | * setting initial seeds to mt[N] using
|
---|
| 81 | * the generator Line 25 of Table 1 in
|
---|
| 82 | * [KNUTH 1981, The Art of Computer Programming
|
---|
| 83 | * Vol. 2 (2nd Ed.), pp102]
|
---|
| 84 | *
|
---|
| 85 | mt(0)= iand(seed,-1)
|
---|
| 86 | do 1000 mti=1,N-1
|
---|
| 87 | mt(mti) = iand(69069 * mt(mti-1),-1)
|
---|
| 88 | 1000 continue
|
---|
| 89 | *
|
---|
| 90 | return
|
---|
| 91 | end
|
---|
| 92 |
|
---|
| 93 | block data twbloks
|
---|
[32289cd] | 94 | integer n,n1,mt,mti
|
---|
[e40e335] | 95 | parameter(N = 624)
|
---|
| 96 | parameter(N1 = N+1)
|
---|
| 97 | dimension mt(0:N-1)
|
---|
| 98 | common /myblock/mti,mt
|
---|
| 99 | save /myblock/
|
---|
| 100 | data mti/N1/
|
---|
| 101 | * mti==N+1 means mt[N] is not initialized
|
---|
| 102 | end
|
---|
| 103 | ************************************************************************
|
---|
| 104 | double precision function grnd()
|
---|
| 105 | *
|
---|
[32289cd] | 106 | integer tmaskb, tmaskc, tshftu, tshfts, tshftt, tshftl, y
|
---|
| 107 |
|
---|
| 108 | integer n, n1, mti, mt, m, mata, lmask, mag01, kk
|
---|
| 109 |
|
---|
[e40e335] | 110 | *
|
---|
| 111 | * Period parameters
|
---|
| 112 | parameter(N = 624)
|
---|
| 113 | parameter(N1 = N+1)
|
---|
| 114 | parameter(M = 397)
|
---|
| 115 | parameter(MATA = -1727483681)
|
---|
| 116 | * constant vector a
|
---|
| 117 | integer UMASK
|
---|
| 118 | * most significant w-r bits
|
---|
| 119 | parameter(LMASK = 2147483647)
|
---|
| 120 | * least significant r bits
|
---|
| 121 | * Tempering parameters
|
---|
| 122 | parameter(TMASKB= -1658038656)
|
---|
| 123 | parameter(TMASKC= -272236544)
|
---|
| 124 | *
|
---|
| 125 | dimension mt(0:N-1)
|
---|
| 126 | * the array for the state vector
|
---|
| 127 | common /myblock/mti,mt
|
---|
| 128 | save /myblock/
|
---|
| 129 | *
|
---|
| 130 | dimension mag01(0:1)
|
---|
| 131 | data mag01/0, MATA/
|
---|
| 132 | save mag01
|
---|
| 133 | * mag01(x) = x * MATA for x=0,1
|
---|
| 134 | *
|
---|
| 135 |
|
---|
| 136 | TSHFTU(y)=ishft(y,-11)
|
---|
| 137 | TSHFTS(y)=ishft(y,7)
|
---|
| 138 | TSHFTT(y)=ishft(y,15)
|
---|
| 139 | TSHFTL(y)=ishft(y,-18)
|
---|
| 140 | *
|
---|
| 141 | UMASK = -2147483647
|
---|
| 142 | UMASK = UMASK - 1
|
---|
| 143 |
|
---|
| 144 | if(mti.ge.N) then
|
---|
| 145 | * generate N words at one time
|
---|
| 146 | if(mti.eq.N+1) then
|
---|
| 147 | * if sgrnd() has not been called,
|
---|
| 148 | call sgrnd(4357)
|
---|
| 149 | * a default initial seed is used
|
---|
| 150 | endif
|
---|
| 151 | *
|
---|
| 152 | do 1000 kk=0,N-M-1
|
---|
| 153 | y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
|
---|
| 154 | mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
|
---|
| 155 | 1000 continue
|
---|
| 156 | do 1100 kk=N-M,N-2
|
---|
| 157 | y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
|
---|
| 158 | mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
|
---|
| 159 | 1100 continue
|
---|
| 160 | y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
|
---|
| 161 | mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
|
---|
| 162 | mti = 0
|
---|
| 163 | endif
|
---|
| 164 | *
|
---|
| 165 | y=mt(mti)
|
---|
| 166 | mti=mti+1
|
---|
| 167 | y=ieor(y,TSHFTU(y))
|
---|
| 168 | y=ieor(y,iand(TSHFTS(y),TMASKB))
|
---|
| 169 | y=ieor(y,iand(TSHFTT(y),TMASKC))
|
---|
| 170 | y=ieor(y,TSHFTL(y))
|
---|
| 171 | *
|
---|
| 172 | if(y.lt.0.0) then
|
---|
| 173 | grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0)
|
---|
| 174 | else
|
---|
| 175 | grnd=dble(y)/(2.0d0**32-1.0d0)
|
---|
| 176 | endif
|
---|
| 177 | *
|
---|
| 178 | return
|
---|
| 179 | end
|
---|
| 180 |
|
---|