c     =======================================================================
c     PROGRAM: 2-d random walk test.
c     -----------------------------------------------------------------------
c     for    a two-dimensional random walk test for testing
c            of pseudorandom numbers
c     by     I. Vattulainen, vattulai@convex.csc.fi
c     alg    random walks, a chi-square test
c     ref    Phys. Rev. Lett. 73, 2513 (1994).
c     title  2drwtest.f
c     size   11kB
c     prec   single/double
c     lang   Fortran77
c     -----------------------------------------------------------------------
c     The author of this software is I. Vattulainen.  Copyright (c) 1994.
c     Permission to use, copy, modify, and distribute this software for
c     any purpose without fee is hereby granted, provided that this entire
c     notice is included in all copies of any software which is or includes
c     a copy or modification of this software and in all copies of the
c     supporting documentation for such software.
c     THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
c     WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKE ANY
c     REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
c     OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
c     -----------------------------------------------------------------------
c     In the two-dimensional random walk test, we consider a sequence of
c     pseudorandom numbers, which determine the directions of jumps of a
c     random walker on a plane. In the test, we calculate the distribution
c     of the final position of a random walker after N steps, and if it
c     deviates (significantly) from the expected behavior, the pseudorandom
c     number generator fails the test. The original reference for this test
c     is I. Vattulainen et al., Phys. Rev. Lett. 73, 2513 (1994).
c     -----------------------------------------------------------------------
c     In this long version of the test program, we calculate three different
c     quantities of interest:
c       (1) the distribution of the random walker's final position to end
c           up in each of four equivalent blocks (separated by axes x = 0
c           and y = 0) on the plane,
c       (2) the mean square displacement as a function of walk length, and
c       (3) the spatial distribution function P(x,y,N) to find the random
c           walker on lattice site (x,y) after N steps.
c     If quantities in (1) and (2) are not needed, parameters IDIST, N2 and
c     ILATT (as well as other lines related to them) may be commented out.
c     -----------------------------------------------------------------------
c     Main parameters in the test are as follows:
c        N          Number of steps in a random walk.
c        ISCANS     Number of samples (random walks).
c        IDIST      Mean square displacement as a function of walk length.
c        ILATT      Distribution to visit lattice site (x,y) after N steps.
c     ------------------------------------------------------
c     Required modification to this code:
c                   initialization and implementation of the pseudorandom
c                   number generator (PRNG), which will be tested. For
c                   initialization, there is a reserved space in the code
c                   below. The PRNG may be appended to the end of the file,
c                   or may be called separately during compiling. In this
c                   version, the default generator is GGL.
c
c     Output files: fort.20    Quantity (1)
c                   fort.21    Quantity (2)
c                   fort.22    Quantity (3)
c     =======================================================================

      IMPLICIT NONE
      INTEGER
     +     N, N2, ISCANS
      PARAMETER( N = 500, N2 = N*N )
      PARAMETER( ISCANS = 10000 )

      INTEGER
     +     I, J, K, ICNT, IX, IY,
     +     IDIST(N2), ILATT(-N:N,-N:N)
      DOUBLE PRECISION
     +     DP1, DP2, DP3
      REAL*8
     +     DSEED, DISEED
      REAL
     +     CHISUM, REXPEC, RBLOCK(4), RAN(N), GGL

c     Initialize some variables:

      DP1 = 0.25D0
      DP2 = 0.5D0
      DP3 = 0.75D0
      DO 25 I=1,4
         RBLOCK(I) = 0.
 25   CONTINUE
      DO 27 I=1,N2
         IDIST(I) = 0
 27   CONTINUE
      DO 29 I=-N,N
         DO 28 J=-N,N
            ILATT(I,J) = 0
 28      CONTINUE
 29   CONTINUE
      CHISUM = 0.
      REXPEC = FLOAT(ISCANS)/4

c     ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
c     First of all, initialize the pseudorandom number generator.
c     In the case of GGL, only one number is needed:

      DSEED   = 12345.D0
c     ```````````````````````````````````````````````````````````
      DISEED = DSEED

c     Begin the sampling loop:

      DO 100 ICNT=1,ISCANS

c     (((((((((((((((((((((((((((((((((((((((((((((((((((((((((((
c     Generate a sequence of N U(0,1] distributed random numbers.
c     If U[0,1), U[0,1], or U(0,1) distributed random numbers are
c     used, the bias (systematic error) is most likely negligible.

         DO 110 I=1,N
            RAN(I) = GGL(DSEED)
 110     CONTINUE

c     )))))))))))))))))))))))))))))))))))))))))))))))))))))))))))

c     Perform a random walk:

         CALL RWALK(N,N2,IX,IY,DP1,DP2,DP3,RAN,IDIST,ILATT)

c     Calculate the distribution for the final position of
c     a random walker:

         IF(IX.LT.0) THEN
            IF(IY.LT.0) THEN
               RBLOCK(1) = RBLOCK(1) + 1
            ELSE
               IF(IY.GT.0) THEN
                  RBLOCK(2) = RBLOCK(2) + 1
               ELSE
                  IF(IY.EQ.0) THEN
                     RBLOCK(1) = RBLOCK(1) + 0.5
                     RBLOCK(2) = RBLOCK(2) + 0.5
                  ENDIF
               ENDIF
            ENDIF
         ELSE
            IF(IX.GT.0) THEN
               IF(IY.LT.0) THEN
                  RBLOCK(3) = RBLOCK(3) + 1
               ELSE
                  IF(IY.GT.0) THEN
                     RBLOCK(4) = RBLOCK(4) + 1
                  ELSE
                     IF(IY.EQ.0) THEN
                        RBLOCK(3) = RBLOCK(3) + 0.5
                        RBLOCK(4) = RBLOCK(4) + 0.5
                     ENDIF
                  ENDIF
               ENDIF
            ELSE
               IF(IY.LT.0) THEN
                  RBLOCK(1) = RBLOCK(1) + 0.5
                  RBLOCK(3) = RBLOCK(3) + 0.5
               ELSE
                  IF(IY.GT.0) THEN
                     RBLOCK(2) = RBLOCK(2) + 0.5
                     RBLOCK(4) = RBLOCK(4) + 0.5
                  ELSE
                     RBLOCK(1) = RBLOCK(1) + 0.25
                     RBLOCK(2) = RBLOCK(2) + 0.25
                     RBLOCK(3) = RBLOCK(3) + 0.25
                     RBLOCK(4) = RBLOCK(4) + 0.25
                  ENDIF
               ENDIF
            ENDIF
         ENDIF

 100  CONTINUE

c     _______________________________________________
c     After the test, calculate the chi-square value:
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

      DO 111 I=1,4
         CHISUM = CHISUM +
     +        ((RBLOCK(I) - REXPEC)**2)/REXPEC
 111  CONTINUE

c     And finally write the main results:

      WRITE(20,*) ' ==============================================='
      WRITE(20,*) ' Parameters in the random walk test:'
      WRITE(20,*) '      Initial seed:                 ', NINT(DISEED)
      WRITE(20,*) '      Walk length (N):              ', N
      WRITE(20,*) '      Number of samples (ISCANS):   ', ISCANS
      WRITE(20,*)
      WRITE(20,*) ' Results:'
      WRITE(20,*) '      Probabilities to end up at each block:'
      WRITE(20,*) '               Block 1: ', RBLOCK(1)/ISCANS
      WRITE(20,*) '               Block 2: ', RBLOCK(2)/ISCANS
      WRITE(20,*) '               Block 3: ', RBLOCK(3)/ISCANS
      WRITE(20,*) '               Block 4: ', RBLOCK(4)/ISCANS
      WRITE(20,*) '      Chi-square value:   ', CHISUM
      WRITE(20,*)
      WRITE(20,*) ' To compare with, with three degree of freedom'
      WRITE(20,*) ' some percentage points of chi-square values'
      WRITE(20,*) ' are as follows:'
      WRITE(20,*) '       1% limit:   0.1148'
      WRITE(20,*) '       5% limit:   0.3518'
      WRITE(20,*) '      95% limit:   7.815'
      WRITE(20,*) '      99% limit:  11.345'
      WRITE(20,*) ' ==============================================='
      CLOSE(20)

      DO 112 I=1,N
         WRITE(21,*) DFLOAT(IDIST(I))/ISCANS
 112  CONTINUE
      CLOSE(21)

      DO 114 J=-N,N
         DO 113 I=-N,N
            IF(ILATT(I,J).NE.0) THEN
               WRITE(22,*) I, J, ILATT(I,J)
            ENDIF
 113     CONTINUE
 114  CONTINUE
      CLOSE(22)

      STOP
      END

c     ======================================================
c     Generate a random walk on a plane. The number of steps
c     is N and the random numbers given in parameter RAN are
c     the main input to this subroutine.
c     ======================================================

      SUBROUTINE RWALK(N,N2,IX,IY,DP1,DP2,DP3,RAN,IDIST,ILATT)

      IMPLICIT NONE
      INTEGER
     +     N, N2, IX, IY, I, IDIST(N2), ID, ILATT(-N:N,-N:N)
      DOUBLE PRECISION
     +     DP1, DP2, DP3
      REAL
     +     RANX, RAN(N)

      IX = 0
      IY = 0

      DO 200 I=1,N
         RANX = RAN(I)
         IF(RANX.GT.DP3) THEN
            IX = IX + 1
         ELSE
            IF(RANX.GT.DP2) THEN
               IX = IX - 1
            ELSE
               IF(RANX.GT.DP1) THEN
                  IY = IY + 1
               ELSE
                  IY = IY - 1
               ENDIF
            ENDIF
         ENDIF
         ID = IX**2 + IY**2
         IDIST(I) = IDIST(I) + ID
 200  CONTINUE
      ILATT(IX,IY) = ILATT(IX,IY) + 1

      RETURN
      END

C     ************************************************************
C     ============================================================
C     A pseudorandom number generator GGL
C     ------------------------------------------------------------

      REAL FUNCTION GGL (DS)

      DOUBLE PRECISION
c     +     D1,
     +     DS, D2
c      DATA D1/2147483648.D0/
      DATA D2/2147483647.D0/
      DS  = DMOD(16807.D0*DS,D2)
c     Generate U(0,1] distributed random numbers:
      GGL = DS/D2

      RETURN
      END

C     ============================================================