C DETERMINE PROBABILITY FOR OCCURRENCE OF SELF TRAPPING
C PATHS WITH GIVEN LENGTH FOR SELF-AVOIDING RANDOM WALKS
C ON THE LATTICE Z4
C
C AUTHOR: HUGO PFOERTNER, HTTP://WWW.PFOERTNER.ORG/
C
C VERSION HISTORY:
C 29.11.2002 SPLIT TEST FOR FINAL POSITION
C 28.11.2002 MODIFIED FOR 4D WALKS
C 14.11.2002 PROBABILITY EXACT SUMMATION MOVED TO SEPARATE PROGRAM
C 12.11.2002 MODIFIED FOR THE 3D WALK
C 07.11.2002 SUM PROBABILITIES IN SUBROUTINE PRBSUM
C 05.11.2002 TRY LONGER WALKS, SORT PROBABILITIES
C 26.05.2002 ACCUMULATE BRANCH PROBABILITIES
C 25.05.2002 ADAPTED FOR 2D
C 15.05.2002 SIZE OF FIELD INCREASED TO MP
C 14.05.2002 INITIAL VERSION
C LENGHT OF 2D WALK
      PARAMETER ( M=18,MP=M+1 )
      DOUBLEPRECISION  PACC
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      INTEGER W(-MP:MP,-MP:MP,-MP:MP,-MP:MP),
     &        IP(0:M), JP(0:M), KP(0:M), LP(0:M)
      INTEGER*8 IPROB(1000000), IV, IPATH, P
      INTEGER*8 PROB(1000)
      INTEGER COUN(1000)
C*** INTEGER*8 PROB(100), COUN(100)
      COMMON /FIELD/ W
      COMMON /WALK/ I, J, K, L, N, IP, JP, KP, LP
C TO EXCLUDE U-TURNS (REVERSE DIRECTION OF PREVIOUS MOVE)
      INTEGER EX(8)
      DATA EX / 2,1, 4,3, 6,5, 8,7 /
C LIST OF THE 8 POSSIBLE STEPS
      INTEGER INC(8), JNC(8), KNC(8), LNC(8)
      DATA INC / -1,1,0,0,0,0,0,0 /, JNC / 0,0,-1,1,0,0,0,0 /,
     &     KNC / 0,0,0,0,-1,1,0,0 /, LNC / 0,0,0,0,0,0,-1,1 /
C CLEAR FIELD
      DO 10 II = -MP,MP
      DO 10 JJ = -MP,MP
      DO 10 KK = -MP,MP
      DO 10 LL = -MP,MP
      W(LL,KK,JJ,II) = 0
10    CONTINUE
C COUNTERS
      IPATH = 0
      ITRAP = 0
C ACCUMULATED PROBABILITY
      PACC = 0.0D0
C BLOCK START POINT AND END POINT OF FIRST STEP
      I1 = 1
      W(0,0,0,0) = 1
      W(1,0,0,0) = 1
      IP(0) = 0
      JP(0) = 0
      KP(0) = 0
      LP(0) = 0
      IP(1) = 1
      JP(1) = 0
      KP(1) = 0
      LP(1) = 0
C POSITION AFTER FIRST STEP (ONLY 1 OF THE FOUR POSSIBLE
C DIRECTIONS IS SELECTED)
      I = 1
      J = 0
      K = 0
      L = 0
C
      DO 20 I2 = 2, 8
      N = 2
      I = I + INC(I2)
      J = J + JNC(I2)
      K = K + KNC(I2)
      L = L + LNC(I2)
      W(I,J,K,L) = 1
      IP(N) = I
      JP(N) = J
      KP(N) = K
      LP(N) = L
      DO 30 I3 = 1, 8
      N = 3
      IF ( I3 .EQ. EX(I2) ) GOTO 30
      CALL STEP ( I3, *30 )
      DO 40 I4 = 1, 8
      N = 4
      IF ( I4 .EQ. EX(I3) ) GOTO 40
      CALL STEP ( I4, *40 )
      DO 50 I5 = 1, 8
      N = 5
      IF ( I5 .EQ. EX(I4) ) GOTO 50
      CALL STEP ( I5, *50 )
C PRINT DIAGNOSTIC INFO TO TRACE PROGRESS
      WRITE (*,1001) I2, I3, I4, I5, IPATH, ITRAP
1001  FORMAT ( 4I2, I16, I9 )
      DO 60 I6 = 1, 8
      N = 6
      IF ( I6 .EQ. EX(I5) ) GOTO 60
      CALL STEP ( I6, *60 )
      DO 70 I7 = 1, 8
      N = 7
      IF ( I7 .EQ. EX(I6) ) GOTO 70
      CALL STEP ( I7, *70 )
      DO 80 I8 = 1, 8
      N = 8
      IF ( I8 .EQ. EX(I7) ) GOTO 80
      CALL STEP ( I8, *80 )
      DO 90 I9 = 1, 8
      N = 9
      IF ( I9 .EQ. EX(I8) ) GOTO 90
      CALL STEP ( I9, *90 )
      DO 100 I10 = 1, 8
      N = 10
      IF ( I10 .EQ. EX(I9) ) GOTO 100
      CALL STEP ( I10, *100 )
C      DO 110 I11 = 1, 8
C      N = 11
C      IF ( I11 .EQ. EX(I10) ) GOTO 110
C      CALL STEP ( I11, *110 )
C      DO 120 I12 = 1, 8
C      N = 12
C      IF ( I12 .EQ. EX(I11) ) GOTO 120
C      CALL STEP ( I12, *120 )
C      DO 130 I13 = 1, 8
C      N = 13
C      IF ( I13 .EQ. EX(I12) ) GOTO 130
C      CALL STEP ( I13, *130 )
C      DO 140 I14 = 1, 8
C      N = 14
C      IF ( I14 .EQ. EX(I13) ) GOTO 140
C      CALL STEP ( I14, *140 )
C      DO 150 I15 = 1, 8
C      N = 15
C      IF ( I15 .EQ. EX(I14) ) GOTO 150
C      CALL STEP ( I15, *150 )
C      DO 160 I16 = 1, 8
C      N = 16
C      IF ( I16 .EQ. EX(I15) ) GOTO 160
C      CALL STEP ( I16, *160 )
C      DO 170 I17 = 1, 8
C      N = 17
C      IF ( I17 .EQ. EX(I16) ) GOTO 170
C      CALL STEP ( I17, *170 )
C TO EXTEND TO WALK OF LENGTH NN REMOVE COMMENTS CXX FOR XX<=NN
C      DO 180 I18 = 1, 8
C      N = 18
C      IF ( I18 .EQ. EX(I17) ) GOTO 180
C      CALL STEP ( I18, *180 )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      IPATH = IPATH + 1
      IF ( W(I-1,J,K,L) .EQ. 0 ) GOTO 400
      IF ( W(I+1,J,K,L) .EQ. 0 ) GOTO 400
      IF ( W(I,J-1,K,L) .EQ. 0 ) GOTO 400
      IF ( W(I,J+1,K,L) .EQ. 0 ) GOTO 400
      IF ( W(I,J,K-1,L) .EQ. 0 ) GOTO 400
      IF ( W(I,J,K+1,L) .EQ. 0 ) GOTO 400
      IF ( W(I,J,K,L-1) .EQ. 0 ) GOTO 400
      IF ( W(I,J,K,L+1) .EQ. 0 ) GOTO 400
C      IF ( MIN ( W(I-1,J,K,L), W(I+1,J,K,L),
C     &           W(I,J-1,K,L), W(I,J+1,K,L),
C     &           W(I,J,K-1,L), W(I,J,K+1,L),
C     &           W(I,J,K,L-1), W(I,J,K,L+1) ) .EQ. 1 )  THEN
          ITRAP = ITRAP + 1
C DETERMINE FACTOR FOR PROBABILITY PRODUCT
          CALL PRODPR ( P )
          PACC = PACC + 1.0D0 / DBLE(P)
          IPROB(ITRAP) = P
C TO PRINT ALL TRAPPED WALKS, ACTIVATE UNTIL LABEL 300
          WRITE (*,*) ITRAP, IPATH, IPROB(ITRAP)
          WRITE (11,*) ITRAP, IPATH, IPROB(ITRAP)
          DO 300 N1 = 1,N
          WRITE (11,1300) IP(N1),JP(N1),KP(N1),LP(N1)
1300      FORMAT ( 4 I3 )
300       CONTINUE
C      ENDIF
C
C CLEAR POSITION
C TO EXTEND TO WALK OF LENGTH NN REMOVE COMMENTS CXX FOR XX<=NN
C      W(IP(18),JP(18),KP(18),LP(18)) = 0
C      I = IP(17)
C      J = JP(17)
C      K = KP(17)
C      L = LP(17)
C 180  CONTINUE
C      W(IP(17),JP(17),KP(17),LP(17)) = 0
C      I = IP(16)
C      J = JP(16)
C      K = KP(16)
C      L = LP(16)
C 170  CONTINUE
C      W(IP(16),JP(16),KP(16),LP(16)) = 0
C      I = IP(15)
C      J = JP(15)
C      K = KP(15)
C      L = LP(15)
C 160  CONTINUE
C      W(IP(15),JP(15),KP(15),LP(15)) = 0
C      I = IP(14)
C      J = JP(14)
C      K = KP(14)
C      L = LP(14)
C 150   CONTINUE
C      W(IP(14),JP(14),KP(14),LP(14)) = 0
C      I = IP(13)
C      J = JP(13)
C      K = KP(13)
C      L = LP(13)
C140   CONTINUE
C      W(IP(13),JP(13),KP(13),LP(13)) = 0
C      I = IP(12)
C      J = JP(12)
C      K = KP(12)
C      L = LP(12)
C130   CONTINUE
C      W(IP(12),JP(12),KP(12),LP(12)) = 0
C      I = IP(11)
C      J = JP(11)
C      K = KP(11)
C      L = LP(11)
C120   CONTINUE
C      W(IP(11),JP(11),KP(11),LP(11)) = 0
C      I = IP(10)
C      J = JP(10)
C      K = KP(10)
C      L = LP(10)
C110   CONTINUE
      W(IP(10),JP(10),KP(10),LP(10)) = 0
      I = IP(9)
      J = JP(9)
      K = KP(9)
      L = LP(9)
100   CONTINUE
      W(IP(9),JP(9),KP(9),LP(9)) = 0
      I = IP(8)
      J = JP(8)
      K = KP(8)
      L = LP(8)
90    CONTINUE
      W(IP(8),JP(8),KP(8),LP(8)) = 0
      I = IP(7)
      J = JP(7)
      K = KP(7)
      L = LP(7)
80    CONTINUE
      W(IP(7),JP(7),KP(7),LP(7)) = 0
      I = IP(6)
      J = JP(6)
      K = KP(6)
      L = LP(6)
70    CONTINUE
      W(IP(6),JP(6),KP(6),LP(6)) = 0
      I = IP(5)
      J = JP(5)
      K = KP(5)
      L = LP(5)
60    CONTINUE
      W(IP(5),JP(5),KP(5),LP(5)) = 0
      I = IP(4)
      J = JP(4)
      K = KP(4)
      L = LP(4)
50    CONTINUE
      W(IP(4),JP(4),KP(4),LP(4)) = 0
      I = IP(3)
      J = JP(3)
      K = KP(3)
      L = LP(3)
40    CONTINUE
      W(IP(3),JP(3),KP(3),LP(3)) = 0
      I = IP(2)
      J = JP(2)
      K = KP(2)
      L = LP(2)
30    CONTINUE
      W(IP(2),JP(2),KP(2),LP(2)) = 0
      I = IP(1)
      J = JP(1)
      K = KP(1)
      L = LP(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 A010575(N)/8)
C ITRAP IS THE NUMBER OF TRAPPED PATHS (WILL BE SUBMITTED TO OEIS)
      WRITE (*,*) 'ITRAP:', ITRAP, ' IPATH:', IPATH
      IF ( ITRAP .EQ. 0 ) STOP
      WRITE (*,*) 'PROB = 1/', 1.0D0/PACC
C SORT BY INCREASING DENOMINATORS OF TRAP PROBABILITIES
CSLATEC      CALL SSORTI ( IPROB, IPROB, ITRAP, 1)
C VSORT8 (N,A) SORTS AN I*8 ARRAY A(N) IN ASCENDING ORDER
      CALL VSORT8 ( ITRAP, IPROB )
      IV = IPROB(1)
      J = 1
C COUNT PATHS HAVING THE SAME TRAP PROBABILITY
      DO 400 I = 2, ITRAP
      IF ( IPROB(I) .NE. IV ) THEN
        IPOT = IPOT + 1
        PROB(IPOT) = IV
        COUN(IPOT) = J
        WRITE (*,*) IV, J
        WRITE (11,*) IV, J
        IV = IPROB(I)
        J = 1
      ELSE
        J = J + 1
      ENDIF
400   CONTINUE
      IPOT = IPOT + 1
      PROB(IPOT) = IV
      COUN(IPOT) = J
      WRITE (*,*) IV, J
      WRITE (11,*) IV, J
C END OF MAIN PROGRAM
      END
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      SUBROUTINE STEP ( MOVE, * )
      PARAMETER ( M=18, MP=M+1)
      INTEGER INC(8), JNC(8), KNC(8), LNC(8)
      INTEGER W(-MP:MP,-MP:MP,-MP:MP,-MP:MP)
      COMMON /FIELD/ W
      INTEGER IP(0:M), JP(0:M), KP(0:M), LP(0:M)
      COMMON /WALK/ I, J, K, L, N, IP, JP, KP, LP
      DATA INC / -1,1,0,0,0,0,0,0 /, JNC / 0,0,-1,1,0,0,0,0 /,
     &     KNC / 0,0,0,0,-1,1,0,0 /, LNC / 0,0,0,0,0,0,-1,1 /
C TRY NEXT STEP
      II = I + INC(MOVE)
      JJ = J + JNC(MOVE)
      KK = K + KNC(MOVE)
      LL = L + LNC(MOVE)
C CHECK TARGET
      IF ( W(II,JJ,KK,LL) .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) -N, II, JJ, KK, LL
        RETURN 1
      ELSE
C ADVANCE WALK BY ONE STEP
        I = II
        J = JJ
        K = KK
        L = LL
C BLOCK VISITED POSITION
        W(I,J,K,L) = 1
C STORE WALK COORDINATES
        IP(N) = I
        JP(N) = J
        KP(N) = K
        LP(N) = L
C        write (11,1000) N, IP(N), JP(N), KP(N), LP(N)
C1000     FORMAT ( I3, 4I3 )
      ENDIF
      RETURN
      END
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      SUBROUTINE PRODPR ( P )
C DETERMINE PROBABILITY FOR STORED PATH
      PARAMETER ( MM=18 )
      INTEGER IP(0:MM), JP(0:MM), KP(0:MM), LP(0:MM)
      COMMON /WALK/ I, J, K, L, N, IP, JP, KP, LP
      INTEGER IN(8), JN(8), KN(8), LN(8), NF(8)
      INTEGER*8 P, Q
      P = 7
      DO 10 M = N-1, 2, -1
      IN(1) = IP(M)-1
      JN(1) = JP(M)
      KN(1) = KP(M)
      LN(1) = LP(M)
      IN(2) = IP(M)+1
      JN(2) = JP(M)
      KN(2) = KP(M)
      LN(2) = LP(M)
      IN(3) = IP(M)
      JN(3) = JP(M)-1
      KN(3) = KP(M)
      LN(3) = LP(M)
      IN(4) = IP(M)
      JN(4) = JP(M)+1
      KN(4) = KP(M)
      LN(4) = LP(M)
      IN(5) = IP(M)
      JN(5) = JP(M)
      KN(5) = KP(M)-1
      LN(5) = LP(M)
      IN(6) = IP(M)
      JN(6) = JP(M)
      KN(6) = KP(M)+1
      LN(6) = LP(M)
      IN(7) = IP(M)
      JN(7) = JP(M)
      KN(7) = KP(M)
      LN(7) = LP(M)-1
      IN(8) = IP(M)
      JN(8) = JP(M)
      KN(8) = KP(M)
      LN(8) = LP(M)+1
      NF(1) = 0
      NF(2) = 0
      NF(3) = 0
      NF(4) = 0
      NF(5) = 0
      NF(6) = 0
      NF(7) = 0
      NF(8) = 0
      DO 20 N1 = M-1, 0, -1
      DO 30 NN = 1, 8
      IF ( IP(N1) .EQ. IN(NN) .AND.
     &     JP(N1) .EQ. JN(NN) .AND.
     &     KP(N1) .EQ. KN(NN) .AND.
     &     LP(N1) .EQ. LN(NN) ) NF(NN) = 1
30    CONTINUE
20    CONTINUE
      Q = 0
      DO 40 NN = 1, 8
      IF (NF(NN) .EQ. 0) Q = Q + 1
40    CONTINUE
      P = P * Q
10    CONTINUE
      END
