/* 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======================================================================== This C language version was written by Jim Butler, and was based on a FORTRAN program posted by David LaSalle of Florida State University. */ #include #include #define TRUE -1 #define FALSE 0 #define boolean int static void rmarin(int ij, int kl); void ranmar(float rvec[], int len); void main() { float temp[101]; int ij, kl, len, i; /* random seeds for the test case: */ ij = 1802; kl = 9373; /* do the initialization */ rmarin(ij,kl); /* generate 20,000 random numbers */ len = 100; for (i=1; i<=200; i++) ranmar(temp, len); /* C If the random number generator is working properly, the next six random C numbers should be: C 6533892.0 14220222.0 7275067.0 C 6172232.0 8354498.0 10633180.0 */ len = 6; ranmar(temp,len); for (i=1; i<=6; i++) printf("%12.1f ", 4096.0*4096.0*temp[i]); } static float u[98], c, cd, cm; static int i97, j97; static boolean test = FALSE; static void rmarin(ij,kl) int ij, kl; { /* 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 */ int i, j, k, l, ii, jj, m; float s, t; if (ij<0 || ij>31328 || kl<0 || kl>30081) { puts("The first random number seed must have a value between 0 and 31328."); puts("The second seed must have a value between 0 and 30081."); exit(1); } 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 += t; t *= 0.5; } u[ii] = s; } c = 362436.0 / 16777216.0; cd = 7654321.0 / 16777216.0; cm = 16777213.0 / 16777216.0; i97 = 97; j97 = 33; test = TRUE; } void ranmar(rvec, len) float rvec[]; /* len random numbers are placed in rvec[1..len] */ int 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. */ { int ivec; float uni; if (test==FALSE) { puts("Call the init routine rmarin() before calling ranmar()."); exit(2); } for (ivec=1; ivec<=len; ivec++) { uni = u[i97] - u[j97]; if (uni < 0.0) uni += 1.0; u[i97] = uni; i97--; if (i97==0) i97 = 97; j97--; if (j97==0) j97 = 97; c -= cd; if (c<0.0) c += cm; uni -= c; if (uni<0.0) uni += 1.0; rvec[ivec] = uni; } }