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 | implicit integer(a-z)
|
---|
67 | *
|
---|
68 | * Period parameters
|
---|
69 | parameter(N = 624)
|
---|
70 | *
|
---|
71 | dimension mt(0:N-1)
|
---|
72 | * the array for the state vector
|
---|
73 | common /myblock/mti,mt
|
---|
74 | save /myblock/
|
---|
75 | *
|
---|
76 | * setting initial seeds to mt[N] using
|
---|
77 | * the generator Line 25 of Table 1 in
|
---|
78 | * [KNUTH 1981, The Art of Computer Programming
|
---|
79 | * Vol. 2 (2nd Ed.), pp102]
|
---|
80 | *
|
---|
81 | mt(0)= iand(seed,-1)
|
---|
82 | do 1000 mti=1,N-1
|
---|
83 | mt(mti) = iand(69069 * mt(mti-1),-1)
|
---|
84 | 1000 continue
|
---|
85 | *
|
---|
86 | return
|
---|
87 | end
|
---|
88 |
|
---|
89 | block data twbloks
|
---|
90 | parameter(N = 624)
|
---|
91 | parameter(N1 = N+1)
|
---|
92 | dimension mt(0:N-1)
|
---|
93 | common /myblock/mti,mt
|
---|
94 | save /myblock/
|
---|
95 | data mti/N1/
|
---|
96 | * mti==N+1 means mt[N] is not initialized
|
---|
97 | end
|
---|
98 | ************************************************************************
|
---|
99 | double precision function grnd()
|
---|
100 | *
|
---|
101 | implicit integer(a-z)
|
---|
102 | *
|
---|
103 | * Period parameters
|
---|
104 | parameter(N = 624)
|
---|
105 | parameter(N1 = N+1)
|
---|
106 | parameter(M = 397)
|
---|
107 | parameter(MATA = -1727483681)
|
---|
108 | * constant vector a
|
---|
109 | integer UMASK
|
---|
110 | * most significant w-r bits
|
---|
111 | parameter(LMASK = 2147483647)
|
---|
112 | * least significant r bits
|
---|
113 | * Tempering parameters
|
---|
114 | parameter(TMASKB= -1658038656)
|
---|
115 | parameter(TMASKC= -272236544)
|
---|
116 | *
|
---|
117 | dimension mt(0:N-1)
|
---|
118 | * the array for the state vector
|
---|
119 | common /myblock/mti,mt
|
---|
120 | save /myblock/
|
---|
121 | *
|
---|
122 | dimension mag01(0:1)
|
---|
123 | data mag01/0, MATA/
|
---|
124 | save mag01
|
---|
125 | * mag01(x) = x * MATA for x=0,1
|
---|
126 | *
|
---|
127 |
|
---|
128 | TSHFTU(y)=ishft(y,-11)
|
---|
129 | TSHFTS(y)=ishft(y,7)
|
---|
130 | TSHFTT(y)=ishft(y,15)
|
---|
131 | TSHFTL(y)=ishft(y,-18)
|
---|
132 | *
|
---|
133 | UMASK = -2147483647
|
---|
134 | UMASK = UMASK - 1
|
---|
135 |
|
---|
136 | if(mti.ge.N) then
|
---|
137 | * generate N words at one time
|
---|
138 | if(mti.eq.N+1) then
|
---|
139 | * if sgrnd() has not been called,
|
---|
140 | call sgrnd(4357)
|
---|
141 | * a default initial seed is used
|
---|
142 | endif
|
---|
143 | *
|
---|
144 | do 1000 kk=0,N-M-1
|
---|
145 | y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
|
---|
146 | mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
|
---|
147 | 1000 continue
|
---|
148 | do 1100 kk=N-M,N-2
|
---|
149 | y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
|
---|
150 | mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
|
---|
151 | 1100 continue
|
---|
152 | y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
|
---|
153 | mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
|
---|
154 | mti = 0
|
---|
155 | endif
|
---|
156 | *
|
---|
157 | y=mt(mti)
|
---|
158 | mti=mti+1
|
---|
159 | y=ieor(y,TSHFTU(y))
|
---|
160 | y=ieor(y,iand(TSHFTS(y),TMASKB))
|
---|
161 | y=ieor(y,iand(TSHFTT(y),TMASKC))
|
---|
162 | y=ieor(y,TSHFTL(y))
|
---|
163 | *
|
---|
164 | if(y.lt.0.0) then
|
---|
165 | grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0)
|
---|
166 | else
|
---|
167 | grnd=dble(y)/(2.0d0**32-1.0d0)
|
---|
168 | endif
|
---|
169 | *
|
---|
170 | return
|
---|
171 | end
|
---|
172 |
|
---|