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 | *
|
---|
49 | * compile with: f90 -i32 -dp -o mt19937 mt19937.f
|
---|
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
|
---|
57 | *
|
---|
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 | *
|
---|
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
|
---|
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
|
---|
94 | integer n,n1,mt,mti
|
---|
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 | *
|
---|
106 | integer tmaskb, tmaskc, tshftu, tshfts, tshftt, tshftl, y
|
---|
107 |
|
---|
108 | integer n, n1, mti, mt, m, mata, lmask, mag01, kk
|
---|
109 |
|
---|
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 |
|
---|