-
Notifications
You must be signed in to change notification settings - Fork 0
/
ziggurat.f90
333 lines (294 loc) · 9.37 KB
/
ziggurat.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
! Marsaglia & Tsang generator for random normals & random exponentials.
! Translated from C by Alan Miller ([email protected])
! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating
! random variables', J. Statist. Software, v5(8).
! This is an electronic journal which can be downloaded from:
! http://www.jstatsoft.org/v05/i08
! N.B. It is assumed that all integers are 32-bit.
! N.B. The value of M2 has been halved to compensate for the lack of
! unsigned integers in Fortran.
! Latest version - 1 January 2001
MODULE Ziggurat
IMPLICIT NONE
PRIVATE
INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND( 12, 60 )
REAL(DP), PARAMETER :: m1=2147483648.0_DP, m2=2147483648.0_DP, &
half=0.5_DP
REAL(DP) :: dn=3.442619855899_DP, tn=3.442619855899_DP, &
vn=0.00991256303526217_DP, &
q, de=7.697117470131487_DP, &
te=7.697117470131487_DP, &
ve=0.003949659822581572_DP
INTEGER, SAVE :: iz, jz, jsr=123456789, kn(0:127), &
ke(0:255), hz
REAL(DP), SAVE :: wn(0:127), fn(0:127), we(0:255), fe(0:255)
LOGICAL, SAVE :: initialized=.FALSE.
PUBLIC :: zigset, shr3, uni, rnor, rexp
CONTAINS
SUBROUTINE zigset( jsrseed )
INTEGER, INTENT(IN) :: jsrseed
INTEGER :: i
! Set the seed
jsr = jsrseed
! Tables for RNOR
q = vn*EXP(half*dn*dn)
kn(0) = (dn/q)*m1
kn(1) = 0
wn(0) = q/m1
wn(127) = dn/m1
fn(0) = 1.0_DP
fn(127) = EXP( -half*dn*dn )
DO i = 126, 1, -1
dn = SQRT( -2.0_DP * LOG( vn/dn + EXP( -half*dn*dn ) ) )
kn(i+1) = (dn/tn)*m1
tn = dn
fn(i) = EXP(-half*dn*dn)
wn(i) = dn/m1
END DO
! Tables for REXP
q = ve*EXP( de )
ke(0) = (de/q)*m2
ke(1) = 0
we(0) = q/m2
we(255) = de/m2
fe(0) = 1.0_DP
fe(255) = EXP( -de )
DO i = 254, 1, -1
de = -LOG( ve/de + EXP( -de ) )
ke(i+1) = m2 * (de/te)
te = de
fe(i) = EXP( -de )
we(i) = de/m2
END DO
initialized = .TRUE.
RETURN
END SUBROUTINE zigset
! Generate random 32-bit integers
FUNCTION shr3( ) RESULT( ival )
INTEGER :: ival
jz = jsr
jsr = IEOR( jsr, ISHFT( jsr, 13 ) )
jsr = IEOR( jsr, ISHFT( jsr, -17 ) )
jsr = IEOR( jsr, ISHFT( jsr, 5 ) )
ival = jz + jsr
RETURN
END FUNCTION shr3
! Generate uniformly distributed random numbers
FUNCTION uni( ) RESULT( fn_val )
REAL(DP) :: fn_val
fn_val = half + 0.2328306e-9_DP * shr3( )
RETURN
END FUNCTION uni
! Generate random normals
FUNCTION rnor( ) RESULT( fn_val )
REAL(DP) :: fn_val
REAL(DP), PARAMETER :: r = 3.442620_DP
REAL(DP) :: x, y
IF( .NOT. initialized ) CALL zigset( jsr )
hz = shr3( )
iz = IAND( hz, 127 )
IF( ABS( hz ) < kn(iz) ) THEN
fn_val = hz * wn(iz)
ELSE
DO
IF( iz == 0 ) THEN
DO
x = -0.2904764_DP* LOG( uni( ) )
y = -LOG( uni( ) )
IF( y+y >= x*x ) EXIT
END DO
fn_val = r+x
IF( hz <= 0 ) fn_val = -fn_val
RETURN
END IF
x = hz * wn(iz)
IF( fn(iz) + uni( )*(fn(iz-1)-fn(iz)) < EXP(-half*x*x) ) THEN
fn_val = x
RETURN
END IF
hz = shr3( )
iz = IAND( hz, 127 )
IF( ABS( hz ) < kn(iz) ) THEN
fn_val = hz * wn(iz)
RETURN
END IF
END DO
END IF
RETURN
END FUNCTION rnor
! Generate random exponentials
FUNCTION rexp( ) RESULT( fn_val )
REAL(DP) :: fn_val
REAL(DP) :: x
IF( .NOT. initialized ) CALL Zigset( jsr )
jz = shr3( )
iz = IAND( jz, 255 )
IF( ABS( jz ) < ke(iz) ) THEN
fn_val = ABS(jz) * we(iz)
RETURN
END IF
DO
IF( iz == 0 ) THEN
fn_val = 7.69711 - LOG( uni( ) )
RETURN
END IF
x = ABS( jz ) * we(iz)
IF( fe(iz) + uni( )*(fe(iz-1) - fe(iz)) < EXP( -x ) ) THEN
fn_val = x
RETURN
END IF
jz = shr3( )
iz = IAND( jz, 255 )
IF( ABS( jz ) < ke(iz) ) THEN
fn_val = ABS( jz ) * we(iz)
RETURN
END IF
END DO
RETURN
END FUNCTION rexp
END MODULE ziggurat
!PROGRAM test_ziggurat
! USE Ziggurat
! IMPLICIT NONE
!
! INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(12, 60)
! REAL :: t1, t2
! REAL(DP) :: x, chisq, expctd
! INTEGER :: i, ten_million=10000000, freq(0:999), j
!
! ! Set the random number seed
! WRITE(*, '(a)', ADVANCE='NO') ' Enter a non-zero integer: '
! READ(*, *) i
!
! CALL zigset( i )
!
! ! First time the generation of uniform, normal & exponential generators
! CALL CPU_TIME( t1 )
! DO i = 1, ten_million
! x = uni( )
! END DO
! CALL CPU_TIME( t2 )
! WRITE(*, '(a, i10, a, f9.2, a)') ' Time to generate ', ten_million, &
! ' uniform random nos. = ', t2-t1, 'sec.'
! CALL CPU_TIME( t1 )
! DO i = 1, ten_million
! x = rnor( )
! END DO
! CALL CPU_TIME( t2 )
! WRITE(*, '(a, i10, a, f9.2, a)') ' Time to generate ', ten_million, &
! ' normal random nos. = ', t2-t1, 'sec.'
!
! CALL CPU_TIME( t1 )
! DO i = 1, ten_million
! x = rexp( )
! END DO
! CALL CPU_TIME( t2 )
! WRITE(*, '(a, i10, a, f9.2, a)') ' Time to generate ', ten_million, &
! ' exponential random nos. = ', t2-t1, 'sec.'
! WRITE(*, *)
!
! ! Now look at the distributions generated
! WRITE(*, *) 'Std. devn. of chi-squared with 999 d.of.f. = 44.70'
! WRITE(*, *) 'Values of chi-squared below should be within about 90. of 999.'
!
! freq = 0
! DO i = 1, ten_million
! j = 1000 * uni( )
! freq(j) = freq(j) +1
! END DO
! chisq = 0.0_DP
! expctd = ten_million / 1000
! DO j = 0, 999
! chisq = chisq + (freq(j) - expctd)**2 / expctd
! END DO
! WRITE(*, '(a, f9.2, a)') ' Uniform distribution, chi-squared = ', chisq, &
! ' with 999 deg. of freedom'
! freq = 0
! DO i = 1, ten_million
! j = 1000 * alnorm( rnor( ), .FALSE. )
! freq(j) = freq(j) +1
! END DO
! chisq = 0.0_DP
! expctd = ten_million / 1000
! DO j = 0, 999
! chisq = chisq + (freq(j) - expctd)**2 / expctd
! END DO
! WRITE(*, '(a, f9.2, a)') ' Normal distribution, chi-squared = ', chisq, &
! ' with 999 deg. of freedom'
! freq = 0
! DO i = 1, ten_million
! j = 1000 * EXP( -rexp( ) )
! IF( j > 999 .OR. j < 0 ) THEN
! WRITE(*, '(a, 2i10)') ' i, j = ', i, j
! ELSE
! freq(j) = freq(j) +1
! END IF
! END DO
! chisq = 0.0_DP
! expctd = ten_million / 1000
! DO j = 0, 999
! chisq = chisq + (freq(j) - expctd)**2 / expctd
! END DO
! WRITE(*, '(a, f9.2, a)') ' Exponential distribution, chi-squared = ', &
! chisq, ' with 999 deg. of freedom'
! STOP
!
!
!
!CONTAINS
!
!
!! Algorithm AS66 Applied Statistics (1973) vol.22, no.3
!
!! Evaluates the tail area of the standardised normal curve
!! from x to infinity if upper is .true. or
!! from minus infinity to x if upper is .false.
!
!! ELF90-compatible version by Alan Miller
!! Latest revision - 29 November 2001
!
!FUNCTION alnorm( x, upper ) RESULT( fn_val )
! REAL(DP), INTENT(IN) :: x
! LOGICAL, INTENT(IN) :: upper
! REAL(DP) :: fn_val
!
! ! Local variables
! REAL(DP), PARAMETER :: zero=0.0_DP, one=1.0_DP, half=0.5_DP, con=1.28_DP
! REAL(DP) :: z, y
! LOGICAL :: up
!
! ! Machine dependent constants
! REAL(DP), PARAMETER :: ltone = 7.0_DP, utzero = 18.66_DP
! REAL(DP), PARAMETER :: p = 0.398942280444_DP, q = 0.39990348504_DP, &
! r = 0.398942280385_DP, a1 = 5.75885480458_DP, &
! a2 = 2.62433121679_DP, a3 = 5.92885724438_DP, &
! b1 = -29.8213557807_DP, b2 = 48.6959930692_DP, &
! c1 = -3.8052E-8_DP, c2 = 3.98064794E-4_DP, &
! c3 = -0.151679116635_DP, c4 = 4.8385912808_DP, &
! c5 = 0.742380924027_DP, c6 = 3.99019417011_DP, &
! d1 = 1.00000615302_DP, d2 = 1.98615381364_DP, &
! d3 = 5.29330324926_DP, d4 = -15.1508972451_DP, &
! d5 = 30.789933034_DP
!
! up = upper
! z = x
! IF( z < zero ) THEN
! up = .NOT. up
! z = -z
! END IF
! IF( z <= ltone .OR. (up .AND. z <= utzero) ) THEN
! y = half*z*z
! IF( z > con ) THEN
! fn_val = r*EXP( -y )/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6))))))
! ELSE
! fn_val = half - z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3))))
! END IF
! ELSE
! fn_val = zero
! END IF
!
! IF( .NOT. up ) fn_val = one - fn_val
! RETURN
!END FUNCTION alnorm
!
!END PROGRAM test_ziggurat