December 2003 - cypherpunks-legacy

(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

Re: (fwd) Re: BSD random() - any good (source included)?
by jamiel＠sybase.com 17 Dec '03

17 Dec '03

At 3:20 PM 07/07/94 -0500, Jim choate wrote 69K of forwards (so far) that I would much rather look at in sci.math or sci.stat.math than in my mailbox, unrequested. Didn't we have a little forwarding talk recently enough for your avarage televangelist viewer to remember it? Kindly cut it out. -j -- "Blah Blah Blah" ___________________________________________________________________ Jamie Lawrence <jamiel(a)sybase.com>

2 1

real live Detweiler messages- anyone?
by jamiel＠sybase.com 17 Dec '03

17 Dec '03

I have a writing assignment coming up (nothing to do with much that is talked about on this list), and some example posts by Detweiler could be really handy as references. I have one from a long time ago, but could use some others- anybody save any for posterity? If so, could you forward me a couple of juicy ones? thanks- -j -- "Blah Blah Blah" ___________________________________________________________________ Jamie Lawrence <jamiel(a)sybase.com>

1 0

cypherpunks mbone this sat.
by Katy Kislitzin 17 Dec '03

17 Dec '03

The monthly meeting of cypherpunks will be broadcast on the mbone this saturday, july 9, from noon to 6. The ports and ID's are as follows: audio: 64264 51421 video: 43185 51863 This is being advertised using sd as well. The main topic will be discussion of Phil Karn's swIPe package, which was announced at Usenix last month. It allows one to set up encrypted IP tunnels over the exisiting IP network. For those of you interested in attending in person, the meeting will be held at Silicon Graphics in Cafe Iris, building 5. To get there, get off 101 at Shorline, and head towards the ampetheater. Turn right into the SGI parking lot after Charlston, and continue down to just before the circle. Building 5 will be on your right. NOTE: THE SGI PARKING LOTS ARE BEING RE-SURFACED THIS WEEKEND. SO YOU WILL HAVE TO PARK ON THE MAIN CAMPUS DRIVE OR ON THE CONCRETE PARKING AREAS. --kt Katy Kislitzin Silicon Graphics I/S Network Software ktk(a)corp.sgi.com

2 1