/* CVS info */ /* \$Date: 2005/07/28 15:31:14 \$ */ /* \$Revision: 1.3 \$ */ /* \$RCSfile: randc.c,v \$ */ /* \$Name: rel_5 \$ */ static char cvs_info[] = "BMARKGRP \$Date: 2005/07/28 15:31:14 \$ \$Revision: 1.3 \$ \$RCSfile: randc.c,v \$ \$Name: rel_5 \$"; #define min(A,B) ((A) < (B) ? (A) : (B)) #ifdef TEST /* TEST PROGRAM TO CALL IRAND */ /******************************/ #define ISEED -1431 main() { unsigned int i,j,ib,k,irand(); long ia[5000]; void vrand(); i = irand(ISEED); ib = irand(0); printf("0 %o\n", i); printf("1 %o\n", ib); for (j=0; j<20; j++) { for (k=0; k< 100; k++) vrand(0, ia, 5000); printf("%d %o\n", j+2, ia[4999]); } } #endif /* This random number generator is as described on page 26 (Algorithm A) of "The Art of Computer Programming, Vol 2", by D. E. Knuth. This algorithm uses a linear congruential generator to fill a 55 long register with 32-bit numbers. The generator steps according to the rule X_N = X_N-24 + X_N-55 mod 2^32. This has a cycle length of (2^55 - 1)*(2^31) or about 2^86, when at least one of the initial 55 fill numbers is odd. Since the largest run of even numbers in the initial fill generator is 28, this condition is guaranteed. Random numbers are then successive entries in the continued sequence. ISEED < 0 initialize the generator to iseed IRAND the returned random deviate between 0 and 2^32 Subroutine VRAND(SEED,ARRAY,LENGTH) is a companion routine which does the same function, except returns a vector of results of length N; regardless of which is called, the same results will be returned. April 19, 1993 */ #define JQ 127773 #define JA 16807 #define JR 2836 #define JM 0x7fffffff #define I32M1 0xffffffff #define JSEED0 967560787 #define ILEN 1024 unsigned long ix[ILEN]; /* The random register */ int iflag = 0; int ipoint; unsigned int irand(int iseed) { int jseed, jh, jl; extern unsigned long ix[ILEN]; /* The random register */ extern int iflag; extern int ipoint; int i; /* printf("iflag = %d, ipoint = %d\n",iflag, ipoint); */ if ((iseed < 0) || !iflag ) { iflag = -1; jseed = -iseed; if ( jseed >= JM || jseed <= 0) jseed = JSEED0; for (i = 0; i < 55; i++) { jh = jseed/JQ; jl = jseed - (jh*JQ); jseed = JA*jl - JR*jh; if(jseed < 0) jseed = jseed + JM; ix[i] = jseed; } for (i = 0; i < 55; i++) { jh = jseed/JQ; jl = jseed - (jh*JQ); jseed = JA*jl - JR*jh; if (jseed < 0) jseed = jseed + JM; ix[i] = jseed ^ (ix[i] << 1); } /* * Run up generator */ for (i = 55; i ILEN) /* Have we used all the available numbers? */ { for (i= 0; i<24; i++) { ix[i] = I32M1 & (ix[i + ILEN -24] + ix[i+ILEN - 55]); } for (i = 24; i<55; i++) { ix[i] = I32M1 & (ix[i-24] + ix[i+ILEN - 55]); } for (i = 55; i= JM || jseed <= 0) jseed = JSEED0; for (i = 0; i < 55; i++) { jh = jseed/JQ; jl = jseed - (jh*JQ); jseed = JA*jl - JR*jh; if(jseed < 0) jseed = jseed + JM; ix[i] = jseed; } for (i = 0; i < 55; i++) { jh = jseed/JQ; jl = jseed - (jh*JQ); jseed = JA*jl - JR*jh; if (jseed < 0) jseed = jseed + JM; ix[i] = jseed ^ (ix[i] << 1); } /* * Run up generator */ for (i = 55; i ILEN) { ranrun(); ipoint = 1; } /* * How many can we move? */ jj = min(ILEN+1-ipoint, n - ii); #ifdef PRAGMA_YES # pragma _CRI ivdep #endif for(i=0; i