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