2019-05-21

### Navigation ### Others

Syntax highlighting by Prism

## MPRNG

### Algorithm creator(s)

George Marsaglia, Arif Zaman, later modified by F. James

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)

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
``````

Mirror provided by Knuth Konrad