December 2003 - cypherpunks-legacy

(fwd) Random Numbers - Request for feedback
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!agate!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd(a)netcom.com Subject: Random Numbers - Request for feedback Message-ID: <deleydCsIB28.KEI(a)netcom.com> Organization: NETCOM On-line Communication Services (408 261-4700 guest) Date: Wed, 6 Jul 1994 06:51:43 GMT Lines: 43 Xref: bga.com sci.stat.math:1315 sci.math:15353 sci.math.num-analysis:3354 RE: Computer Generated Random Numbers A few closing comments and requests for further information: 1. All my tests on random number generators were performed on VAX/VMS computers. VAX uses a 32-bit architecture, so the random number generators I tested were ones which used a word size of 32 bits or less. I would be interested in anybody's test results of a random number generator utilizing a larger word size, such as xrand() using SIZE=63. 2. Anyone know of some good references on primitive polynomials mod 2 and their applications? They're used in additive congruential random number generators like the xrand() one tested here. They're also used by file transfer programs such as xmodem to insure error free transmission, and they're used in cryptography too. Anyone know of a good book on Abstract Algebra? (The ones I have just briefly touch the topic and then move on.) 3. Resolution: Usually the random number generator is set up to return a floating point value between 0 and 1. A typical floating point variable R can only represent a finite number of different values between 0 and 1. If you magnify the result too much the discreetness of the floating point datum will become obvious. For example, in VAX architecture the F-floating datum has a precision of approximately one part in 2**23. Multiplying R by a very large number N to create a random variable between 0 and N will fail if N is too large because some of the values between 0 and N have no corresponding R value which maps to them (i.e. the mapping is no longer a surjection or onto map). For an F_floating datum, N above 2**23 is obviously too large. But even below 2**23 there's still a problem of some bins having 2 R values which map to them while other bins have only 1. We need to get N small enough so that the number of R values which maps to any bin is about the same, close enough so that differences aren't noticed when we test the random number generator. -David Deley deleyd(a)netcom.com

1 0

(fwd) Random Numbers - Results of testing BSD random()
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd(a)netcom.com Subject: Random Numbers - Results of testing BSD random() Message-ID: <deleydCsIAyv.K6n(a)netcom.com> Organization: NETCOM On-line Communication Services (408 261-4700 guest) Date: Wed, 6 Jul 1994 06:49:42 GMT Lines: 119 Xref: bga.com sci.stat.math:1314 sci.math:15352 sci.math.num-analysis:3353 BSD random() Here are the partial results. Further tests were not performed due to lack of time. So far the generator appears to be comparable to a shuffled linear congruential generator. DEFINITION: Generating polynomial: x^31 + x^3 + 1 (primitive polynomial) Initialize circular queue of 31 elements using ANSI C linear congruential generator. Recursion formula: a[i] = a[i] + a[i-3] RATING: 1-D FAILS above 800,000 bpd (bins per dimension) 2-D FAILS above 3000 bpd 3-D FAILS above 210 bpd 4-D PASSES at 50 bpd (highest tested so far) 5-D not tested 6-D not tested 7-D not tested 8-D not tested This is an additive congruential type random number generator. An array table[31] is initially filled with random numbers using the ANSI C linear congruential random number generator. Random numbers are then generated using the recursion formula: table[k] = (table[k-31] + table[k-3]) mod 32 (Note that x**31 + x**3 + 1 is a primitive polynomial mod 2, which is being used here as a generator.) Since we are using the array table[] as a circular queue with 31 elements then table[k-31] is just table[k] before it gets replaced with the new value. The recursion formula becomes: table[k] = table[k] + table[k-3] The generator works well in practice. Knuth claims the sequence will have period 2**31 - 1. Knuth also claims there is very little theory to prove that this generator does or does not have desirable random properties. I would be interested if anyone knows of any recent developments in this area. -David Deley deleyd(a)netcom.com (So sorry, I lost the name of the original person who posted this code below which was used in the tests. -D.D.) /*** Code to implement random() & srandom() of BSD Unix. It was taken (though coded somewhat differently) from the Gnu BSD implementation. ***/ #include <stdio.h> #include <stdlib.h> #define LONG31 #ifdef LONG31 /* x^31 + x^3 + 1 */ #define SIZE 31 #define SIZE1 30 #define P1 3 #define P2 0 #else /* LONG63: x^63 + x + 1 */ #define SIZE 63 #define SIZE1 62 #define P1 1 #define P2 0 #endif #define LONG_MAX 0x7fffffff int p1=P1, p2=P2; long table[SIZE]; /*** return a "random" number in range [0, LONG_MAX] */ long xrand () { int r; table[p1] = table[p1] + table[p2]; /* add two table elements */ r = (table[p1] >> 1) & LONG_MAX; /* throw least significant bit away */ if (p1 == SIZE1) { /* increment the table indexes */ p1 = 0; p2 = p2 + 1; } else if (p2 == SIZE1) { p1 = p1 + 1; p2 = 0; } else { p1 = p1 + 1; p2 = p2 + 1; } return (r); } /*** use a linear congruential type generator to seed the state table & cycle the entire table 10 times */ void sxrand (seed) long seed; { int i; table[0] = seed; for (i=1; i<SIZE; ++i) table[i] = (table[i-1] * 1103515145) + 12345; /* lousy */ for (i=0; i<10*SIZE; ++i) (void) xrand(); }

1 0

(fwd) Random Numbers - CHIKSN.FOR
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd(a)netcom.com Subject: Random Numbers - CHIKSN.FOR Message-ID: <deleydCsIBAp.KuF(a)netcom.com> Organization: NETCOM On-line Communication Services (408 261-4700 guest) Date: Wed, 6 Jul 1994 06:56:49 GMT Lines: 526 Xref: bga.com sci.stat.math:1316 sci.math:15354 sci.math.num-analysis:3355 {Approx. 520 lines} C CHIKSN.FOR C C This is the program which impliments the chi-square test used to test C random number generators. Presented here if you would wish to play C with it yourself, maybe do some testing of your own. See the paper C "Computer Generated Random Numbers" sections 4 and 6 for an explanation C of what this program does. C C This is not polished code you put on a shelf and admire, this is code C you dig your hands into and work with to make it do what you want it to. C C The main routine is the one to tinker with. This code is meant to be C modified to suit your needs. The goal is to fill up the bins with C balls. Here's a brief outline: C C 1. ASK USER INPUT: C a. Number of dimensions? NDIM C b. Number of bins per dimension? NBINSPD C c. Total number of balls to throw at bins? NBALLS C d. Number of tests to run? NCHITESTS C e. Random number generator to use (if more than one defined) C f. SEED value to initialize generator with C C 2. CREATE ARRAY BINS and ZERO ARRAY C C 3. THROW THE BALLS AT THE BINS and CALCULATE PROBABILITY: C C LOOP (do CHITEST=1 to NCHITESTS) C zero BIN array C LOOP (do J=1 to NBALLS) C get random numbers r1,r2,...,rn) C increment BIN(r1,r2,...,rn) C ENDLOOP C CALL CHSONE to calculate chi-square probability C ENDLOOP C CALL KSONE to calculate Kolmogorov-Smirnov probability C C The main routine here is where you put a call to your random number C generator, or for speed you can attempt a direct implimentation of your C random number generator to save the overhead of a call. (It can make a C difference when you call the random number generator 100 million times C in one test). C C The main routine defines a large one-dimensional array called BINS, the C maximum size of which would depend on your account quotas and machine C specific limitations. The array BINS keeps track of how many balls have C fallen into each bin. The size of array BINS determines the maximum C number of bins a user may select for a test. (10,000,000 bins is a C typical number you may want to use if possible.) C C So steps are: C C 1. Check definition of one-dimensional array NBINS C that it's not too large for your account quotas C or system limitations. C C 2. Place your random number generator to be tested where it C bluntly says "PLACE YOUR RANDOM NUMBER GENERATOR HERE". C The output is an integer IRANDOM between 0 and NBINS-1 C (NBINS is the number of bins chosen by the user). C C Currently the program is set up to use subroutine RAN1, a portable C random number generator from the book "NUMERICAL RECIPES: The Art C of Scientific Computing". I've had trouble on our UNIX system C not making array R(97) static even though the code says to. C Compiling with the -static qualifier works. C C The random number generator being tested is used to "randomly" select a C bin for the ball to fall in, and the counter for that bin is C incremented. Note for a multi-dimensional test we calculate the C appropriate index into the linear array BINS by hand. After all the C balls are thrown we call the subroutines to do the heavy math. C C All the subroutines should be fairly standard FORTRAN-77 modified C versions of routines from the book "NUMERICAL RECIPES: The Art of C Scientific Computing" by William H. Press, Brian P. Flannery, Saul A. C Teukolsky, and William T. Vetterling, and you should look there for C further reference as to what the routines are doing. (Note: the book C comes in several programming language forms including C, PASCAL, BASIC, C as well as FORTRAN, so you can take your pick and rewrite this code in C any language you please.) C C Note: C CHIKSN.FOR is currently set up to be run by a process with a very C large page file quota (pgflquo). If you get a 'exceed quota' error C attempting to run this then all you need to do is change the line C which reads: C C INTEGER*2 BINS(20 000 000) !The bins. C C to something smaller like: C C INTEGER*2 BINS(1 000 000) !The bins. C C C To compile: C $ FORTRAN CHIKSN C $ LINK CHIKSN C or C % f77 chiksn.f !(Some UNIX F77 compilers require -save option) C C Sample run: C Test MTH$RANDOM in 3-D with 10 bins per dimension and 10 balls per bin: C C $ RUN CHIKSN C Input number of dimensions NDIM: 3 C Input number of bins per dimension NBINSPD: 10 C Total number of bins = NBINSPD**NDIM = 1000 C Minimum number of balls = 5*NBINS = 5000 C Input total number of balls NBALLS: 10000 C Input number of Chi-Square tests NCHITESTS (min=2) : 2 C Choose random number generator to test C (1) MTH$RANDOM, C (2) RANDU, C (3) ANSI C, C (4) Microsoft C C (5) Turbo Pascal C (6) DES C : 1 C Input starting SEED value: 1 C C BALLS= 10000 CHISQ= 993.0002441 PROB= 0.4524292 C BALLS= 10000 CHISQ= 974.0001831 PROB= 0.2915459 C KS D= 0.5475708 PROB= 0.5863269 C C----------------------------------------------------------------------- PROGRAM CHIKSN C Perform a CHI-SQUARE test on a sequence of sets of N random numbers C NDIM = number of dimensions C NBINSPD = number of bins per dimension C NBINS = total number of bins. NBINS = NBINSPD**NDIM C NBALLS = total number of balls. Should be at least 5*(NBINS**NDIM) C NCHITESTS = Number of chi-square tests to do. Must be 2 or more. C EBINS = Expected value for each bin. EBINS = NBALLS/NBINS C SEED = Initial seed value for random number generator C C Note 1: The maximum size of array NBINS may be determined by the users C page file quota (pgflquo in AUTHORIZE). Also, it is recommended C the user have a very large working set quota (wsquo,wsextent) C to reduce page faulting. This can greatly improve speed. C We use INTEGER*2 array here to save space. C C Note 2: The maximum number of Chi-Square tests that can be saved is C arbitrary (array SAVEPROB). The user may choose any value. C C Note 3: The user may choose any starting seed value. Some restrictions C may apply depending upon the particular random number generator C being used. For example, RANDU should always be started with C an odd value of SEED. MTH$RANDOM may be started with any value C of SEED. C C Note 4: The MTH$RANDOM generator is used by the VAX FORTRAN intrinsic C function RAN and the VAX BASIC function RND. It is defined as: C C SEED = 69069*SEED + 1 mod 2**32 C X = SEED/2**32 C C Note 5: The RANDU generator is obsolite due to very strong correlation C in 3d space. ( Prove to yourself using 65539 = 2**16 + 3 that C SEED[i+2] = 6*SEED[i+1] - 9*SEED[i] ). It is defined as: C C SEED = 65539*SEED mod 2**31 C X = SEED/2**31 C C The RANDU generator should be started with an odd value of SEED. C C Note 6: The C standard library function rand() is defined as: C C SEED = 1103515245*SEED + 12345 mod 2**32 C IX = SEED mod 2**31 C C This standard random number generator is defined in the book: C The C Programming Language C Brian W. Kernighan and Dennis M. Ritchie C Prentice Hall, 1978 C C The same generator is defined in the ANSI C version by the same C authors above, and the same generator is used in VAX C. C C Note 7: The Microsoft C version 4.0 library function rand() impliments C the following: C C SEED = 214013*SEED + 2531011 mod 2**32 C IX = bits 16-31 of SEED C C Note 8: The Turbo Pascal version 6.0 function impliments the following: C C SEED = 134775813*SEED + 1 mod 2**32 C IX = bits 16-32 of SEED C IMPLICIT NONE INTEGER NDIM !Number of dimensions INTEGER NBINS !Number of bins INTEGER NBINSPD !Number of bins per dimension INTEGER NBALLS !Number of random numbers per !chi-square test. INTEGER NCHITESTS !Number of chi-square tests to do C INTEGER*2 BINS(20 000 000) !The bins. (see note 1) INTEGER*2 BINS(200 000) !Less bins. (see note 1) REAL EBINS !Expected number of balls per bin REAL SAVEPROB(100) !Array to save results of !chi-square tests (see note 2) INTEGER*4 SEED(2) !Only SEED(1) result is ever used. INTEGER*2 W(4) !Seeds for RANDU EQUIVALENCE(SEED,W) !for RANDU COMMON / SEEDSTORE / SEED INTEGER I,J,K,MRANDO,NBYTES,CLEAR,CHITEST,INDEX,IRANDOM,NCLEAR REAL*4 FRANDOM, RRANDOM CHARACTER*8 TIMEBUF EQUIVALENCE (IRANDOM,FRANDOM) REAL FOR$IRAN !The RANDU random number generator REAL RAND !UNIX rand() INTEGER*4 xrand !BSD random() REAL RAN1 !test generator supplied REAL D INTEGER JISHFT,IRANDOM2,COUNT REAL RANDES !DES FUNCTION (not supplied) INTEGER KEY(2) REAL CHSQ,PROB !Chi-square value, !chi-square probability REAL*4 FNBINSPD !float(NBINSPD) REAL*4 TWO31F REAL*4 TWO16F REAL*4 TWO15F TWO31F = 2.0**31.0 TWO16F = 2.0**16.0 TWO15F = 2.0**15.0 C*DES KEY(1) = 12345 !Choose any number you want C*DES KEY(2) = 678901 !to initialize DES with C*DES CALL DES_INIT(KEY) !DES code not included. 104 FORMAT(' Input number of dimensions NDIM: ',$) 100 FORMAT(' Input number of bins per dimension NBINSPD: ',$) 105 FORMAT(' Total number of bins = NBINSPD**NDIM = ',I) 106 FORMAT(' Minimum number of balls = 5*NBINS = ',I) 101 FORMAT(' Input total number of balls NBALLS: ',$) 103 FORMAT(' Input number of Chi-Square tests NCHITESTS (min=2) : ',$) 102 FORMAT(' Choose random number generator to test'/, 1 ' /*(1)*/ xrand(),'/ 1 ' /*(2)*/ UNIX rand(),'/ 1 ' /*(3)*/ MTH$RANDOM,'/ 2 ' /*(4)*/ RANDU,'/ 3 ' /*(5)*/ ANSI C,'/ 4 ' /*(6)*/ Microsoft C'/ 5 ' /*(7)*/ Turbo Pascal'/ 7 ' /*(8)*/ DES'/ 8 ' (9) another random number generator (choose this one)'/ 6 ' : ',$) 107 FORMAT(' Input starting SEED value: ',$) 200 FORMAT(BN,I) C ***GET USER INPUT*** 10 WRITE(6,104) !Input number of dimensions READ(5,200) NDIM WRITE(6,100) !Input number of bins per dimension READ(5,200) NBINSPD FNBINSPD = FLOAT(NBINSPD) NBINS = NBINSPD**NDIM !Calculate total number of bins WRITE(6,105) NBINS !Total number of bins is... WRITE(6,106) 5*NBINS !Minimum number of balls is... WRITE(6,101) !Input total number of balls READ(5,200) NBALLS WRITE(6,103) !Input number of chi-square tests to do READ(5,200) NCHITESTS WRITE(6,102) !Choose random number generator to test READ(5,200) MRANDO WRITE(6,107) !Starting SEED value READ(5,200) SEED(1) SEED(2) = 1 !Used only if random number generator !uses bigger than 32 bits C INITIALIZE GENERATOR IF NEEDED C*XRAND CALL sxrand(SEED(1)) !Initialize xrand() CALL RAN1(-SEED(1)) !Initialize RAN1 generator C Calculate expected average number of balls for each bin EBINS = FLOAT(NBALLS)/FLOAT(NBINS) C CALL TIME(TIMEBUF) C WRITE(6,201) TIMEBUF C201 FORMAT(1X,A8) DO CHITEST=1,NCHITESTS C *** ZERO BIN ARRAY *** DO I=1,NBINS BINS(I) = 0 ENDDO C*VMS !Quickly set BINS(k) = 0, k=1,...NBINS C*VMS !Does the equivalent of above C*VMS !but a lot faster. C*VMS K = 1 C*VMS NBYTES = NBINS*2 !total number of bytes to zero C*VMS DO WHILE (NBYTES .GT. 0) C*VMS IF (NBYTES .LE. 65534) THEN !maximum number of bytes we can clear C*VMS NCLEAR = NBYTES !in one call to LIB$MOVC5 is 65535 C*VMS ELSE C*VMS NCLEAR = 65534 !max that LIB$MOVC3 can do in one call C*VMS ENDIF !(make nclear an even number so we can divide evenly by 2) C*VMS CALL LIB$MOVC5(0,0,0,NCLEAR,BINS(K)) !Clear a block of memory C*VMS NBYTES = NBYTES - NCLEAR !Number of bytes still left to clear C*VMS K = K + NCLEAR/2 !Number of bytes cleared so far + 1 C*VMS ENDDO C Main Loop DO J=1,NBALLS INDEX = 1 DO I=0,NDIM-1 C ***PLACE YOUR RANDOM NUMBER GENERATOR HERE*** C Set IRANDOM using whatever random number generator you choose C IRANDOM = integer between 0 and NBINS-1 c IF (MRANDO .EQ. 1) THEN c IRANDOM = INT( ( float( xrand() ) /TWO31F ) *FNBINSPD) c ELSEIF (MRANDO .EQ. 2) THEN c IRANDOM = INT( RAND(SEED(1)) *FNBINSPD) !UNIX rand() c ELSEIF (MRANDO .EQ. 3) THEN c IRANDOM = INT( RAN(SEED(1)) *FNBINSPD) !VMS mth$random c ELSEIF (MRANDO .EQ. 4) THEN c IRANDOM = INT( FOR$IRAN(W(2),W(1)) *FNBINSPD) !Infamous randu c ELSEIF (MRANDO .EQ. 5) THEN c CALL LIB$EMUL(1103515245,SEED,12345,SEED) !ANSI C c IRANDOM = SEED(1) .AND. '7FFFFFFF'X c IRANDOM = INT( FLOAT(IRANDOM)/(TWO31F) *FNBINSPD) c ELSEIF (MRANDO .EQ. 6) THEN c CALL LIB$EMUL(214013,SEED,2531011,SEED) !Microsoft C 4.0 c IRANDOM = W(2) .AND. '7FFF'X c IRANDOM = INT( FLOAT(IRANDOM)/(TWO15F) *FNBINSPD) c ELSEIF (MRANDO .EQ. 7) THEN c CALL LIB$EMUL(134775813,SEED,1,SEED) !Turbo Pascal 6.0 c IRANDOM = SEED(1) .AND. 'FFFF0000'X c IRANDOM = JISHFT(IRANDOM,-16) c IRANDOM = INT( FLOAT(IRANDOM)/(TWO16F) * FNBINSPD) c ELSEIF (MRANDO .EQ. 8) THEN c IRANDOM = INT( RANDES() * FNBINSPD ) !DES (not supplied) c ELSEIF (MRANDO .EQ. 9) THEN IRANDOM = INT( RAN1(SEED(1)) * FNBINSPD ) c ENDIF C Calculate index by hand. INDEX = INDEX + IRANDOM*(NBINSPD**I) ENDDO BINS(INDEX) = BINS(INDEX) + 1 !ball fell in this bin C IF ( MOD(J, 1 000 000) .EQ. 0 ) THEN C CALL TIME(TIMEBUF) C WRITE(6,302) J, TIMEBUF 302 FORMAT(1X,'AT BALL:',I,3X,'TIME=',A8) C WRITE(6,303) SEED(2), SEED(1) 303 FORMAT(1X,'HEX: SEED(2)= ',Z,' SEED(1)= ',Z) C WRITE(6,304) SEED(2), SEED(1) 304 FORMAT(1X,'DEC: SEED(2)= ',I,' SEED(1)= ',I) C ENDIF ENDDO 400 CALL CHSONE(BINS,EBINS,NBINS,CHSQ,PROB) SAVEPROB(CHITEST) = PROB WRITE(6,1) NBALLS,CHSQ,PROB 1 FORMAT(' BALLS=',I,' CHISQ=',F,' PROB=',F) ENDDO C Now see if all the chi-square values are chi-square distributed: IF (NCHITESTS .GT. 1) THEN CALL KSONE(SAVEPROB,NCHITESTS,D,PROB) WRITE(6,2) D,PROB 2 FORMAT(1X,'KS D=',F,' PROB=',F) ENDIF END C============================================================================ C From book NUMERICAL RECIPES: The Art of Scientific Computing C Here for demonstration purposes C Replace this with whatever random number generator you want to test C Initialize with negative number FUNCTION RAN1(IDUM) REAL R(97) SAVE R !(Some UNIX F77 compilers require -save option on compile) PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6) PARAMETER (M3=243000,IA3=4561,IC3=51349) DATA IFF /0/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 IX1=MOD(IC1-IDUM,M1) IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IX1,M2) IX1=MOD(IA1*IX1+IC1,M1) IX3=MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 11 CONTINUE IDUM=1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE write(1,100) R write(1,102) R(J) 100 format(f) 102 format(1x,'RAN1 = ', F) RAN1=R(J) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 RETURN END C---------------------------------------------------------------------------- C CALCULATE THE CHI-SQUARE PROBABILITY. SINCE NBINS IS LARGE, IT IS JUST C THE CUMULATIVE GAUSSIAN DISTRIBUTION AFTER WE NORMALIZE THE VARIABLES. C OR ERROR FUNCTION. FUNCTION CHIPROB(NBINS,CHISQ) C Formula is the inverse of one given in Knuth for going the other way. INTEGER NBINS,DF REAL*4 CHISQ,Z DF = NBINS-1 Z = ( SQRT(24.0*CHISQ - 6.0*DF + 16.0) - 3*SQRT(2.0*DF) ) / 4.0 CHIPROB = ERF(Z) RETURN END FUNCTION ERF(X) C Return approximation to the complimentary error function erfc(X). C Return is not normalized err function. See book for details. C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C Modified to return normalized error function erf(X) C (It's a polynomial approximation) REAL ERFCC,Z,T Z=ABS(X/1.414213) !Normalize T=1./(1.+0.5*Z) ERFCC=T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+ * T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+ * T*(1.48851587+T*(-.82215223+T*.17087277))))))))) IF (X.LT.0.) ERFCC=2.-ERFCC ERF = 1.0 - ERFCC/2.0 !Normalize and compliment RETURN END C---------------------------------------------------------------------------- C THE FOLLOWING SUBROUTINES CALCULATE THE CHI-SQUARE VALUE: SUBROUTINE CHSONE(BINS,EBINS,NBINS,CHSQ,PROB) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing INTEGER NBINS INTEGER*2 BINS(NBINS) REAL EBINS,CHSQ,PROB CHSQ=0. IF(EBINS.LE.0.) PAUSE 'CHSONE: EBINS must be > 0' DO 11 J=1,NBINS CHSQ=CHSQ+(BINS(J)-EBINS)**2/EBINS 11 CONTINUE PROB=CHIPROB(NBINS,CHSQ) RETURN END C============================================================================ C THE FOLLOWING SUBROUTINES CALCULATE THE KOLMOGOROV-SMIRNOV PROBABILITY SUBROUTINE KSONE(DATA,N,D,PROB) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C DF - degrees of freedom. Passsed to FUNC INTEGER N REAL DATA(N) REAL D,PROB CALL PIKSRT(N,DATA) EN=N D=0. FO=0. DO 11 J=1,N FN=J/EN FF=DATA(J) DT=AMAX1(ABS(FO-FF),ABS(FN-FF)) IF(DT.GT.D)D=DT FO=FN 11 CONTINUE PROB=PROBKS(SQRT(EN)*D) RETURN END C---------------------------------------------------------------------------- FUNCTION PROBKS(ALAM) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C Note the routine in the Numerical Recipes book erronously returns C 1 instead of 0 for large values of ALAM. PARAMETER (EPS1=0.001, EPS2=1.E-8) A2=-2.*ALAM**2 FAC=2. PROBKS=0. TERMBF=0. DO 11 J=1,100 TERM=FAC*EXP(A2*J**2) PROBKS=PROBKS+TERM C Error in Numerical Recipes book. Terminate if TERM underflows. C** IF(ABS(TERM).LT.EPS1*TERMBF.OR.ABS(TERM).LT.EPS2*PROBKS)RETURN IF(ABS(TERM).LE.EPS1*TERMBF.OR.ABS(TERM).LE.EPS2*PROBKS)RETURN FAC=-FAC TERMBF=ABS(TERM) 11 CONTINUE PROBKS=1.0 RETURN END C---------------------------------------------------------------------------- SUBROUTINE PIKSRT(N,ARR) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C See book for details. INTEGER N REAL ARR(N) DO 12 J=2,N A=ARR(J) DO 11 I=J-1,1,-1 IF(ARR(I).LE.A)GO TO 10 ARR(I+1)=ARR(I) 11 CONTINUE I=0 10 ARR(I+1)=A 12 CONTINUE RETURN END

2 1

(fwd) Random Numbers - CORELA.FOR
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!agate!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd(a)netcom.com Subject: Random Numbers - CORELA.FOR Message-ID: <deleydCsIBCu.KzB(a)netcom.com> Organization: NETCOM On-line Communication Services (408 261-4700 guest) Date: Wed, 6 Jul 1994 06:58:06 GMT Lines: 211 Xref: bga.com sci.stat.math:1317 sci.math:15355 sci.math.num-analysis:3356 {Approx. 200 lines} PROGRAM CORELA C Perform a KS test comparing the first 100 elements of a random C number generator, starting with SEED values of 1..10 C C This is not polished code you put on a shelf and admire, this is code C you dig your hands into and work with to make it do what you want it to. C C Place your random number generator to be tested where it C bluntly says "PLACE YOUR RANDOM NUMBER GENERATOR HERE". C The output is a floating point between 0 (inclusive) and 1 (exclusive). C C Currently the program is set up to use subroutine RAN1, a portable C random number generator from the book "NUMERICAL RECIPES: The Art C of Scientific Computing". I've had trouble on our UNIX system C not making array R(97) static even though the code says to. C Compiling with the -static qualifier works. C IMPLICIT NONE INTEGER SINC,I,J REAL SEQ(100,10) REAL AR(10) REAL D,PROB INTEGER MRANDO, SEEDINIT, IRANDOM INTEGER*4 SEED(2) !Only SEED(1) result is ever used. INTEGER*2 W(4) !Seeds for RANDU EQUIVALENCE(SEED,W) !for RANDU COMMON / SEEDSTORE / SEED REAL*4 FRANDOM REAL FOR$IRAN !The RANDU random number generator REAL RAN1 !test generator supplied INTEGER JISHFT,IRANDOM2,COUNT REAL*4 FNBINSPD !float(NBINSPD) REAL*4 TWO31F REAL*4 TWO16F REAL*4 TWO15F TWO31F = 2.0**31.0 TWO16F = 2.0**16.0 TWO15F = 2.0**15.0 102 FORMAT(' Choose random number generator to test'/, 1 ' /*(1)*/ MTH$RANDOM,'/ 2 ' /*(2)*/ RANDU,'/ 3 ' /*(3)*/ ANSI C,'/ 4 ' /*(4)*/ Microsoft C'/ 5 ' /*(5)*/ Turbo Pascal'/ 8 ' (9) another random number generator (choose this one)'/ 6 ' : ',$) 107 FORMAT(' Input starting SEED value: ',$) 108 FORMAT(' Input increment between SEED values: ',$) 200 FORMAT(BN,I) 10 CONTINUE WRITE(6,102) !Choose random number generator to test READ(5,200) MRANDO WRITE(6,107) !Starting SEED value READ(5,200) SEED(1) SEEDINIT = SEED(1) SEED(2) = 1 WRITE(6,108) !INCREMENT VALUE READ(5,200) SINC C Main Loop DO J=1,10 !10 sequences DO I=1,100 !sequence length of first 100 numbers C ***PLACE YOUR RANDOM NUMBER GENERATOR HERE*** C Set FRANDOM using whatever random number generator you choose C to a floating point value in the range [0,1) FRANDOM = RAN1(SEED(1)) C IF (MRANDO .EQ. 1) THEN C FRANDOM = RAN(SEED(1)) !mth$random C ELSEIF (MRANDO .EQ. 2) THEN C FRANDOM = FOR$IRAN(W(2),W(1)) !randu C ELSEIF (MRANDO .EQ. 3) THEN C CALL LIB$EMUL(1103515245,SEED,12345,SEED) !VAX C C IRANDOM = SEED(1) .AND. '7FFFFFFF'X C FRANDOM = FLOAT(IRANDOM)/(TWO31F) C ELSEIF (MRANDO .EQ. 4) THEN C CALL LIB$EMUL(214013,SEED,2531011,SEED) !Microsoft C 4.0 C IRANDOM = W(2) .AND. '7FFF'X C FRANDOM = FLOAT(IRANDOM)/(TWO15F) C ELSEIF (MRANDO .EQ. 5) THEN C CALL LIB$EMUL(134775813,SEED,1,SEED) !Turbo Pascal 6.0 C IRANDOM = SEED(1) .AND. 'FFFF0000'X C IRANDOM = JISHFT(IRANDOM,-16) C FRANDOM = FLOAT(IRANDOM)/(TWO16F) C ENDIF SEQ(I,J) = FRANDOM ENDDO SEEDINIT = SEEDINIT + SINC !calculate new initial seed SEED(1) = SEEDINIT !set new initial seed ENDDO C Do a KS test on each edlement comparing the 10 sequences DO I=1,100 DO J=1,10 AR(J) = SEQ(I,J) !Transfer to short array ENDDO CALL KSONE(AR,10,D,PROB) WRITE(6,2) I,PROB 2 FORMAT(1X,'I=',I4,' KS PROB=',F) ENDDO END C============================================================================ C From book NUMERICAL RECIPES: The Art of Scientific Computing C Here for demonstration purposes C Replace this with whatever random number generator you want to test C Initialize with negative number FUNCTION RAN1(IDUM) DIMENSION R(97) SAVE R !(Some UNIX F77 compilers require -save option on compile) PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6) PARAMETER (M3=243000,IA3=4561,IC3=51349) DATA IFF /0/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 IX1=MOD(IC1-IDUM,M1) IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IX1,M2) IX1=MOD(IA1*IX1+IC1,M1) IX3=MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 11 CONTINUE IDUM=1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE RAN1=R(J) R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1 RETURN END C============================================================================== C THE FOLLOWING SUBROUTINES CALCULATE THE KOLMOGOROV-SMIRNOV PROBABILITY SUBROUTINE KSONE(DATA,N,D,PROB) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C DF - degrees of freedom. Passsed to FUNC INTEGER N REAL DATA(N) REAL D,PROB CALL PIKSRT(N,DATA) EN=N D=0. FO=0. DO 11 J=1,N FN=J/EN FF=DATA(J) DT=AMAX1(ABS(FO-FF),ABS(FN-FF)) IF(DT.GT.D)D=DT FO=FN 11 CONTINUE PROB=PROBKS(SQRT(EN)*D) RETURN END C------------------------------------------------------------------------------ FUNCTION PROBKS(ALAM) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C Note the routine in the Numerical Recipes book erronously returns C 1 instead of 0 for large values of ALAM. PARAMETER (EPS1=0.001, EPS2=1.E-8) A2=-2.*ALAM**2 FAC=2. PROBKS=0. TERMBF=0. DO 11 J=1,100 TERM=FAC*EXP(A2*J**2) PROBKS=PROBKS+TERM C Error in Numerical Recipes book. Terminate if TERM underflows. C** IF(ABS(TERM).LT.EPS1*TERMBF.OR.ABS(TERM).LT.EPS2*PROBKS)RETURN IF(ABS(TERM).LE.EPS1*TERMBF.OR.ABS(TERM).LE.EPS2*PROBKS)RETURN FAC=-FAC TERMBF=ABS(TERM) 11 CONTINUE PROBKS=1.0 RETURN END C------------------------------------------------------------------------------ SUBROUTINE PIKSRT(N,ARR) C Adapted from book NUMERICAL RECIPES: The Art of Scientific Computing C See book for details. INTEGER N REAL ARR(N) DO 12 J=2,N A=ARR(J) DO 11 I=J-1,1,-1 IF(ARR(I).LE.A)GO TO 10 ARR(I+1)=ARR(I) 11 CONTINUE I=0 10 ARR(I+1)=A 12 CONTINUE RETURN END

1 0

(fwd) Random Numbers - SPECTRAL.FOR
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.stat.math,sci.math,sci.math.num-analysis Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!swrinde!ihnp4.ucsd.edu!agate!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd(a)netcom.com Subject: Random Numbers - SPECTRAL.FOR Message-ID: <deleydCsIBEG.L2B(a)netcom.com> Organization: NETCOM On-line Communication Services (408 261-4700 guest) Date: Wed, 6 Jul 1994 06:59:04 GMT Lines: 426 Xref: bga.com sci.stat.math:1318 sci.math:15356 sci.math.num-analysis:3357 {Approx. 420 lines} PROGRAM SPECTRAL ! Performs the spectral test for a linear congruential random number generator. ! ! This program adapted from: ! ! ALGORYTHM AS 193 APPLIED STATISTICS, (1983) VOL. 32, NO.3 PG. 328-335 ! T. R. Hopkins ! Modified to run on VAX/VMS systems using REAL*16 (64 bit) variables. ! The original is FORTRAN-66 compliant. ! ! Consider linear congruential generators of the form: ! SEED = (A*SEED + C) mod M ! ! Given A, and M, the spectral test calculates NUSQ (NU**2), LOGNU (base 2), ! and MU. As a guide, Knuth suggests a multiplier A may be considered ! adequate if the values of MU returned by the spectral test are > 0.1 . ! For an exceptionally good multiplier, these values will all be greater ! than unity. ! ! The spectral test may be applied if: ! 1. The sequence has maximal period, or ! 2. M is prime and C = 0 and the period length is M-1, or ! 3. M = 2**e and A mod 8 = 5 or A mod 8 = 3. ! In this third case the spectral test is applied using ! A = A and M = 2**(e-2). For example, in analyzing RANDU, ! use A = 65539 and M = 536870912 (2**29). ! ! Further information on the spectral test is in: ! ! Knuth, Donald E. "The Art of Computer Programming Vol. 2: Seminumerical ! algorithms, 2nd edition. Reading, Mass.: Addison-Wesley. 1981 ! ! The value of parameter BIGT determines how many dimensions are calculated. ! Higher dimensions may be obtained by changing this parameter and recompiling. ! Note that 12 is about the highest feasible. Above 12 the program may take ! days to complete. C Example: C MTH$RANDOM is defined as C C SEED = (69069*SEED + 1) MOD 2**32 C C Here A = 69069 C and M = 2**32 = 4294967296 C C $ RUN SPECTRAL C INPUT A: 69069 C INPUT M: 4294967296 C C A= 69069.0 C M= 4294967296.0 C BIGT= 6 C NUSQ= C NUSQ ( 2)= 4243209856.000000 C NUSQ ( 3)= 2072544.000000 C NUSQ ( 4)= 52804.000000 C NUSQ ( 5)= 6990.000000 C NUSQ ( 6)= 242.000000 C LOGNU= C LOGNU( 2)= 15.991254 C LOGNU( 3)= 10.491486 C LOGNU( 4)= 7.844180 C LOGNU( 5)= 6.385538 C LOGNU( 6)= 3.959432 C MU= C MU= ( 2)= 3.103734 C MU= ( 3)= 2.909942 C MU= ( 4)= 3.203639 C MU= ( 5)= 5.006469 C MU= ( 6)= 0.017052 C C Now examine the MU values. All values are above 1 except the very last C value MU(6) is 0.01, indicating MTH$RANDOM may not perform as well in a C 6-D test. C C C Run spectral again this time trying the values for the bad RANDU generator: C C MTH$RANDOM is defined as C C SEED = (65539*SEED) MOD 2**31 C C Here A = 65539 C and M = 2**31 but we use M = 2**29 for reasons discussed above C C C $ RUN SPECTRAL C INPUT A: 65539 C INPUT M: 536870912 !(2**29) C C A= 65539.0 C M= 536870912.0 C BIGT= 6 C NUSQ= C NUSQ ( 2)= 536936458.000000 C NUSQ ( 3)= 118.000000 C NUSQ ( 4)= 116.000000 C NUSQ ( 5)= 116.000000 C NUSQ ( 6)= 116.000000 C LOGNU= C LOGNU( 2)= 14.500088 C LOGNU( 3)= 3.441322 C LOGNU( 4)= 3.428990 C LOGNU( 5)= 3.428990 C LOGNU( 6)= 3.428990 C MU= C MU= ( 2)= 3.141976 C MU= ( 3)= 0.000010 C MU= ( 4)= 0.000124 C MU= ( 5)= 0.001421 C MU= ( 6)= 0.015025 C C Notice here the MU values for dimensions 2 through 6 are all extremely C small. This generator does horribly on these dimensions. The spectral C test noticed it right away. PARAMETER BIGT = 6 !Number of dimensions to go up to. Max is 12. PARAMETER IU = BIGT !(Beyond 12 program may take days to run.) PARAMETER IV = BIGT INTEGER*4 IFAULT REAL*16 A, M, MU(BIGT), NUSQ(BIGT), LOGNU(BIGT), U(IU,BIGT), 2 V(IV,BIGT), Z(BIGT) 100 FORMAT(' INPUT A: ',$) 101 FORMAT(' INPUT M: ',$) 200 FORMAT(BN,G33.0) 201 WRITE(6,100) READ(5,200) A !MTH$RANDOM example: A = 69069.0 WRITE(6,101) READ(5,200) M !MTH$RANDOM example: M = 4294967296.0 (2**32) CALL SPECT(A,M,BIGT,MU,NUSQ,LOGNU,U,IU,V,IV,Z,IFAULT) IF (IFAULT .GT. 0) THEN IF (IFAULT .EQ. 1) THEN PRINT*, ' BIGT < 2' ELSEIF (IFAULT .EQ. 2) THEN PRINT*, ' A .GE. M .OR. A .LE. 0 .OR. M .LE. 0' ELSEIF (IFAULT .EQ. 3) THEN PRINT*, ' M > Mmax' ELSEIF (IFAULT .EQ. 4) THEN PRINT*, ' A and M not relatively prime' ELSEIF (IFAULT .EQ. 5) THEN PRINT*, ' Intermediate result > Mmax * Mmax' ELSE PRINT*, ' IFAULT .GT. 5' ENDIF STOP ENDIF WRITE(6,1) A WRITE(6,2) M WRITE(6,3) BIGT WRITE(6,41) DO I=2,BIGT WRITE(6,4) I,NUSQ(I) ENDDO WRITE(6,51) DO I=2,BIGT WRITE(6,5) I,LOGNU(I) ENDDO WRITE(6,61) DO I=2,BIGT WRITE(6,6) I,MU(I) ENDDO 1 FORMAT(' A=',F33.1) 2 FORMAT(' M=',F33.1) 3 FORMAT(' BIGT=',I) 41 FORMAT(' NUSQ=') 4 FORMAT(' NUSQ (',I,')=',F33.6) 51 FORMAT(' LOGNU=') 5 FORMAT(' LOGNU(',I,')=',F33.6) 61 FORMAT(' MU=') 6 FORMAT(' MU= (',I,')=',F33.6) C GOTO 201 END SUBROUTINE SPECT(A, M, BIGT, MU, NUSQ, LOGNU, U, IU, V, IV, Z, * IFAULT) C C ALGORYTHM AS 193 APPLIED STATISTICS, (1983) VOL. 32, NO.3 PG. 328-335 C T. R. Hopkins C C A REVISED ALGORITHM FOR THE SPECTRAL TEST C Modified to use REAL*16 variables for VAX/VMS C IMPLICIT NONE INTEGER*4 I, I2, J, K INTEGER*4 BIGT, IU, IV, T, T1, IFAULT REAL*16 A, M, MU(BIGT), NUSQ(BIGT), LOGNU(BIGT), * U(IU, BIGT), V(IV, BIGT), Z(BIGT), * H, HPRIME, MMAX, MMAX2, MSQ, P, PI, PPRIME, Q, * QTEMP, R, S, SIGN, UC, VC, VIJ, VJJ, W, ZERO, ONE, TWO, FOUR, * DINT, DNINT, VPROD DATA ZERO /0.0Q0/, ONE /1.0Q0/, TWO /2.0Q0/, FOUR /4.0Q0/ C C SUITABLE VALUES FOR C 1) IBM REAL*8 C DATA MMAX/33554432.0D0/ C 2) IBM REAL*16 C 3) CDC 7600 DOUBLE PRECISION C DATA MMAX/35184372088832.0D0/ C DATA MMAX /9007199254740992.0D0/ C C A VAX/VMS REAL*16 has precision approximately one part in 2**112 C Knuth claims values rarely if ever exceed M**2 C So Hopkins takes maxval = 8*m**2 and solves 2**112 = 8*m**2 for M C giving Mmax = 2**(112/2)/8 DATA MMAX /9.0Q15/ C C TEST THE VALIDITY OF THE INPUT PARAMETERS C MMAX2 = MMAX * MMAX IFAULT = 0 IF (BIGT .LT. 2) IFAULT = 1 IF (A .GE. M .OR. A .LE. ZERO .OR. M .LE. ZERO) IFAULT = 2 IF (M .GT. MMAX) IFAULT = 3 IF (IFAULT .GT. 0) RETURN C C CHECK A AND M ARE RELATIVELY PRIME C NEED VALID A AND M C USE EUCLIDS ALGORITHM C H = A HPRIME = M 10 R = QMOD(HPRIME, H) IF (R .EQ. ZERO) GOTO 20 HPRIME = H H = R GOTO 10 20 IF (H .NE. ONE) IFAULT = 4 ! A and M not relatively prime IF (IFAULT .NE. 0) RETURN MSQ = M * M C C ALL STEPS REFER TO THOSE IN KNUTHS ALGORITHM C STEP 1 - INITIALIZATION C H = A HPRIME = M P = ONE PPRIME = ZERO R = A S = ONE + A * A C C STEP 2 - EUCLIDEAN STEP C 30 Q = QINT(HPRIME / H) UC = HPRIME - Q * H VC = PPRIME - Q * P W = UC * UC + VC * VC IF (W .GE. S) GOTO 40 S = W HPRIME = H H = UC PPRIME = P P = VC GOTO 30 C C STEP 3 - COMPUTE NU(2) C 40 UC = UC - H VC = VC - P W = UC * UC + VC * VC IF (W .GE. S) GOTO 50 S = W HPRIME = UC PPRIME = VC 50 NUSQ(2) = S C C INITIALIZE U AND V MATRICES C NOTE WE STORE BY COLUMNS WHEREAS KNUTH STORES BY ROWS C T = 2 U(1, 1) = -H U(1, 2) = -HPRIME U(2, 1) = P U(2, 2) = PPRIME SIGN = ONE IF (PPRIME .GT. ZERO) SIGN = -ONE V(1, 1) = SIGN * PPRIME V(1, 2) = -SIGN * P V(2, 1) = SIGN * HPRIME V(2, 2) = -SIGN * H C C STEP 4 - ADVANCE T C 60 IF (T .EQ. BIGT) GOTO 200 T1 = T T = T + 1 R = QMOD(A * R, M) U(1, T) = -R U(T, T) = ONE U(T, 1) = ZERO V(1, T) = ZERO V(T, T) = M DO 70 I = 2, T1 U(I, T) = ZERO U(T, I) = ZERO V(I, T) = ZERO 70 CONTINUE DO 90 I = 1, T1 QTEMP = V(1, I) * R Q = QNINT(QTEMP / M) V(T, I) = QTEMP - Q * M DO 80 I2 = 1, T 80 U(I2, T) = U(I2, T) + Q * U(I2, I) 90 CONTINUE S = QMIN1(S, VPROD(U(1, T), U(1, T), T)) K = T J = 1 C C STEP 5 - TRANSFORM C 100 DO 120 I = 1, T IF (I .EQ. J) GOTO 120 VIJ = VPROD(V(1, I), V(1, J), T) VJJ = VPROD(V(1, J), V(1, J), T) IF (TWO * QABS(VIJ) .LE. VJJ) GOTO 120 Q = QNINT(VIJ / VJJ) DO 110 I2 = 1, T V(I2, I) = V(I2, I) - Q * V(I2, J) U(I2, J) = U(I2, J) + Q * U(I2, I) 110 CONTINUE K = J 120 CONTINUE C C STEP 6 - EXAMINE NEW BOUND C IF (K .EQ. J) S = QMIN1(S, VPROD(U(1, J), U(1, J), T)) C C STEP 7 - ADVANCE J C J = J + 1 IF (J .EQ. T + 1) J = 1 IF (J .NE. K) GOTO 100 C C STEP 8 - PREPARE FOR SEARCH C C MU AND LOGNU ARE USED TO STORE KNUTHS X AND Y RESPECTIVELY C DO 130 I = 1, T MU(I) = ZERO LOGNU(I) = ZERO QTEMP = VPROD(V(1, I), V(1, I), T) IF (QTEMP .GT. MMAX2) GOTO 240 !Intermediate result > Mmax * Mmax QTEMP = QTEMP / MSQ Z(I) = QINT(QSQRT(QINT(QTEMP * S))) 130 CONTINUE K = T C C STEP 9 - ADVANCE XK C 140 IF (MU(K) .EQ. Z(K)) GOTO 190 MU(K) = MU(K) + ONE DO 150 I = 1, T 150 LOGNU(I) = LOGNU(I) + U(I, K) C C STEP 10 - ADVANCE K C 160 K = K + 1 IF (K .GT. T) GOTO 180 MU(K) = -Z(K) DO 170 I = 1, T 170 LOGNU(I) = LOGNU(I) - TWO * Z(K) * U(I, K) GOTO 160 180 S = QMIN1(S, VPROD(LOGNU, LOGNU, T)) C C STEP 11 - DECREASE K C 190 K = K - 1 IF (K .GE. 1) GOTO 140 NUSQ(T) = S GOTO 60 C C CALCULATE NU AND LOG(NU) C 200 DO 210 I = 2, BIGT MU(I) = QSQRT(NUSQ(I)) LOGNU(I) = QLOG(MU(I)) / QLOG(TWO) 210 CONTINUE C C CALCULATE TRANSFORMED MU VALUES C PI = 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37511 Q = ONE DO 220 T = 2, BIGT, 2 Q = Q * PI * TWO / QEXT(T) MU(T) = Q * MU(T) ** T / M 220 CONTINUE IF (BIGT .EQ. 2) RETURN Q = TWO DO 230 T = 3, BIGT, 2 Q = Q * PI * TWO / QEXT(T) MU(T) = Q * MU(T) ** T / M 230 CONTINUE RETURN 240 IFAULT = 5 !Intermediate result > Mmax * Mmax RETURN END REAL*16 FUNCTION VPROD(U, V, T) C C ALGORYTHM AS 193 APPLIED STATISTICS, (1983) VOL. 32, NO.3 PG. 328-335 C C AUXILIARY FUNCTION TO CALCULATE THE INNER PRODUCT OF C THE TWO VECTORS U AND V OF LENGTH T. C Modified to REAL*16 C INTEGER T REAL*16 U(T), V(T), SUM, ZERO DATA ZERO /0.0Q0/ C SUM = ZERO DO 10 I = 1, T 10 SUM = SUM + U(I) * V(I) VPROD = SUM RETURN END

1 0

(fwd) Re: BSD random() - any good (source included)?
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.math,sci.stat.math Path: bga.com!news.sprintlink.net!hookup!europa.eng.gtefsd.com!howland.reston.ans.net!spool.mu.edu!agate!library.ucla.edu!csulb.edu!csus.edu!netcom.com!deleyd From: deleyd(a)netcom.com Subject: Re: BSD random() - any good (source included)? Message-ID: <deleydCs1L9A.2L6(a)netcom.com> Followup-To: sci.math,sci.stat.math Organization: NETCOM On-line Communication Services (408 261-4700 guest) X-Newsreader: TIN [version 1.2 PL1] References: <2ua9ln$4lv(a)news.tamu.edu> Date: Mon, 27 Jun 1994 06:12:45 GMT Lines: 13 Xref: bga.com sci.math:14916 sci.stat.math:1219 I did a research paper on Computer Generated Random Number Sequences in 1991. Included are the results of testing numerous popular generators. The code used for testing the generators is also available if one is so inclined to do some testing of a particular generator. (The only thing is a thorough test to determine the limits of the generator can take many hours of CPU time). Perhaps later this week I'll post the paper and see what the response is. I'm always a bit apprehensive to post. Never sure what the response will be. Maybe someone will think it's interesting. David Deley deleyd(a)netcom.com

1 0

(fwd) BSD random() - any good (source included)?
by Jim choate 17 Dec '03

17 Dec '03

Path: bga.com!news.sprintlink.net!news.onramp.net!convex!cs.utexas.edu!bcm!news.tamu.edu!henrik From: henrik(a)stat.tamu.edu (Henrik Schmiediche) Newsgroups: sci.math,sci.stat.math Subject: BSD random() - any good (source included)? Date: 22 Jun 1994 21:15:35 GMT Organization: Department of Statistics, Texas A&M University Lines: 140 Message-ID: <2ua9ln$4lv(a)news.tamu.edu> NNTP-Posting-Host: picard.tamu.edu Xref: bga.com sci.math:14740 sci.stat.math:1193 Hello, the BSD random() function returns a pseudo random number. I would like to know if anyone knows how good this random number generator is and if it has been thouroughly tested. Below are two descriptions of the generator for two different sources. Looking at the source code it is obvious that this generator is seeded using a linear congruetial generator that leaves much to be desired (low bits alternate). I remember reading somewhere that the trinomials used by random() are not optimal but I can't remember the source. The generator does have some great advantages like being very fast and having a very long period, but both these advantages are meaningless if the random numbers it produces are not very good. Anyone know more about random() and if it is any good? I have include a source code implementation below (I wrote it originally so I could inline the code into my own application which spend a significant amount of time generating random numbers). - henrik According to the SunOS doc's: "random () uses a non-linear additive feedback random number generator employing a default table of size 31 long integers to return successive pseudo-random numbers in the range from 0 to (2**31)-1. The period of this random number generator is very large, approximately 16*((2**31)-1)." The BSD source code (from glibc) says: "The random number generation technique is a linear feedback shift register approach, employing trinomials (since there are fewer terms to sum up that way). In this approach, the least significant bit of all the numbers in the state table will act as a linear feedback shift register, and will have period 2^deg - 1 (where deg is the degree of the polynomial being used, assuming that the polynomial is irreducible and primitive). The higher order bits will have longer periods, since their values are also influenced by pseudo-random carries out of the lower bits. The total period of the generator is approximately deg*(2**deg - 1); thus doubling the amount of state information has a vast influence on the period of the generator." For table size of 31 long ints random() use the trinomial: x**31 + x**3 + 1. For 63 long ints it uses the trinomial x**63 + x + 1. ***************************************************************************** /*** Code to implement random() & srandom() of BSD Unix. It was taken (though coded somewhat differently) from the Gnu BSD implementation. ***/ #include <stdio.h> #include <stdlib.h> #ifdef LONG31 /* x^31 + x^3 + 1 */ #define SIZE 31 #define SIZE1 30 #define P1 3 #define P2 0 #else /* LONG63: x^63 + x + 1 */ #define SIZE 63 #define SIZE1 62 #define P1 1 #define P2 0 #endif #define LONG_MAX 0x7fffffff int p1=P1, p2=P2; long table[SIZE]; /*** return a "random" number in range [0, LONG_MAX] */ long xrand () { int r; table[p1] = table[p1] + table[p2]; /* add two table elements */ r = (table[p1] >> 1) & LONG_MAX; /* throw least significant bit away */ if (p1 == SIZE1) { /* increment the table indexes */ p1 = 0; p2 = p2 + 1; } else if (p2 == SIZE1) { p1 = p1 + 1; p2 = 0; } else { p1 = p1 + 1; p2 = p2 + 1; } return (r); } /*** use a linear congruential type generator to seed the state table & cycle the entire table 10 times */ void sxrand (seed) long seed; { int i; table[0] = seed; for (i=1; i<SIZE; ++i) table[i] = (table[i-1] * 1103515145) + 12345; /* lousy */ for (i=0; i<10*SIZE; ++i) (void) xrand(); } /*** a small test program */ void main () { int i; sxrand (1); /* BSD default */ for (i=1; i<=40; ++i) printf ("%ld", xrand() % 10 ); /* least random bits ? */ /* 6714066113586447326208220248220881760069 (cc -DLONG63) */ /* 9418752338157675324663485137890734831064 (cc -DLONG31) */ printf ("\n"); } -- Henrik Schmiediche, Dept. of Statistics, Texas A&M, College Station, TX 77843 E-mail: henrik(a)stat.tamu.edu | Tel: (409) 862-1764 | Fax: (409) 845-3144 Finger for pgp 2.5 key, fingerprint: CE8F BD6C 59FC DA85 376A BB96 2E83 FF5E

1 0

(fwd) Re: BSD random() - any good (source included)?
by Jim choate 17 Dec '03

17 Dec '03

Path: bga.com!news.sprintlink.net!hookup!europa.eng.gtefsd.com!newsxfer.itd.umich.edu!nntp.cs.ubc.ca!mala.bc.ca!hakatac!rbursey Newsgroups: sci.math,sci.stat.math Subject: Re: BSD random() - any good (source included)? Message-ID: <iFwLoc3w165w(a)hakatac.almanac.bc.ca> From: rbursey(a)hakatac.almanac.bc.ca (robert bursey) Date: Mon, 27 Jun 94 18:24:29 PDT References: <deleydCs1L9A.2L6(a)netcom.com> Distribution: na Organization: Sir HackAlot's UNIX BBS, Port Alberni, B.C. Lines: 18 Xref: bga.com sci.math:14978 sci.stat.math:1239 deleyd(a)netcom.com writes: > I did a research paper on Computer Generated Random Number Sequences in > 1991. Included are the results of testing numerous popular generators. > The code used for testing the generators is also available if one > is so inclined to do some testing of a particular generator. (The only > thing is a thorough test to determine the limits of the generator can > take many hours of CPU time). > > Perhaps later this week I'll post the paper and see what the response > is. I'm always a bit apprehensive to post. Never sure what the > response will be. Maybe someone will think it's interesting. > > David Deley > deleyd(a)netcom.com Does anybody know of a good test for randomness? I would definitely like to know how good computer RNG's are. Post away!

1 0

(fwd) Re: BSD random() - any good (source included)?
by Jim choate 17 Dec '03

17 Dec '03

Newsgroups: sci.math,sci.stat.math Path: bga.com!news.sprintlink.net!hookup!ames!lll-winken.llnl.gov!overload.lbl.gov!dog.ee.lbl.gov!ihnp4.ucsd.edu!usc!howland.reston.ans.net!europa.eng.gtefsd.com!MathWorks.Com!news.kei.com!ssd.intel.com!carr From: carr(a)ssd.intel.com (George Carr) Subject: Re: BSD random() - any good (source included)? Message-ID: <Cs4B6w.IqE(a)SSD.intel.com> Sender: usenet(a)SSD.intel.com Nntp-Posting-Host: shiva Organization: Supercomputer Systems Divison, Intel Corp. References: <deleydCs1L9A.2L6(a)netcom.com> <iFwLoc3w165w(a)hakatac.almanac.bc.ca> Distribution: na Date: Tue, 28 Jun 1994 17:28:08 GMT Lines: 34 Xref: bga.com sci.math:14998 sci.stat.math:1241 In article <iFwLoc3w165w(a)hakatac.almanac.bc.ca>, rbursey(a)hakatac.almanac.bc.ca (robert bursey) writes: |> deleyd(a)netcom.com writes: |> |> > I did a research paper on Computer Generated Random Number Sequences in |> > 1991. Included are the results of testing numerous popular generators. |> > The code used for testing the generators is also available if one |> > is so inclined to do some testing of a particular generator. (The only |> > thing is a thorough test to determine the limits of the generator can |> > take many hours of CPU time). |> > |> > Perhaps later this week I'll post the paper and see what the response |> > is. I'm always a bit apprehensive to post. Never sure what the |> > response will be. Maybe someone will think it's interesting. |> > |> > David Deley |> > deleyd(a)netcom.com |> |> Does anybody know of a good test for randomness? I would definitely like to |> know how good computer RNG's are. Post away! The classic reference is Volume 2 of Donald Knuth's The Art of Computer Programming, Second Edition, Seminumerical Algorithms. I highly recommend it to anyone wanting to know what "random" is all about. If you really need to know whether your generator is random-enough for your application you should expect to do your own testing and yes it will require many hours of your time in addition to that of your computer. -- George R. Carr, Jr. Intel Supercomputer Systems Division Parallel Systems Engineer NOAA Forecast Systems Laboratory carr(a)ssd.intel.com Research Lab 3, Rm A227 carr(a)neko.fsl.noaa.gov 3100 Marine St, R/E/FS5 voice: 303-497-6130 Boulder, CO 80303 fax: 303-497-6821

1 0

(fwd) Re: BSD random() - any good (source included)?
by Jim choate 17 Dec '03

17 Dec '03

Path: bga.com!news.sprintlink.net!hookup!ames!lll-winken.llnl.gov!noc.near.net!pad-thai.aktis.com!la-jiao.aktis.com!not-for-mail From: don(a)cam.ov.com (Donald T. Davis) Newsgroups: sci.math,sci.stat.math Subject: Re: BSD random() - any good (source included)? Date: 28 Jun 1994 17:52:48 -0400 Organization: OpenVision Technologies, Inc. Lines: 65 Distribution: na Message-ID: <2uq63g$g5c(a)la-jiao.aktis.com> References: <deleydCs1L9A.2L6(a)netcom.com> <iFwLoc3w165w(a)hakatac.almanac.bc.ca> <Cs4B6w.IqE(a)SSD.intel.com> NNTP-Posting-Host: la-jiao.aktis.com Xref: bga.com sci.math:15008 sci.stat.math:1243 (George Carr) writes: >(robert bursey) writes: >|> deleyd(a)netcom.com writes: >|> >|> > I did a research paper on Computer Generated Random Number Sequences in >|> > 1991. Included are the results of testing numerous popular generators. >|> > The code used for testing the generators is also available if one >|> > is so inclined to do some testing of a particular generator. (The only >|> > thing is a thorough test to determine the limits of the generator can >|> > take many hours of CPU time). >|> > >|> > Perhaps later this week I'll post the paper and see what the response >|> > is. I'm always a bit apprehensive to post. Never sure what the >|> > response will be. Maybe someone will think it's interesting. >|> > >|> > David Deley >|> > deleyd(a)netcom.com >|> >|> Does anybody know of a good test for randomness? I would definitely like to >|> know how good computer RNG's are. Post away! > >The classic reference is Volume 2 of Donald Knuth's The Art of Computer >Programming, Second Edition, Seminumerical Algorithms. I highly recommend >it to anyone wanting to know what "random" is all about. > >If you really need to know whether your generator is random-enough for >your application you should expect to do your own testing and yes it will >require many hours of your time in addition to that of your computer. >-- knuth's chapter's practical results are about linear-congruential rngs, their optimization and testing. though these rngs are still distressingly common, nonlinear rngs are the way to go for two burgeoning areas that consume random numbers: graphics and cryptography. both areas are concerned with getting extremely long periods, but cryptography is also concerned with proving unpredictability of secure rngs. that is, knowing some outputs of an rng as applied to a given seed, it should be impossible to deduce or predict other outputs' values. so, you see, the "good test for randomness" depends strongly on which features of a random variable you want to use. if you're careful, knuth's approach will work fine for some statistical applications, like monte-carlo techniques. but knuth's is by no means the last word on the subject. btw, for cryptographic purposes, the received wisdom is that there is NO adequate test for randomness; if an rng passes lots of tests, that's very nice, but the presumption is that the variable's deterministic structure is simply hidden, and that the clever-enough test was not yet applied or devised. nevertheless, in the cryptographic field, the list of tests used is long. typically, you design the test to probe the weaknesses of a specific rng algorithm. period tests, runs tests, and substring-interarrival tests are common, and some people like entropy estimates. one of the trusted names in the crypto-rng literature is marsaglia; he has published extensively on the subject of rng-testing. i don't know what the graphics literature on rngs is like; i only know that if you want to simulate textured surfaces, like grass, a bad rng makes a striped texture. be forwarned: the rng literature is amazingly vast, with a low signal- to-noise ratio. it seems that everyone thinks he can design a "good" rng. btw, i favor hardware rngs. -don davis openvision technologies cambridge, ma

1 0