From: Jim choate <ravage@bga.com>
Date: Thu, 7 Jul 94 13:11:14 PDT
To: cypherpunks@toad.com
Subject: (fwd) Random Numbers - SPECTRAL.FOR
Message-ID: <199407072011.PAA29181@ivy.bga.com>
MIME-Version: 1.0
Content-Type: text/plain


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