1 #Region "Microsoft.VisualBasic::120bc0a3db9510757289f89883c14448, 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 Overloads) Sub 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 ' ' * MersenneTwisterFast.java
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     ''' <summary>
80     ''' MersenneTwisterFast:
81     '''  
82     ''' A simulation quality fast random number generator (MT19937) with the same
83     ''' public methods as java.util.Random.
84     '''  
85     '''  
86     ''' About the Mersenne Twister. This is a Java version of the C-program for
87     ''' MT19937: Integer version. next(32) generates one pseudorandom unsigned
88     ''' integer (32bit) which is uniformly distributed among 0 to 2^32-1 for each
89     ''' call. next(int bits) >>>'s by (32-bits) to get a value ranging between 0 and
90     ''' 2^bits-1 long inclusive; hope that's correct. setSeed(seed) set initial
91     ''' values to the working area of 624 words. For setSeed(seed), seed is any
92     ''' 32-bit integer except for 0.
93     '''  
94     ''' Reference. M. Matsumoto and T. Nishimura, "Mersenne Twister: A
95     ''' 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator",
96     ''' <i>ACM Transactions on Modeling and Computer Simulation,</i> Vol. 8, No. 1,
97     ''' January 1998, pp 3--30.
98     '''  
99     '''  
100     ''' Bug Fixes. This implementation implements the bug fixes made in Java 1.2's
101     ''' version of Random, which means it can be used with earlier versions of Java.
102     ''' See <a href=
103     ''' "http://www.javasoft.com/products/jdk/1.2/docs/api/java/util/Random.html">
104     ''' the JDK 1.2 java.util.Random documentation</a> for further documentation on
105     ''' the random-number generation contracts made. Additionally, there's an
106     ''' undocumented bug in the JDK java.util.Random.nextBytes() method, which this
107     ''' code fixes.
108     '''  
109     '''  
110     ''' Important Note. Just like java.util.Random, this generator accepts a long
111     ''' seed but doesn't use all of it. java.util.Random uses 48 bits. The Mersenne
112     ''' Twister instead uses 32 bits (int size). So it's best if your seed does not
113     ''' exceed the int range.
114     '''  
115     '''  
116     ''' <a href="http://www.cs.umd.edu/users/seanl/">Sean Luke's web page</a>
117     '''  
118     '''  
119     ''' - added shuffling method (Alexei Drummond)
120     '''  
121     ''' - added gamma RV method (Marc Suchard)
122     '''  
123     ''' This is now package private - it should be accessed using the instance in
124     ''' Random
125     ''' </summary>
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         ''' <summary>
172         ''' Constructor using the time of day as default seed.
173         ''' </summary>
174         Public Sub New()
175             Me.New(Now.ToBinary + seedAdditive_)
176             seedAdditive_ += nextInt()
177         End Sub
178
179         ''' <summary>
180         ''' Constructor using a given seed. Though you pass this seed in as a long,
181         ''' it's best to make sure it's actually an integer.
182         ''' </summary>
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         ''' <summary>
194         ''' Initalize the pseudo random number generator. The Mersenne Twister only
195         ''' uses an integer for its seed; It's best that you don't pass in a long
196         ''' that's bigger than an int.
197         ''' </summary>
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         ''' <summary>
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         ''' </summary>
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         ''' <summary>
766         ''' Returns a uniform random permutation of int objects in array
767         ''' </summary>
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         ''' <summary>
779         ''' Shuffles an array.
780         ''' </summary>
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         ''' <summary>
792         ''' Shuffles an array. Shuffles numberOfShuffles times
793         ''' </summary>
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         ''' <summary>
808         ''' Returns an array of shuffled indices of length l.
809         ''' </summary>
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         ''' <summary>
826         ''' Returns a uniform random permutation of ints 0,...,l-1
827         ''' </summary>
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         ''' <summary>
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 &lt;
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         ''' </summary>
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