C DETERMINE DISTRIBUTION OF TERMINATION POINTS OF
C PATHS WITH GIVEN NUMBER OF STEPS FOR SELF-AVOIDING RANDOM WALKS
C ON THE CUBIC LATTICE Z3
C
C AUTHOR: HUGO PFOERTNER, HTTP://WWW.PFOERTNER.ORG/
C
C VERSION HISTORY:
C 09.12.2002 COMMENT IMPROVED
C 05.12.2002 OUTPUT OF SUM OF SQUARES, INTEGER*8 AVERAGES
C 04.12.2002 AVERAGES DOUBLE PRECISION
C 03.12.2002 INITIAL VERSION
C MAXIMUM NUMBER OF STEPS OF 2D WALK
      PARAMETER ( M=19,MP=M+1,M2=3*M*M )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      INTEGER W(-MP:MP,-MP:MP,-MP:MP), IP(0:M), JP(0:M), KP(0:M)
      INTEGER*8 IPATH, P
      INTEGER*8 SSASUM, MAASUM, SSTSUM, MATSUM
      REAL*16 EUASUM, EUTSUM
      INTEGER*8 EUDALL(M2), MANALL(M2), EUDTRP(M2), MANTRP(M2)
      COMMON /FIELD/ W
      COMMON /WALK/ I, J, K, L, IP, JP, KP
C TO EXCLUDE U-TURNS (REVERSE DIRECTION OF PREVIOUS MOVE)
      INTEGER EX(6)
      DATA EX / 2,1, 4,3, 6,5 /
C CLEAR FIELD
      DO 10 II = -MP,MP
      DO 10 JJ = -MP,MP
      DO 10 KK = -MP,MP
      W(KK,JJ,II) = 0
10    CONTINUE
C RESET COUNTERS
      DO 11 II = 1, M2
      MANALL(II) = 0
      MANTRP(II) = 0
11    CONTINUE
      DO 12 II = 1, M2
      EUDALL(II) = 0
      EUDTRP(II) = 0
12    CONTINUE
C COUNTERS
      IPATH = 0
      ITRAP = 0
C BLOCK START POINT AND END POINT OF FIRST STEP
      I1 = 1
      W(0,0,0) = 1
      W(1,0,0) = 1
      IP(0) = 0
      JP(0) = 0
      KP(0) = 0
      IP(1) = 1
      JP(1) = 0
      KP(1) = 0
C POSITION AFTER FIRST STEP (ONLY 1 OF THE SIX POSSIBLE
C DIRECTIONS IS SELECTED)
      I = 1
      J = 0
      K = 0
C
      DO 20 I2 = 1, 6
      L = 2
      IF ( I2 .EQ. EX(I1) ) GOTO 20
      CALL STEP ( I2, *20 )
      DO 30 I3 = 1, 6
      L = 3
      IF ( I3 .EQ. EX(I2) ) GOTO 30
      CALL STEP ( I3, *30 )
      DO 40 I4 = 1, 6
      L = 4
      IF ( I4 .EQ. EX(I3) ) GOTO 40
      CALL STEP ( I4, *40 )
      DO 50 I5 = 1, 6
      L = 5
      IF ( I5 .EQ. EX(I4) ) GOTO 50
      CALL STEP ( I5, *50 )
C PRINT DIAGNOSTIC INFO TO TRACE PROGRESS
C      WRITE (*,1001) I2, I3, I4, I5, IPATH, ITRAP
C1001  FORMAT ( 4I2, I16, I12 ) 
      DO 60 I6 = 1, 6
      L = 6
      IF ( I6 .EQ. EX(I5) ) GOTO 60
      CALL STEP ( I6, *60 )
      DO 70 I7 = 1, 6
      L = 7
      IF ( I7 .EQ. EX(I6) ) GOTO 70
      CALL STEP ( I7, *70 )
      DO 80 I8 = 1, 6
      L = 8
      IF ( I8 .EQ. EX(I7) ) GOTO 80
      CALL STEP ( I8, *80 )
      DO 90 I9 = 1, 6
      L = 9
      IF ( I9 .EQ. EX(I8) ) GOTO 90
      CALL STEP ( I9, *90 )
      DO 100 I10 = 1, 6
      L = 10
      IF ( I10 .EQ. EX(I9) ) GOTO 100
      CALL STEP ( I10, *100 )
      DO 110 I11 = 1, 6
      L = 11
      IF ( I11 .EQ. EX(I10) ) GOTO 110
      CALL STEP ( I11, *110 )
      DO 120 I12 = 1, 6
      L = 12
      IF ( I12 .EQ. EX(I11) ) GOTO 120
      CALL STEP ( I12, *120 )
      DO 130 I13 = 1, 6
      L = 13
      IF ( I13 .EQ. EX(I12) ) GOTO 130
      CALL STEP ( I13, *130 )
      DO 140 I14 = 1, 6
      L = 14
      IF ( I14 .EQ. EX(I13) ) GOTO 140
      CALL STEP ( I14, *140 )
      DO 150 I15 = 1, 6
      L = 15
      IF ( I15 .EQ. EX(I14) ) GOTO 150
      CALL STEP ( I15, *150 )
C TO EXTEND TO WALK OF LENGTH NN REMOVE COMMENTS CXX FOR XX<=NN
      DO 160 I16 = 1, 6
      L = 16
      IF ( I16 .EQ. EX(I15) ) GOTO 160
      CALL STEP ( I16, *160 )
      DO 170 I17 = 1, 6
      L = 17
      IF ( I17 .EQ. EX(I16) ) GOTO 170
      CALL STEP ( I17, *170 )
C18      DO 180 I18 = 1, 6
C18      L = 18
C18      IF ( I18 .EQ. EX(I17) ) GOTO 180
C18      CALL STEP ( I18, *180 )
C19      DO 190 I19 = 1, 6
C19      L = 19
C19      IF ( I19 .EQ. EX(I18) ) GOTO 190
C19      CALL STEP ( I19, *190 )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      IPATH = IPATH + 1
      MAND = ABS(I) + ABS(J) + ABS(K)
      MANALL(MAND) = MANALL(MAND) + 1
      NEUC = I*I + J*J + K*K
      EUDALL(NEUC) = EUDALL(NEUC) + 1
      IF ( MIN ( W(I-1,J,K), W(I+1,J,K),
     &           W(I,J-1,K), W(I,J+1,K),
     &           W(I,J,K-1), W(I,J,K+1) ) .EQ. 1 )  THEN
          ITRAP = ITRAP + 1
          MANTRP(MAND) = MANTRP(MAND) + 1
          EUDTRP(NEUC) = EUDTRP(NEUC) + 1
      ENDIF
C
C CLEAR POSITION
C TO EXTEND TO WALK OF LENGTH NN REMOVE COMMENTS CXX FOR XX<=NN
C19      W(IP(19),JP(19),KP(19)) = 0
C19      I = IP(18)
C19      J = JP(18)
C10      K = KP(18)
C19 190  CONTINUE
C18      W(IP(18),JP(18),KP(18)) = 0
C18      I = IP(17)
C18      J = JP(17)
C18      K = KP(17)
C18 180  CONTINUE
      W(IP(17),JP(17),KP(17)) = 0
      I = IP(16)
      J = JP(16)
      K = KP(16)
 170  CONTINUE
      W(IP(16),JP(16),KP(16)) = 0
      I = IP(15)
      J = JP(15)
      K = KP(15)
 160  CONTINUE
      W(IP(15),JP(15),KP(15)) = 0
      I = IP(14)
      J = JP(14)
      K = KP(14)
150   CONTINUE
      W(IP(14),JP(14),KP(14)) = 0
      I = IP(13)
      J = JP(13)
      K = KP(13)
140   CONTINUE
      W(IP(13),JP(13),KP(13)) = 0
      I = IP(12)
      J = JP(12)
      K = KP(12)
130   CONTINUE
      W(IP(12),JP(12),KP(12)) = 0
      I = IP(11)
      J = JP(11)
      K = KP(11)
120   CONTINUE
      W(IP(11),JP(11),KP(11)) = 0
      I = IP(10)
      J = JP(10)
      K = KP(10)
110   CONTINUE
      W(IP(10),JP(10),KP(10)) = 0
      I = IP(9)
      J = JP(9)
      K = KP(9)
100   CONTINUE
      W(IP(9),JP(9),KP(9)) = 0
      I = IP(8)
      J = JP(8)
      K = KP(8)
90    CONTINUE
      W(IP(8),JP(8),KP(8)) = 0
      I = IP(7)
      J = JP(7)
      K = KP(7)
80    CONTINUE
      W(IP(7),JP(7),KP(7)) = 0
      I = IP(6)
      J = JP(6)
      K = KP(6)
70    CONTINUE
      W(IP(6),JP(6),KP(6)) = 0
      I = IP(5)
      J = JP(5)
      K = KP(5)
60    CONTINUE
      W(IP(5),JP(5),KP(5)) = 0
      I = IP(4)
      J = JP(4)
      K = KP(4)
50    CONTINUE
      W(IP(4),JP(4),KP(4)) = 0
      I = IP(3)
      J = JP(3)
      K = KP(3)
40    CONTINUE
      W(IP(3),JP(3),KP(3)) = 0
      I = IP(2)
      J = JP(2)
      K = KP(2)
30    CONTINUE
      W(IP(2),JP(2),KP(2)) = 0
      I = IP(1)
      J = JP(1)
      K = KP(1)
20    CONTINUE
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      WRITE (*,*) '#LENGTH OF WALK:', L
C
C IPATH IS THE NUMBER OF ALL POSSIBLE PATHS OF GIVEN LENGTH
C (OEIS A001412/6)
C ITRAP IS THE NUMBER OF TRAPPED PATHS (=2*A077817(N) IN OEIS)
      WRITE (*,*)  '#IPATH:', IPATH, ' ITRAP:', ITRAP
      SSASUM = 0
      MAASUM = 0
      SSTSUM = 0
      MATSUM = 0
      EUASUM = 0.0Q0
      EUTSUM = 0.0Q0
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      WRITE (*,*)  '#DIST    SAWMAN    SAWEUR2   STWMAN  STWEUR2'
      DO 400 I = 1, M2
      IF ( MANALL(I) + EUDALL(I) .GT. 0 )
     &  WRITE (*,1400) I, MANALL(I), EUDALL(I), MANTRP(I), EUDTRP(I)
C SUM OF ALL WALK LENGTHS EITHER MANHATTAN OR EUCLID
        SSASUM = SSASUM + I*EUDALL(I)
        EUASUM = EUASUM + QSQRT(QFLOAT(I)) * EUDALL(I)
        MAASUM = MAASUM + I*MANALL(I)
        SSTSUM = SSTSUM + I*EUDTRP(I)
        EUTSUM = EUTSUM + QSQRT(QFLOAT(I)) * EUDTRP(I)
        MATSUM = MATSUM + I*MANTRP(I)
400   CONTINUE
1400  FORMAT ( I4, 2I11, 2I9 )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      WRITE (*,1401) DBLE(MAASUM)/DBLE(IPATH), EUASUM/IPATH,
     & DBLE(MATSUM)/DBLE(MAX(1,ITRAP)), EUTSUM/MAX(1,ITRAP)
1401  FORMAT ('#Averages Manhattan all, SumSquares all, same trapped:',
     &        /,'#', 4F13.8 )
C
C TO BE SUBMITTED TO OEIS
      WRITE (*,*) '#ManhSums: ', MAASUM, MATSUM
      WRITE (*,*) '#SquarSums:', SSASUM, SSTSUM
C END OF MAIN PROGRAM
      END
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      SUBROUTINE STEP ( MOVE, * )
      PARAMETER ( M=19, MP=M+1)
      INTEGER INC(6), JNC(6), KNC(6)
      INTEGER W(-MP:MP,-MP:MP,-MP:MP)
      COMMON /FIELD/ W
      INTEGER IP(0:M), JP(0:M), KP(0:M)
      COMMON /WALK/ I, J, K, L, IP, JP, KP
      DATA INC / 1,-1,0,0,0,0 /, JNC / 0,0,1,-1,0,0 /
     &     KNC / 0,0,0,0,1,-1 /
C TRY NEXT STEP
      II = I + INC(MOVE)
      JJ = J + JNC(MOVE)
      KK = K + KNC(MOVE)
C CHECK TARGET
      IF ( W(II,JJ,KK) .NE. 0 ) THEN
C TARGET ALREADY VISITED BEFORE, USE ALTERNATE RETURN TO
C SIGNAL IMPOSSIBLE MOVE DIRECTION TO THE CALLING PROGRAM
C        WRITE (11,1000) -L, II, JJ, KK
        RETURN 1
      ELSE
C ADVANCE WALK BY ONE STEP
        I = II
        J = JJ
        K = KK
C BLOCK VISITED POSITION
        W(I,J,K) = 1
C STORE WALK COORDINATES
        IP(L) = I
        JP(L) = J
        KP(L) = K
C        write (11,1000) L, IP(L), JP(L), KP(L)
C1000     FORMAT ( I2, 3I3 )
      ENDIF
      RETURN
      END
