1 |
#Region "Microsoft.VisualBasic::120bc0a3db9510757289f89883c14448, Microsoft.VisualBasic.Core\Language\Language\Java\MersenneTwisterFast.vb"
|
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 |
#End Region
|
50 |
|
51 |
Imports sys = System.Math
|
52 |
|
53 |
|
54 |
|
55 |
Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut
|
56 |
|
57 |
This file is part of BEAST.
|
58 |
See the NOTICE file distributed with this work for additional
|
59 |
information regarding copyright ownership and licensing.
|
60 |
|
61 |
BEAST is free software; you can redistribute it and/or modify
|
62 |
it under the terms of the GNU Lesser General Public License as
|
63 |
published by the Free Software Foundation; either version 2
|
64 |
of the License, or (at your option) any later version.
|
65 |
|
66 |
BEAST is distributed in the hope that it will be useful,
|
67 |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
68 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
69 |
GNU Lesser General Public License for more details.
|
70 |
|
71 |
You should have received a copy of the GNU Lesser General Public
|
72 |
License along with BEAST; if not, write to the
|
73 |
Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
|
74 |
Boston, MA 02110-1301 USA
|
75 |
|
76 |
|
77 |
Namespace Language.Java
|
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 |
<Serializable>
|
127 |
Public Class MersenneTwisterFast
|
128 |
|
129 |
Private Const serialVersionUID As Long = 6185086957226269797L
|
130 |
Period parameters
|
131 |
Private Const N As Integer = 624
|
132 |
Private Const M As Integer = 397
|
133 |
Private Const MATRIX_A As Integer = &H9908B0DF private static final *
|
134 |
constant vector a
|
135 |
Private Const UPPER_MASK As Integer = &H80000000 most significant w-r
|
136 |
bits
|
137 |
Private Const LOWER_MASK As Integer = &H7FFFFFFF least significant r
|
138 |
bits
|
139 |
|
140 |
Tempering parameters
|
141 |
Private Const TEMPERING_MASK_B As Integer = &H9D2C5680
|
142 |
Private Const TEMPERING_MASK_C As Integer = &HEFC60000
|
143 |
|
144 |
#define TEMPERING_SHIFT_U(y) (y >>> 11)
|
145 |
#define TEMPERING_SHIFT_S(y) (y << 7)
|
146 |
#define TEMPERING_SHIFT_T(y) (y << 15)
|
147 |
#define TEMPERING_SHIFT_L(y) (y >>> 18)
|
148 |
|
149 |
Private mt As Integer() the array for the state vector
|
150 |
Private mti As Integer mti==N+1 means mt[N] is not initialized
|
151 |
Private mag01 As Integer()
|
152 |
|
153 |
a good initial seed (of int size, though stored in a long)
|
154 |
Private Const GOOD_SEED As Long = 4357
|
155 |
|
156 |
Private nextNextGaussian As Double
|
157 |
Private haveNextNextGaussian As Boolean
|
158 |
|
159 |
The following can be accessed externally by the static accessor methods
|
160 |
which
|
161 |
inforce synchronization
|
162 |
Public Shared ReadOnly DEFAULT_INSTANCE As New MersenneTwisterFast
|
163 |
|
164 |
Added to curernt time in default constructor, and then adjust to allow
|
165 |
for programs that construct
|
166 |
multiple MersenneTwisterFast in a short amount of time.
|
167 |
Private Shared seedAdditive_ As Long = 0
|
168 |
|
169 |
Private initializationSeed As Long
|
170 |
|
171 |
|
172 |
Constructor using the time of day as default seed.
|
173 |
|
174 |
Public Sub New()
|
175 |
Me.New(Now.ToBinary + seedAdditive_)
|
176 |
seedAdditive_ += nextInt()
|
177 |
End Sub
|
178 |
|
179 |
|
180 |
Constructor using a given seed. Though you pass this seed in as a long,
|
181 |
it
|
182 |
|
183 |
<param name="seed">
|
184 |
generator starting number, often the time of day. </param>
|
185 |
Private Sub New(seed As Long)
|
186 |
If seed = 0 Then
|
187 |
seed = GOOD_SEED
|
188 |
Else
|
189 |
seed = seed
|
190 |
End If
|
191 |
End Sub
|
192 |
|
193 |
|
194 |
Initalize the pseudo random number generator. The Mersenne Twister only
|
195 |
uses an integer for its seed; It
|
196 |
that
|
197 |
|
198 |
Public Property Seed As Long
|
199 |
Set(seed As Long)
|
200 |
If seed = 0 Then Throw New System.ArgumentException("Non zero random seed required.")
|
201 |
initializationSeed = seed
|
202 |
haveNextNextGaussian = False
|
203 |
|
204 |
mt = New Integer(N - 1) {}
|
205 |
|
206 |
setting initial seeds to mt[N] using
|
207 |
the generator Line 25 of Table 1 in
|
208 |
[KNUTH 1981, The Art of Computer Programming
|
209 |
Vol. 2 (2nd Ed.), pp102]
|
210 |
|
211 |
the 0xffffffff is commented out because in Java
|
212 |
ints are always 32 bits; hence i & 0xffffffff == i
|
213 |
|
214 |
mt(0) = (CInt(seed)) & 0xffffffff;
|
215 |
|
216 |
For mti = 1 To N - 1
|
217 |
mt(mti) = (69069 * mt(mti - 1)) & 0xffffffff;
|
218 |
Next mti
|
219 |
|
220 |
mag01[x] = x * MATRIX_A for x=0,1
|
221 |
mag01 = New Integer(1) {}
|
222 |
mag01(0) = &H0
|
223 |
mag01(1) = MATRIX_A
|
224 |
End Set
|
225 |
Get
|
226 |
Return initializationSeed
|
227 |
End Get
|
228 |
End Property
|
229 |
|
230 |
|
231 |
Public Function nextInt() As Integer
|
232 |
Dim y As Integer
|
233 |
|
234 |
If mti >= N Then generate N words at one time
|
235 |
Dim kk As Integer
|
236 |
|
237 |
For kk = 0 To N - M - 1
|
238 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
239 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
240 |
Next kk
|
241 |
Do While kk < N - 1
|
242 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
243 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
244 |
kk += 1
|
245 |
Loop
|
246 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
247 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
248 |
|
249 |
mti = 0
|
250 |
End If
|
251 |
|
252 |
y = mt(mti)
|
253 |
mti += 1
|
254 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
255 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
256 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
257 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
258 |
|
259 |
Return y
|
260 |
End Function
|
261 |
|
262 |
Public Function nextShort() As Short
|
263 |
Dim y As Integer
|
264 |
|
265 |
If mti >= N Then generate N words at one time
|
266 |
Dim kk As Integer
|
267 |
|
268 |
For kk = 0 To N - M - 1
|
269 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
270 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
271 |
Next kk
|
272 |
Do While kk < N - 1
|
273 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
274 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
275 |
kk += 1
|
276 |
Loop
|
277 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
278 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
279 |
|
280 |
mti = 0
|
281 |
End If
|
282 |
|
283 |
y = mt(mti)
|
284 |
mti += 1
|
285 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
286 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
287 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
288 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
289 |
|
290 |
Return CShort(CInt(CUInt(y) >> 16))
|
291 |
End Function
|
292 |
|
293 |
Public Function nextChar() As Char
|
294 |
Dim y As Integer
|
295 |
|
296 |
If mti >= N Then generate N words at one time
|
297 |
Dim kk As Integer
|
298 |
|
299 |
For kk = 0 To N - M - 1
|
300 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
301 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
302 |
Next kk
|
303 |
Do While kk < N - 1
|
304 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
305 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
306 |
kk += 1
|
307 |
Loop
|
308 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
309 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
310 |
|
311 |
mti = 0
|
312 |
End If
|
313 |
|
314 |
y = mt(mti)
|
315 |
mti += 1
|
316 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
317 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
318 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
319 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
320 |
|
321 |
Return ChrW(CInt(CUInt(y) >> 16))
|
322 |
End Function
|
323 |
|
324 |
Public Function nextBoolean() As Boolean
|
325 |
Dim y As Integer
|
326 |
|
327 |
If mti >= N Then generate N words at one time
|
328 |
Dim kk As Integer
|
329 |
|
330 |
For kk = 0 To N - M - 1
|
331 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
332 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
333 |
Next kk
|
334 |
Do While kk < N - 1
|
335 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
336 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
337 |
kk += 1
|
338 |
Loop
|
339 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
340 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
341 |
|
342 |
mti = 0
|
343 |
End If
|
344 |
|
345 |
y = mt(mti)
|
346 |
mti += 1
|
347 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
348 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
349 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
350 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
351 |
|
352 |
Return ((CInt(CUInt(y) >> 31)) <> 0)
|
353 |
End Function
|
354 |
|
355 |
Public Function nextByte() As SByte
|
356 |
Dim y As Integer
|
357 |
|
358 |
If mti >= N Then generate N words at one time
|
359 |
Dim kk As Integer
|
360 |
|
361 |
For kk = 0 To N - M - 1
|
362 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
363 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
364 |
Next kk
|
365 |
Do While kk < N - 1
|
366 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
367 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
368 |
kk += 1
|
369 |
Loop
|
370 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
371 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
372 |
|
373 |
mti = 0
|
374 |
End If
|
375 |
|
376 |
y = mt(mti)
|
377 |
mti += 1
|
378 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
379 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
380 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
381 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
382 |
|
383 |
Return CByte(CInt(CUInt(y) >> 24))
|
384 |
End Function
|
385 |
|
386 |
Public Sub nextBytes(bytes As SByte())
|
387 |
Dim y As Integer
|
388 |
|
389 |
For x As Integer = 0 To bytes.Length - 1
|
390 |
If mti >= N Then generate N words at one time
|
391 |
Dim kk As Integer
|
392 |
|
393 |
For kk = 0 To N - M - 1
|
394 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
395 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
396 |
Next kk
|
397 |
Do While kk < N - 1
|
398 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
399 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
400 |
kk += 1
|
401 |
Loop
|
402 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
403 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
404 |
|
405 |
mti = 0
|
406 |
End If
|
407 |
|
408 |
y = mt(mti)
|
409 |
mti += 1
|
410 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
411 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
412 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
413 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
414 |
|
415 |
bytes(x) = CByte(CInt(CUInt(y) >> 24))
|
416 |
Next x
|
417 |
End Sub
|
418 |
|
419 |
Public Function nextLong() As Long
|
420 |
Dim y As Integer
|
421 |
Dim z As Integer
|
422 |
|
423 |
If mti >= N Then generate N words at one time
|
424 |
Dim kk As Integer
|
425 |
|
426 |
For kk = 0 To N - M - 1
|
427 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
428 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
429 |
Next kk
|
430 |
Do While kk < N - 1
|
431 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
432 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
433 |
kk += 1
|
434 |
Loop
|
435 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
436 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
437 |
|
438 |
mti = 0
|
439 |
End If
|
440 |
|
441 |
y = mt(mti)
|
442 |
mti += 1
|
443 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
444 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
445 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
446 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
447 |
|
448 |
If mti >= N Then generate N words at one time
|
449 |
Dim kk As Integer
|
450 |
|
451 |
For kk = 0 To N - M - 1
|
452 |
z = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
453 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
454 |
Next kk
|
455 |
Do While kk < N - 1
|
456 |
z = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
457 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
458 |
kk += 1
|
459 |
Loop
|
460 |
z = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
461 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
462 |
|
463 |
mti = 0
|
464 |
End If
|
465 |
|
466 |
z = mt(mti)
|
467 |
mti += 1
|
468 |
z = z Xor CInt(CUInt(z) >> 11) TEMPERING_SHIFT_U(z)
|
469 |
z = z Xor (z << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(z)
|
470 |
z = z Xor (z << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(z)
|
471 |
z = z Xor (CInt(CUInt(z) >> 18)) TEMPERING_SHIFT_L(z)
|
472 |
|
473 |
Return ((CLng(y)) << 32) + CLng(z)
|
474 |
End Function
|
475 |
|
476 |
Public Function nextDouble() As Double
|
477 |
Dim y As Integer
|
478 |
Dim z As Integer
|
479 |
|
480 |
If mti >= N Then generate N words at one time
|
481 |
Dim kk As Integer
|
482 |
|
483 |
For kk = 0 To N - M - 1
|
484 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
485 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
486 |
Next kk
|
487 |
Do While kk < N - 1
|
488 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
489 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
490 |
kk += 1
|
491 |
Loop
|
492 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
493 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
494 |
|
495 |
mti = 0
|
496 |
End If
|
497 |
|
498 |
y = mt(mti)
|
499 |
mti += 1
|
500 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
501 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
502 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
503 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
504 |
|
505 |
If mti >= N Then generate N words at one time
|
506 |
Dim kk As Integer
|
507 |
|
508 |
For kk = 0 To N - M - 1
|
509 |
z = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
510 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
511 |
Next kk
|
512 |
Do While kk < N - 1
|
513 |
z = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
514 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
515 |
kk += 1
|
516 |
Loop
|
517 |
z = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
518 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
519 |
|
520 |
mti = 0
|
521 |
End If
|
522 |
|
523 |
z = mt(mti)
|
524 |
mti += 1
|
525 |
z = z Xor CInt(CUInt(z) >> 11) TEMPERING_SHIFT_U(z)
|
526 |
z = z Xor (z << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(z)
|
527 |
z = z Xor (z << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(z)
|
528 |
z = z Xor (CInt(CUInt(z) >> 18)) TEMPERING_SHIFT_L(z)
|
529 |
|
530 |
derived from nextDouble documentation in jdk 1.2 docs, see top
|
531 |
Return (((CLng(CInt(CUInt(y) >> 6))) << 27) + (CInt(CUInt(z) >> 5))) / CDbl(1L << 53)
|
532 |
End Function
|
533 |
|
534 |
Public Function nextGaussian() As Double
|
535 |
If haveNextNextGaussian Then
|
536 |
haveNextNextGaussian = False
|
537 |
Return nextNextGaussian
|
538 |
Else
|
539 |
Dim v1, v2, s As Double
|
540 |
Do
|
541 |
Dim y As Integer
|
542 |
Dim z As Integer
|
543 |
Dim a As Integer
|
544 |
Dim b As Integer
|
545 |
|
546 |
If mti >= N Then generate N words at one time
|
547 |
Dim kk As Integer
|
548 |
|
549 |
For kk = 0 To N - M - 1
|
550 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
551 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
552 |
Next kk
|
553 |
Do While kk < N - 1
|
554 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
555 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
556 |
kk += 1
|
557 |
Loop
|
558 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
559 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
560 |
|
561 |
mti = 0
|
562 |
End If
|
563 |
|
564 |
y = mt(mti)
|
565 |
mti += 1
|
566 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
567 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
568 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
569 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
570 |
|
571 |
If mti >= N Then generate N words at one time
|
572 |
Dim kk As Integer
|
573 |
|
574 |
For kk = 0 To N - M - 1
|
575 |
z = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
576 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
577 |
Next kk
|
578 |
Do While kk < N - 1
|
579 |
z = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
580 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
581 |
kk += 1
|
582 |
Loop
|
583 |
z = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
584 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(z) >> 1)) Xor mag01(z And &H1)
|
585 |
|
586 |
mti = 0
|
587 |
End If
|
588 |
|
589 |
z = mt(mti)
|
590 |
mti += 1
|
591 |
z = z Xor CInt(CUInt(z) >> 11) TEMPERING_SHIFT_U(z)
|
592 |
z = z Xor (z << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(z)
|
593 |
z = z Xor (z << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(z)
|
594 |
z = z Xor (CInt(CUInt(z) >> 18)) TEMPERING_SHIFT_L(z)
|
595 |
|
596 |
If mti >= N Then generate N words at one time
|
597 |
Dim kk As Integer
|
598 |
|
599 |
For kk = 0 To N - M - 1
|
600 |
a = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
601 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(a) >> 1)) Xor mag01(a And &H1)
|
602 |
Next kk
|
603 |
Do While kk < N - 1
|
604 |
a = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
605 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(a) >> 1)) Xor mag01(a And &H1)
|
606 |
kk += 1
|
607 |
Loop
|
608 |
a = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
609 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(a) >> 1)) Xor mag01(a And &H1)
|
610 |
|
611 |
mti = 0
|
612 |
End If
|
613 |
|
614 |
a = mt(mti)
|
615 |
mti += 1
|
616 |
a = a Xor CInt(CUInt(a) >> 11) TEMPERING_SHIFT_U(a)
|
617 |
a = a Xor (a << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(a)
|
618 |
a = a Xor (a << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(a)
|
619 |
a = a Xor (CInt(CUInt(a) >> 18)) TEMPERING_SHIFT_L(a)
|
620 |
|
621 |
If mti >= N Then generate N words at one time
|
622 |
Dim kk As Integer
|
623 |
|
624 |
For kk = 0 To N - M - 1
|
625 |
b = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
626 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(b) >> 1)) Xor mag01(b And &H1)
|
627 |
Next kk
|
628 |
Do While kk < N - 1
|
629 |
b = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
630 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(b) >> 1)) Xor mag01(b And &H1)
|
631 |
kk += 1
|
632 |
Loop
|
633 |
b = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
634 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(b) >> 1)) Xor mag01(b And &H1)
|
635 |
|
636 |
mti = 0
|
637 |
End If
|
638 |
|
639 |
b = mt(mti)
|
640 |
mti += 1
|
641 |
b = b Xor CInt(CUInt(b) >> 11) TEMPERING_SHIFT_U(b)
|
642 |
b = b Xor (b << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(b)
|
643 |
b = b Xor (b << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(b)
|
644 |
b = b Xor (CInt(CUInt(b) >> 18)) TEMPERING_SHIFT_L(b)
|
645 |
|
646 |
|
647 |
* derived from nextDouble documentation in jdk 1.2 docs, see
|
648 |
* top
|
649 |
|
650 |
v1 = 2 * ((((CLng(CInt(CUInt(y) >> 6))) << 27) + (CInt(CUInt(z) >> 5))) / CDbl(1L << 53)) - 1
|
651 |
v2 = 2 * ((((CLng(CInt(CUInt(a) >> 6))) << 27) + (CInt(CUInt(b) >> 5))) / CDbl(1L << 53)) - 1
|
652 |
s = v1 * v1 + v2 * v2
|
653 |
Loop While s >= 1
|
654 |
Dim multiplier As Double = sys.Sqrt(-2 * sys.Log(s) / s)
|
655 |
nextNextGaussian = v2 * multiplier
|
656 |
haveNextNextGaussian = True
|
657 |
Return v1 * multiplier
|
658 |
End If
|
659 |
End Function
|
660 |
|
661 |
Public Function nextFloat() As Single
|
662 |
Dim y As Integer
|
663 |
|
664 |
If mti >= N Then generate N words at one time
|
665 |
Dim kk As Integer
|
666 |
|
667 |
For kk = 0 To N - M - 1
|
668 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
669 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
670 |
Next kk
|
671 |
Do While kk < N - 1
|
672 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
673 |
mt(kk) = mt(kk + (M - N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
674 |
kk += 1
|
675 |
Loop
|
676 |
y = (mt(N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
677 |
mt(N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
678 |
|
679 |
mti = 0
|
680 |
End If
|
681 |
|
682 |
y = mt(mti)
|
683 |
mti += 1
|
684 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
685 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
686 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
687 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
688 |
|
689 |
Return (CInt(CUInt(y) >> 8)) / (CSng(1 << 24))
|
690 |
End Function
|
691 |
|
692 |
|
693 |
Returns an integer drawn uniformly from 0 to n-1. Suffice it to say, n
|
694 |
must be > 0, or an IllegalArgumentException is raised.
|
695 |
|
696 |
Public Overridable Function nextInt(n As Integer) As Integer
|
697 |
If n <= 0 Then Throw New System.ArgumentException("n must be positive")
|
698 |
|
699 |
If (n And -n) = n Then i.e., n is a power of 2
|
700 |
Dim y As Integer
|
701 |
|
702 |
If mti >= MersenneTwisterFast.N Then generate N words at one time
|
703 |
Dim kk As Integer
|
704 |
|
705 |
For kk = 0 To MersenneTwisterFast.N - M - 1
|
706 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
707 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
708 |
Next kk
|
709 |
Do While kk < MersenneTwisterFast.N - 1
|
710 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
711 |
mt(kk) = mt(kk + (M - MersenneTwisterFast.N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
712 |
kk += 1
|
713 |
Loop
|
714 |
y = (mt(MersenneTwisterFast.N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
715 |
mt(MersenneTwisterFast.N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
716 |
|
717 |
mti = 0
|
718 |
End If
|
719 |
|
720 |
y = mt(mti)
|
721 |
mti += 1
|
722 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
723 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
724 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
725 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
726 |
|
727 |
Return CInt(Fix((n * CLng(CInt(CUInt(y) >> 1))) >> 31))
|
728 |
End If
|
729 |
|
730 |
Dim bits, val As Integer
|
731 |
Do
|
732 |
Dim y As Integer
|
733 |
|
734 |
If mti >= MersenneTwisterFast.N Then generate N words at one time
|
735 |
Dim kk As Integer
|
736 |
|
737 |
For kk = 0 To MersenneTwisterFast.N - M - 1
|
738 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
739 |
mt(kk) = mt(kk + M) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
740 |
Next kk
|
741 |
Do While kk < MersenneTwisterFast.N - 1
|
742 |
y = (mt(kk) And UPPER_MASK) Or (mt(kk + 1) And LOWER_MASK)
|
743 |
mt(kk) = mt(kk + (M - MersenneTwisterFast.N)) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
744 |
kk += 1
|
745 |
Loop
|
746 |
y = (mt(MersenneTwisterFast.N - 1) And UPPER_MASK) Or (mt(0) And LOWER_MASK)
|
747 |
mt(MersenneTwisterFast.N - 1) = mt(M - 1) Xor (CInt(CUInt(y) >> 1)) Xor mag01(y And &H1)
|
748 |
|
749 |
mti = 0
|
750 |
End If
|
751 |
|
752 |
y = mt(mti)
|
753 |
mti += 1
|
754 |
y = y Xor CInt(CUInt(y) >> 11) TEMPERING_SHIFT_U(y)
|
755 |
y = y Xor (y << 7) And TEMPERING_MASK_B TEMPERING_SHIFT_S(y)
|
756 |
y = y Xor (y << 15) And TEMPERING_MASK_C TEMPERING_SHIFT_T(y)
|
757 |
y = y Xor (CInt(CUInt(y) >> 18)) TEMPERING_SHIFT_L(y)
|
758 |
|
759 |
bits = (CInt(CUInt(y) >> 1))
|
760 |
val = bits Mod n
|
761 |
Loop While bits - val + (n - 1) < 0
|
762 |
Return val
|
763 |
End Function
|
764 |
|
765 |
|
766 |
Returns a uniform random permutation of int objects in array
|
767 |
|
768 |
Public Sub permute(array As Integer())
|
769 |
Dim l As Integer = array.Length
|
770 |
For i As Integer = 0 To l - 1
|
771 |
Dim index As Integer = nextInt(l - i) + i
|
772 |
Dim temp As Integer = array(index)
|
773 |
array(index) = array(i)
|
774 |
array(i) = temp
|
775 |
Next i
|
776 |
End Sub
|
777 |
|
778 |
|
779 |
Shuffles an array.
|
780 |
|
781 |
Public Sub shuffle(array As Integer())
|
782 |
Dim l As Integer = array.Length
|
783 |
For i As Integer = 0 To l - 1
|
784 |
Dim index As Integer = nextInt(l - i) + i
|
785 |
Dim temp As Integer = array(index)
|
786 |
array(index) = array(i)
|
787 |
array(i) = temp
|
788 |
Next i
|
789 |
End Sub
|
790 |
|
791 |
|
792 |
Shuffles an array. Shuffles numberOfShuffles times
|
793 |
|
794 |
Public Sub shuffle(array As Integer(), numberOfShuffles As Integer)
|
795 |
Dim i As Integer, j As Integer, temp As Integer, l As Integer = array.Length
|
796 |
For shuffle As Integer = 0 To numberOfShuffles - 1
|
797 |
Do
|
798 |
i = nextInt(l)
|
799 |
j = nextInt(l)
|
800 |
Loop While i <> j
|
801 |
temp = array(j)
|
802 |
array(j) = array(i)
|
803 |
array(i) = temp
|
804 |
Next shuffle
|
805 |
End Sub
|
806 |
|
807 |
|
808 |
Returns an array of shuffled indices of length l.
|
809 |
|
810 |
<param name="l">
|
811 |
length of the array required. </param>
|
812 |
Public Overridable Function shuffled(l As Integer) As Integer()
|
813 |
|
814 |
Dim array As Integer() = New Integer(l - 1) {}
|
815 |
|
816 |
initialize array
|
817 |
For i As Integer = 0 To l - 1
|
818 |
array(i) = i
|
819 |
Next i
|
820 |
shuffle(array)
|
821 |
|
822 |
Return array
|
823 |
End Function
|
824 |
|
825 |
|
826 |
Returns a uniform random permutation of ints 0,...,l-1
|
827 |
|
828 |
<param name="l">
|
829 |
length of the array required. </param>
|
830 |
Public Overridable Function permuted(l As Integer) As Integer()
|
831 |
|
832 |
Dim array As Integer() = New Integer(l - 1) {}
|
833 |
|
834 |
initialize array
|
835 |
For i As Integer = 0 To l - 1
|
836 |
array(i) = i
|
837 |
Next i
|
838 |
permute(array)
|
839 |
|
840 |
Return array
|
841 |
End Function
|
842 |
|
843 |
|
844 |
|
845 |
* Gamma Distribution - Acceptance Rejection combined with *
|
846 |
Acceptance Complement * *
|
847 |
******************************************************************
|
848 |
* FUNCTION: - gds samples a random number from the standard * gamma
|
849 |
distribution with parameter a > 0. * Acceptance Rejection gs for a <
|
850 |
1 , * Acceptance Complement gd for a >= 1 . * REFERENCES: - J.H.
|
851 |
Ahrens, U. Dieter (1974): Computer methods * for sampling from gamma,
|
852 |
beta, Poisson and * binomial distributions, Computing 12, 223-246. *
|
853 |
- J.H. Ahrens, U. Dieter (1982): Generating gamma * variates by a
|
854 |
modified rejection technique, * Communications of the ACM 25, 47-54.
|
855 |
* SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with *
|
856 |
unsigned long integer *seed * - NORMAL(seed) ... Normal generator
|
857 |
N(0,1). * *
|
858 |
*****************************************************************
|
859 |
|
860 |
Public Overridable Function nextGamma(alpha As Double, lambda As Double) As Double
|
861 |
Dim a As Double = alpha
|
862 |
Dim aa As Double = -1.0, aaa As Double = -1.0, b As Double = 0.0, c As Double = 0.0, d As Double = 0.0, e As Double, r As Double, s As Double = 0.0, si As Double = 0.0, ss As Double = 0.0, q0 As Double = 0.0, q1 As Double = 0.0416666664, q2 As Double = 0.0208333723, q3 As Double = 0.0079849875, q4 As Double = 0.0015746717, q5 As Double = -0.0003349403, q6 As Double = 0.0003340332, q7 As Double = 0.0006053049, q8 As Double = -0.0004701849, q9 As Double = 0.000171032, a1 As Double = 0.333333333, a2 As Double = -0.249999949, a3 As Double = 0.199999867, a4 As Double = -0.166677482, a5 As Double = 0.142873973, a6 As Double = -0.124385581, a7 As Double = 0.11036831, a8 As Double = -0.112750886, a9 As Double = 0.104089866, e1 As Double = 1.0, e2 As Double = 0.499999994, e3 As Double = 0.166666848, e4 As Double = 0.041664508, e5 As Double = 0.008345522, e6 As Double = 0.001353826, e7 As Double = 0.000247453
|
863 |
|
864 |
Dim gds, p, q, t, sign_u, u, v, w, x As Double
|
865 |
Dim v1, v2, v12 As Double
|
866 |
|
867 |
Check for invalid input values
|
868 |
|
869 |
If a <= 0.0 Then Throw New System.ArgumentException
|
870 |
If lambda <= 0.0 Then Dim TempIllegalArgumentException As System.ArgumentException = New System.ArgumentException
|
871 |
|
872 |
If a < 1.0 Then CASE A: Acceptance rejection algorithm gs
|
873 |
b = 1.0 + 0.36788794412 * a Step 1
|
874 |
Do
|
875 |
p = b * nextDouble()
|
876 |
If p <= 1.0 Then Step 2. Case gds <= 1
|
877 |
gds = sys.Exp(sys.Log(p) / a)
|
878 |
If sys.Log(nextDouble()) <= -gds Then
|
879 |
Return (gds / lambda)
|
880 |
End If Step 3. Case gds > 1
|
881 |
Else
|
882 |
gds = -sys.Log((b - p) / a)
|
883 |
If sys.Log(nextDouble()) <= ((a - 1.0) * sys.Log(gds)) Then Return (gds / lambda)
|
884 |
End If
|
885 |
Loop CASE B: Acceptance complement algorithm gd (gaussian
|
886 |
Else
|
887 |
distribution, box muller transformation)
|
888 |
If a <> aa Then Step 1. Preparations
|
889 |
aa = a
|
890 |
ss = a - 0.5
|
891 |
s = sys.Sqrt(ss)
|
892 |
d = 5.656854249 - 12.0 * s
|
893 |
End If
|
894 |
Step 2. Normal deviate
|
895 |
Do
|
896 |
v1 = 2.0 * nextDouble() - 1.0
|
897 |
v2 = 2.0 * nextDouble() - 1.0
|
898 |
v12 = v1 * v1 + v2 * v2
|
899 |
Loop While v12 > 1.0
|
900 |
t = v1 * sys.Sqrt(-2.0 * sys.Log(v12) / v12)
|
901 |
x = s + 0.5 * t
|
902 |
gds = x * x
|
903 |
If t >= 0.0 Then Return (gds / lambda) Immediate acceptance
|
904 |
|
905 |
u = nextDouble() Step 3. Uniform random number
|
906 |
If d * u <= t * t * t Then Return (gds / lambda) Squeeze acceptance
|
907 |
|
908 |
If a <> aaa Then Step 4. Set-up for hat case
|
909 |
aaa = a
|
910 |
r = 1.0 / a
|
911 |
q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r
|
912 |
If a > 3.686 Then
|
913 |
If a > 13.022 Then
|
914 |
b = 1.77
|
915 |
si = 0.75
|
916 |
c = 0.1515 / s
|
917 |
Else
|
918 |
b = 1.654 + 0.0076 * ss
|
919 |
si = 1.68 / s + 0.275
|
920 |
c = 0.062 / s + 0.024
|
921 |
End If
|
922 |
Else
|
923 |
b = 0.463 + s - 0.178 * ss
|
924 |
si = 1.235
|
925 |
c = 0.195 / s - 0.079 + 0.016 * s
|
926 |
End If
|
927 |
End If
|
928 |
If x > 0.0 Then Step 5. Calculation of q
|
929 |
v = t / (s + s) Step 6.
|
930 |
If sys.Abs(v) > 0.25 Then
|
931 |
q = q0 - s * t + 0.25 * t * t + (ss + ss) * sys.Log(1.0 + v)
|
932 |
Else
|
933 |
q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v
|
934 |
End If Step 7. Quotient acceptance
|
935 |
If sys.Log(1.0 - u) <= q Then Return (gds / lambda)
|
936 |
End If
|
937 |
|
938 |
Do Step 8. Double exponential deviate t
|
939 |
Do
|
940 |
e = -sys.Log(nextDouble())
|
941 |
u = nextDouble()
|
942 |
u = u + u - 1.0
|
943 |
sign_u = If(u > 0, 1.0, -1.0)
|
944 |
t = b + (e * si) * sign_u
|
945 |
Loop While t <= -0.71874483771719 Step 9. Rejection of t
|
946 |
v = t / (s + s) Step 10. New q(t)
|
947 |
If sys.Abs(v) > 0.25 Then
|
948 |
q = q0 - s * t + 0.25 * t * t + (ss + ss) * sys.Log(1.0 + v)
|
949 |
Else
|
950 |
q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v
|
951 |
End If
|
952 |
If q <= 0.0 Then Continue Do Step 11.
|
953 |
If q > 0.5 Then
|
954 |
w = sys.Exp(q) - 1.0
|
955 |
Else
|
956 |
w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * q + e1) * q
|
957 |
End If Step 12. Hat acceptance
|
958 |
If c * u * sign_u <= w * sys.Exp(e - 0.5 * t * t) Then
|
959 |
x = s + 0.5 * t
|
960 |
Return (x * x / lambda)
|
961 |
End If
|
962 |
Loop
|
963 |
End If
|
964 |
End Function
|
965 |
End Class
|
966 |
End Namespace
|