random number generators
- To: mathgroup at smc.vnet.net
- Subject: [mg4758] random number generators
- From: wilson figueroa <wfigueroa at mosquito.com>
- Date: Thu, 5 Sep 1996 02:00:08 -0400
- Organization: Mosquito Net, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
Hello group,
I was wondering if anyone out there has done any of the following
implementations of random number generators in Mma?
I could really use some help here!
The first one is the linear congruential method
Question 1: Linear Congruential Generator (LCG)
The algorithm is:
x[n+1] = (7^7*x[n])Mod[2^31-1]
and uniform variates are given by:
u[n] = x[n+1]/(2^31-1).
I implemented this but it does not seem to work properly.
The implementation I used was:
$RecursionLimit = 100000;
x[0] = 1;
x[n+1] := Mod[(7^7*x[n]), 2^31 -1]//N
u[n] := x[n]/(2^31-1)
Table[x[n], {n,1,10000}]
Now when I do a Table of say 10000 of these, the results do not seem
correct. Is there a more efficient way of implementing this rather
straightforward rng that is fast and gives correct results?
Question 2: Marsaglia generator (subtract with borrow)
In an article "Some portable very long period rng's" Marsaglia and
Zaman (Computers in Physics, Jan/FEb 1994) algorithms are given
which allow combination generators to be used.
The first question I have is on the Subtract with borrow (which were
recently found to fail a statistical test).
The algorithm is
a[n] = (a[n-2] - a[n-5] - "C") mod (2^32-10)
the period of this generator is 2^160
uniform deviates are given by: u[n] = a[n]/(2^32-1)
The "C" here is a carry bit, not a constant!
This implementation has modulus m = 2^32 - 10, five current
integer values, say, v, w, x, y, z, and current c.
The generating procedure is:
If y < v + c then s = v + c + - y - 10 and c = 1
Else s = y - v - c and c = 0.
Then s is the output ans promotions provide a new
set of five: v = w, w = x, x = y, y = s.
The period is 2^160 for any initial values, not all zero.
Has anyone implemented such an algorithm using Mathematica?
I have several concerns including floating point versus integer
machine precision.
I would really appreciate any help for what appeared to me to be
a simple algorithm!
Question 3: Mix a LCG with a SWB
An algorithm is given in which a LCG and SWB are combined.
The algorithms are as follows:
LCG: x[n] = (69069*x[n-1]+1013904243)mod(2^32)
SWB: x[n] = (x[n-2]-x[n-3]-"C")mod(2^32-18)
The "C" is a carry bit again.
A "C-program" is given as follows:
typedef unsigned long int unlong;
unlong x=521288629, y=362436069, z=16163801, c=1, n=1131199209;
ulong mzran13()
{ long int s;
if (y>x+c) {s=y-(x+c); c=0; )
else {s=y-(x+c)-18; c=1; }
x=y; y=z; return (z=s) + (n-60069*n+1013904243);
};
void ran13set(ulong xx, ulong yy, ulong zz, ulong nn)
{ x = xx; y = yy; z = zz; c = y > z; }
My first problem is I can't read most of this code and secondly, I
would like to implement this algorithm in Mathematica.
Question 4:
Has anyone implemented the Inverse Pseudorandom Number Generator
discussed in an article by Peter Hellelalek (Univ of Salzburg)
using Mathematica?
The algorithm is:
y[n+1] = (a(not y[n]) + b)mod p
Note the not y[n] is (1/y[n]) provided y[n] does not equal zero,
if y[n] equals zero, set not y[n] to zero.
See http;//random.mat.sbg.ac.at. for more details.
I know there are many questions here and I apologize if this is bad
netiquette.
If there is anyone out there with more experience in these areas that
can help me, it would be greatly appreciated.
Thank you.
Wilson Figueroa
"Why is there so much to explore in a finite amount of time?"
wfiguero at mosquito.com
==== [MESSAGE SEPARATOR] ====