MINUIT%FSU.MFENET@CCC.MFECC.LLNL.GOV (06/30/89)
Mail message for mkahn@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@scri1.scri.fsu.edu (128.186.4.1) or minuit%fsu@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@DWIFH1.BITNET> and Helmut Weber <WEBER@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 *************************** */
blbates@AERO4.LARC.NASA.GOV ("Brent L. Bates TAD/TAB ms294 x42854") (06/30/89)
Thanks for the C version. -- Brent L. Bates NASA-Langley Research Center M.S. 294 Hampton, Virginia 23665-5225 (804) 864-2854 E-mail: blbates@aero4.larc.nasa.gov or blbates@aero2.larc.nasa.gov