MPRNG
Algorithm creator(s)
George Marsaglia, Arif Zaman, later modified by F. James
PB author(s)
Torsten Rienow
Description
'The best known random number generator available', it passes all PRNG tests and has an enormous period of 2^144. It is a combination of Fibonacci sequence (using subtraction plus one, modulo one) and arithmetic sequence (using subtraction).
Note
See Paul Bourke's analysis of the algorithm. (Mirror note: that site's no longer available)
Source
n/a
See also
n/a
Source Code
Download source code file mprng.bas (Right-click -> "Save as ...")
'The code below covers some aspects from the thread found at http://www.powerbasic.com/support/forums/Forum6/HTML/003084.html
'The code below is a port from C to PB from the following source http://astronomy.swin.edu.au/~pbourke/analysis/random/
'My code is Public Domain
'(c) by Torsten Rienow
'The code below is Public Domain
#Compile Exe
Type g_Type 'covers the global/module variables from c implementation
u(96) As Double
c As Double
CD As Double
cm As Double
i97 As Dword
j97 As Dword
test As Dword
End Type
Sub RandomInitialise(ByVal ij As Dword, ByVal kl As Dword, g_Data As g_Type)
Dim s As Double
Dim t As Double
Dim ii As Dword
Dim i As Dword
Dim j As Dword
Dim k As Dword
Dim l As Dword
Dim jj As Dword
Dim m As Dword
'/*
'Handle the seed range errors
'First random number seed must be between 0 and 31328
'Second seed must have a value between 0 and 30081
'*/
If (ij < 0 Or ij > 31328 Or kl < 0 Or kl > 30081) Then
ij = 1802
kl = 9373
End If
i = (ij / 177) Mod 177 + 2
j = (ij Mod 177) + 2
k = (kl / 169) Mod 178 + 1
l = (kl Mod 169)
ii = 0
While ii<97
s = 0.0
t = 0.5
jj = 0
While jj<24
m = (((i * j) Mod 179) * k) Mod 179
i = j
j = k
k = m
l = (53 * l + 1) Mod 169
If (((l * m Mod 64)) >= 32) Then
s = s + t
End If
t = t * 0.5
Incr jj
Wend
g_data.u(ii) = s
Incr ii
Wend
g_data.c = 362436.0 / 16777216.0
g_data.cd = 7654321.0 / 16777216.0
g_data.cm = 16777213.0 / 16777216.0
g_data.i97 = 97
g_data.j97 = 33
g_data.test = 1
End Sub
'/*
' This is the random number generator proposed by George Marsaglia in
' Florida State University Report: FSU-SCRI-87-50
'*/
Function RandomUniform(g_Data As g_Type) As Double
Dim uni As Double
'/* Make sure the initialisation routine has been called */
If (g_data.test = 0) Then
RandomInitialise 1802,9373, g_Data
End If
uni = g_Data.u(g_Data.i97 - 1) - g_Data.u(g_Data.j97 - 1)
If (uni <= 0.0) Then
Incr uni
End If
g_Data.u(g_Data.i97 - 1) = uni
Decr g_Data.i97
If (g_Data.i97 = 0) Then
g_Data.i97 = 97
End If
Decr g_Data.j97
If (g_Data.j97 = 0) Then
g_Data.j97 = 97
End If
g_Data.c = g_Data.c - g_data.cd
If (g_Data.c < 0.0) Then
g_Data.c = g_Data.c + g_Data.cm
End If
uni = uni - g_Data.c
If (uni < 0.0) Then
Incr uni
End If
Function = uni
End Function
'/*
' ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
' THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
' VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
' The function returns a normally distributed pseudo-random number
' with a given mean and standard devaiation. Calls are made to a
' function subprogram which must return independent random
' numbers uniform in the interval (0,1).
' The algorithm uses the ratio of uniforms method of A.J. Kinderman
' and J.F. Monahan augmented with quadratic bounding curves.
'*/
Function RandomGaussian(ByVal mean As Double, ByVal stddev As Double) As Double
Dim q As Double
Dim u As Double
Dim v As Double
Dim x As Double
Dim y As Double
Dim g_data As g_Type
'/*
'Generate P = (u,v) uniform in rect. enclosing acceptance region
'Make sure that any random numbers <= 0 are rejected, since
'gaussian() requires uniforms > 0, but RandomUniform() delivers >= 0.
'*/
Do
u = RandomUniform(g_Data)
v = RandomUniform(g_Data)
If (u <= 0.0 Or v <= 0.0) Then
u = 1.0
v = 1.0
End If
v = 1.7156 * (v - 0.5)
'/* Evaluate the quadratic form */
x = u - 0.449871
y = Abs(v) + 0.386595
q = x * x + y * (0.19600 * y - 0.25472 * x)
'/* Accept P if inside inner ellipse */
If q < 0.27597 Then
Exit Do
End If
'/* Reject P if outside outer ellipse, or outside acceptance region */
Loop Until 0 = ((q > 0.27846) Or (v * v > -4.0 * Log(u) * u * u))
'/* Return ratio of P's coordinates as the normal deviate */
Function = (mean + stddev * v / u)
End Function
'/*
' Return random integer within a range, lower -> upper INCLUSIVE
'*/
Function RandomInt(ByVal lower As Dword, ByVal upper As Dword) As Dword
Dim g_Data As g_Type
Function = (RandomUniform(g_Data) * (upper - lower + 1)) + lower
End Function
'/*
' Return random float within a range, lower -> upper
'*/
Function RandomDouble(ByVal lower As Double, ByVal upper As Double) As Double
Dim g_Data As g_Type
Function = ((upper - lower) * RandomUniform(g_Data) + lower)
End Function
Function GeneralTest() As String
Dim g_Data As g_Type
Dim s As String
Dim i As Dword
Dim r As Double
'/* Generate 20000 random numbers */
RandomInitialise 1802, 9373, g_Data
i = 0
While i < 20000
r = RandomUniform(g_Data)
Incr i
Wend
'/*
'If the random number generator is working properly,
'the next six random numbers should be:
'6533892.0 14220222.0 7275067.0
'6172232.0 8354498.0 10633180.0
'*/
i = 0
While i < 6
s = s + $Tab + Format$(4096.0 * 4096.0 * RandomUniform(g_Data), "##########.0") + $CrLf
Incr i
Wend
Function = s
End Function
Function PbMain
Dim s As String
s = "General Test" + $CrLf + GeneralTest + $CrLf + _
"Random Int" + $CrLf + $Tab + Format$(RandomInt(0, 100)) + $CrLf + $CrLf + _
"Random Double" + $CrLf + $Tab + Format$(RandomDouble(0, 100))
MsgBox s
End Function