Wichmann-Hill Pseudo Random Number Generator: an alternative for VB Rnd() function

[anhn]

The native VB Rnd() function can generate a sequence of upto only 16,777,216 numbers before they are repeated. This is adequate in most cases but some projects may require many more numbers.

Below is what I wrote in Excel (can be used in VB except the TestSub) based on the algorithm developed by B.A. Wichman and I.D. Hill that is implemented in Excel 2003 for WorksheetFunction RAND().

As said, this algorithm guarantees to provide a sequence of more than 10,000,000,000,000 (10^13) numbers before they are repeated.

I am happy with everyone's input to make this function better and faster.

Please point out any problem you see.

On my PC, withExcel 2003, it can generate 1 million numbers in every around 550 milliseconds (including time for storing values.)

Put the top 3 Code blocks in one module and the test Sub in another module.
Have fun!

'====================================================================================
'http://support.microsoft.com/kb/828795
'
'The basic idea is that if you take three random numbers on [0,1] and sum them,
'the fractional part of the sum is itself a random number on [0,1].
'The critical statements in the Fortran code listing from the original
'Wichman and Hill article are:
'
'C  IX, IY, IZ SHOULD BE SET TO INTEGER VALUES BETWEEN 1 AND 30000 BEFORE FIRST ENTRY
'IX = MOD(171 * IX, 30269)
'IY = MOD(172 * IY, 30307)
'IZ = MOD(170 * IZ, 30323)
'RANDOM = AMOD(FLOAT(IX) / 30269.0 + FLOAT(IY) / 30307.0 + FLOAT(IZ) / 30323.0, 1.0)
'====================================================================================
Option Explicit
Private ix As Long, iy As Long, iz As Long

Sub RandomizeX(Optional ByVal Number)
   Const MaxLong As Double = 2 ^ 31 - 1
   Dim n As Long
   Dim d As Double
   
   If IsMissing(Number) Then
      n = Timer * 60
      '-- Timer is only updated in every 1/60th second.
      '-- Multiply by 60 to reduce the chance of seed
      '-- to be repeated in subsequence calls of RandomizedX
   Else
      d = Abs(Int(Val(Number)))
      If d > MaxLong Then '-- prevent Long overflow
         d = d - Int(d / MaxLong) * MaxLong
      End If
      n = d
   End If
   ix = (n Mod 30269)
   iy = (n Mod 30307)
   iz = (n Mod 30323)
   '-- ix, iy, iz cannot be 0
   If ix = 0 Then ix = 171
   If iy = 0 Then ix = 172
   If iz = 0 Then ix = 170
End Sub

Function RndX(Optional ByVal Number As Long = 1) As Double
   Dim r As Double
   
   If ix = 0 Then '-- ix, iy, iz cannot be 0.
      ix = 171    '-- Initial values of ix, iy and iz are 0.
      iy = 172    '-- If any of these equal 0,
      iz = 170    '-- it will be stucked with 0 forever.
   End If
   If Number <> 0 Then
      If Number < 0 Then
         RandomizeX Number
      End If
      ix = (171 * ix) Mod 30269 '-- This has been tested:
      iy = (172 * iy) Mod 30307 '-- ix, iy, iz will never be 0.
      iz = (170 * iz) Mod 30323 '-- ---------------------------
   End If
   r = ix / 30269# + iy / 30307# + iz / 30323#
   RndX = r - Int(r)
End Function

Run TestRndX(). If you want to stop, use Ctrl-Break then un-remark 'Exit Do
then let the Sub continues to exit properly to reset ScreenUpdating.

Option Explicit
Private Declare Function GetTickCount Lib "kernel32" () As Long

Sub TestRndX()
   Dim r0 As Double
   Dim r As Double
   Dim rMin As Double
   Dim rMax As Double
   Dim rSum As Double
   Dim rAvg As Double
   Dim TenPctCount(0 To 9) As Double
   Dim p As Integer
   Dim n As Long
   Dim m As Long
   Dim t As Double
   Dim aCell As Range
   Dim CellRow As Long
   Const BottomRow = 65536
   Const Million = 1000000
   Dim Headers As Variant
   Dim ColWidths As Variant
   
   '----- Prepare Sheet for report ---------------------------------------
   Headers = Array("Count", "Cur Rnd", "millisecs", "Min", "Max", "Avg", _
                   "<0.1", "<0.2", "<0.3", "<0.4", "<0.5", _
                   "<0.6", "<0.7", "<0.8", "<0.9", "<1.0")
   ColWidths = Array(15, 18, 8, 18, 18, 18, _
                     13, 13, 13, 13, 13, 13, 13, 13, 13, 13)
   With Sheet1
      For p = 0 To UBound(Headers)
         .Cells(1, p + 1).ColumnWidth = ColWidths(p)
      Next
      .Range(.Cells(1, 1), .Cells(1, 16)) = Headers
      .Range("A:A,C:C,G:P").NumberFormat = "#,##0"
      .Range("B:B,D:F").NumberFormat = "0.000000000000000"
   '-----------------------------------------------------------------------
      RandomizeX
      r0 = RndX '---- FIRST Random Number
      rMin = r0
      rMax = r0
      rSum = r0
      rAvg = r0
      CellRow = 2
      .Cells(CellRow, 1) = 0
      .Cells(CellRow, 2) = r0
      Application.ScreenUpdating = False
      Do While CellRow < BottomRow
         '=========================
         'Exit Do '-- Un-remark this after Control-Break then let
                  '   the Sub continue to reset ScreenUpdating
         '=========================
         t = GetTickCount
         For n = 1 To Million
            r = RndX
            If rMin > r Then rMin = r
            If rMax < r Then rMax = r
            rSum = rSum + r
            p = Int(r * 10)
            TenPctCount(p) = TenPctCount(p) + 1
            If r = r0 Then
               .Cells(CellRow, 1) = m * Million + n
               .Cells(CellRow, 2) = r
               .Cells(CellRow, 3) = "Repeated!"
               Exit Do
            End If
         Next
         t = GetTickCount - t
         m = m + 1
         CellRow = CellRow + 1
         .Cells(CellRow, 1) = CDbl(m) * Million
         .Cells(CellRow, 2) = r
         .Cells(CellRow, 3) = t
         .Cells(CellRow, 4) = rMin
         .Cells(CellRow, 5) = rMax
         .Cells(CellRow, 6) = rSum / (CDbl(m) * Million) '-- rAvg
         For p = 0 To 9
            .Cells(CellRow, 7 + p) = TenPctCount(p)
         Next
         If m Mod 10 = 0 Then
            Application.ScreenUpdating = True
            ActiveWindow.ScrollRow = CellRow - 10
            DoEvents
            Application.ScreenUpdating = False
         End If
      Loop
   End With
   Application.ScreenUpdating = True
End Sub