C version of Marsaglia's ran()
MINUIT%FSU.MFENET at CCC.MFECC.LLNL.GOV
MINUIT%FSU.MFENET at CCC.MFECC.LLNL.GOV
Fri Jun 30 03:21:41 AEST 1989
Mail message for mkahn at tripost.com (my mailer couldn't find him)
-----------------------------------------------------------------------
Here is the C version of Marsaglia's random number generator. If you'll
send me your US mail address, I will send you some papers on the generator
and the tests that were used to verify it.
- David LaSalle
SCRI
Florida State University
Tallahassee, FL 32306-4052
(904)644-8532
arpanet: minuit at scri1.scri.fsu.edu (128.186.4.1)
or
minuit%fsu at nmfecc.arpa
---------------------------- cut here -------------------------------------
/*
This Random Number Generator is based on the algorithm in a FORTRAN
version published by George Marsaglia and Arif Zaman, Florida State
University; ref.: see original comments below.
At the fhw (Fachhochschule Wiesbaden, W.Germany), Dept. of Computer
Science, we have written sources in further languages (C, Modula-2
Turbo-Pascal(3.0, 5.0), Basic and Ada) to get exactly the same test
results compared with the original FORTRAN version.
April 1989
Karl-L. Noell <NOELL at DWIFH1.BITNET>
and Helmut Weber <WEBER at DWIFH1.BITNET>
*/
#define FALSE 0
#define TRUE 1
/*
C This random number generator originally appeared in "Toward a Universal
C Random Number Generator" by George Marsaglia and Arif Zaman.
C Florida State University Report: FSU-SCRI-87-50 (1987)
C
C It was later modified by F. James and published in "A Review of Pseudo-
C random Number Generators"
C
C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
C (However, a newly discovered technique can yield
C a period of 10^600. But that is still in the development stage.)
C
C It passes ALL of the tests for random number generators and has a period
C of 2^144, is completely portable (gives bit identical results on all
C machines with at least 24-bit mantissas in the floating point
C representation).
C
C The algorithm is a combination of a Fibonacci sequence (with lags of 97
C and 33, and operation "subtraction plus one, modulo one") and an
C "arithmetic sequence" (using subtraction).
C
C========================================================================
*/
main()
{
float temp[100];
int ij, kl, len, i;
void rmarin(),ranmar();
/* These are the seeds needed to produce the test case results */
ij = 1802;
kl = 9373;
/* Do the initialization */
rmarin(ij,kl);
/* Generate 20000 random numbers */
len = 100;
for ( i = 1; i<=200; i++)
ranmar (temp, len);
/* 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
*/
len = 6;
ranmar(temp, len);
for (i=1; i<=6; i++)
printf ("%12.1f\n", 4096.0*4096.0*temp[i-1]);
}
/* ***************************** End of main ************************** */
/*
C This is the initialization routine for the random number generator RANMAR()
C NOTE: The seed variables can have values between: 0 <= IJ <= 31328
C 0 <= KL <= 30081
C The random number sequences created by these two seeds are of sufficient
C length to complete an entire calculation with. For example, if sveral
C different groups are working on different parts of the same calculation,
C each group could be assigned its own IJ seed. This would leave each group
C with 30000 choices for the second seed. That is to say, this random
C number generator can create 900 million different subsequences -- with
C each subsequence having a length of approximately 10^30.
C
C Use IJ = 1802 & KL = 9373 to test the random number generator. The
C subroutine RANMAR should be used to generate 20000 random numbers.
C Then display the next six random numbers generated multiplied by 4096*4096
C If the random number generator is working properly, the random numbers
C should be:
C 6533892.0 14220222.0 7275067.0
C 6172232.0 8354498.0 10633180.0
*/
/* Globals (accessible from rmarin() and ranmar() : */
float u[97],c,cd,cm;
int i97,j97;
int test=FALSE;
void rmarin(ij,kl)
int ij,kl;
{
float s, t ;
int ii, i, j, k, l, jj, m ;
if((ij < 0) || (ij > 31328) || (kl < 0) || (kl > 30081))
printf (" The first random number seed must have a value between ",
" 0 and 31328\n",
" The second seed must have a value between 0 and 30081\n");
i = (ij/177)%177 + 2 ;
j = (ij%177) + 2 ;
k = (kl/169)%178 + 1 ;
l = (kl%169) ;
for ( ii = 1; ii <= 97; ii++) {
s = 0.0 ;
t = 0.5 ;
for ( jj = 1; jj <= 24; jj++) {
m = (((i*j)%179)*k)%179 ;
i = j ;
j = k ;
k = m ;
l = (53*l+1)%169 ;
if (((l*m%64)) >= 32) s = s + t ;
t = 0.5 * t ;
}
u[ii-1] = s ;
}
c = 362436.0 / 16777216.0 ;
cd = 7654321.0 / 16777216.0 ;
cm = 16777213.0 /16777216.0 ;
i97 = 97 ;
j97 = 33 ;
test = TRUE ;
}
/* *************************** End of rmarin *************************** */
/* subroutine ranmar(RVEC, LEN)
C This is the random number generator proposed by George Marsaglia in
C Florida State University Report: FSU-SCRI-87-50
C It was slightly modified by F. James to produce an array of pseudorandom
C numbers.
*/
void ranmar (rvec,len)
float rvec [];
int len ;
{
/* common /raset1/ U, C, CD, CM, I97, J97, TEST */
int ivec;
float uni;
if(!test) {
printf (" Call the init routine (RMARIN) before calling ranmar\n");
/* here you may insert an appropriate break, exit, abort etc. */
}
for (ivec = 1; ivec <= len; ivec++) {
uni = u[i97-1] - u[j97-1];
if( uni <= 0.0 ) uni = uni + 1.0 ;
u[i97-1] = uni ;
i97 = i97 - 1 ;
if(i97 == 0) i97 = 97 ;
j97 = j97 - 1 ;
if(j97 == 0) j97 = 97 ;
c = c - cd ;
if( c < 0.0 ) c = c + cm ;
uni = uni - c ;
if( uni < 0.0 ) uni = uni + 1.0 ;
rvec[ivec-1] = uni ;
}
}
/* *************************** End of ranmar *************************** */
More information about the Comp.sys.sgi
mailing list