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@netcom.com Subject: Random Numbers - SPECTRAL.FOR Message-ID: <deleydCsIBEG.L2B@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