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