C DETERMINE PROBABILITY FOR OCCURRENCE OF SELF TRAPPING C PATHS WITH GIVEN LENGTH FOR SELF-AVOIDING RANDOM WALKS C ON THE SQUARE LATTICE C C AUTHOR: HUGO PFOERTNER, HTTP://WWW.ABOUTHUGO.DE/ C C VERSION HISTORY: C 27.11.2002 PRINT POWERS OF 2 AND 3 FOR DENOMINATOR, C USE NORMALIZED FORM OF DENOMINATOR FOR OUTPUT C 22.11.2002 PRINT INFO ABOUT PROGRESS OF COMPUTATION 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 MAXIMUM LENGHT OF 2D WALK PARAMETER ( M=23,MP=M+1 ) DOUBLEPRECISION PACC C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. INTEGER W(-MP:MP,-MP:MP), IP(0:M), JP(0:M), COUN(500) INTEGER*8 IPATH, P, PROB(500), IV, IPROB(36000000), NPOT8 COMMON /FIELD/ W COMMON /WALK/ I, J, L, IP, JP C TO EXCLUDE U-TURNS (REVERSE DIRECTION OF PREVIOUS MOVE) INTEGER EX(4) EXTERNAL NPOT8 DATA EX / 2,1, 4,3 / C CLEAR FIELD DO 10 II = -MP,MP DO 10 JJ = -MP,MP W(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) = 1 W(1,0) = 1 IP(0) = 0 JP(0) = 0 IP(1) = 1 JP(1) = 0 C POSITION AFTER FIRST STEP (ONLY 1 OF THE FOUR POSSIBLE C DIRECTIONS IS SELECTED) I = 1 J = 0 C DO 20 I2 = 1, 4 L = 2 IF ( I2 .EQ. EX(I1) ) GOTO 20 CALL STEP ( I2, *20 ) DO 30 I3 = 1, 4 L = 3 IF ( I3 .EQ. EX(I2) ) GOTO 30 CALL STEP ( I3, *30 ) DO 40 I4 = 1, 4 L = 4 IF ( I4 .EQ. EX(I3) ) GOTO 40 CALL STEP ( I4, *40 ) DO 50 I5 = 1, 4 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, 4 L = 6 IF ( I6 .EQ. EX(I5) ) GOTO 60 CALL STEP ( I6, *60 ) DO 70 I7 = 1, 4 L = 7 IF ( I7 .EQ. EX(I6) ) GOTO 70 CALL STEP ( I7, *70 ) DO 80 I8 = 1, 4 L = 8 IF ( I8 .EQ. EX(I7) ) GOTO 80 CALL STEP ( I8, *80 ) DO 90 I9 = 1, 4 L = 9 IF ( I9 .EQ. EX(I8) ) GOTO 90 CALL STEP ( I9, *90 ) DO 100 I10 = 1, 4 L = 10 IF ( I10 .EQ. EX(I9) ) GOTO 100 CALL STEP ( I10, *100 ) DO 110 I11 = 1, 4 L = 11 IF ( I11 .EQ. EX(I10) ) GOTO 110 CALL STEP ( I11, *110 ) DO 120 I12 = 1, 4 L = 12 IF ( I12 .EQ. EX(I11) ) GOTO 120 CALL STEP ( I12, *120 ) DO 130 I13 = 1, 4 L = 13 IF ( I13 .EQ. EX(I12) ) GOTO 130 CALL STEP ( I13, *130 ) DO 140 I14 = 1, 4 L = 14 IF ( I14 .EQ. EX(I13) ) GOTO 140 CALL STEP ( I14, *140 ) DO 150 I15 = 1, 4 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 C16 DO 160 I16 = 1, 4 C16 L = 16 C16 IF ( I16 .EQ. EX(I15) ) GOTO 160 C16 CALL STEP ( I16, *160 ) C17 DO 170 I17 = 1, 4 C17 L = 17 C17 IF ( I17 .EQ. EX(I16) ) GOTO 170 C17 CALL STEP ( I17, *170 ) C18 DO 180 I18 = 1, 4 C18 L = 18 C18 IF ( I18 .EQ. EX(I17) ) GOTO 180 C18 CALL STEP ( I18, *180 ) C19 DO 190 I19 = 1, 4 C19 L = 19 C19 IF ( I19 .EQ. EX(I18) ) GOTO 190 C19 CALL STEP ( I19, *190 ) C20 DO 200 I20 = 1, 4 C20 L = 20 C20 IF ( I20 .EQ. EX(I19) ) GOTO 200 C20 CALL STEP ( I20, *200 ) C21 DO 210 I21 = 1, 4 C21 L = 21 C21 IF ( I21 .EQ. EX(I20) ) GOTO 210 C21 CALL STEP ( I21, *210 ) C22 DO 220 I22 = 1, 4 C22 L = 22 C22 IF ( I22 .EQ. EX(I21) ) GOTO 220 C22 CALL STEP ( I22, *220 ) C23 DO 230 I23 = 1, 4 C23 L = 23 C23 IF ( I23 .EQ. EX(I22) ) GOTO 230 C23 CALL STEP ( I23, *230 ) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. IPATH = IPATH + 1 IF ( MIN ( W(I-1,J), W(I+1,J), & W(I,J-1), W(I,J+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 C WRITE (*,*) ITRAP, IPATH, IPROB(ITRAP) C DO 300 N = 1,L C WRITE (*,1300) IP(N),JP(N) C1300 FORMAT ( 2 I3 ) C300 CONTINUE ENDIF C C CLEAR POSITION C TO EXTEND TO WALK OF LENGTH NN REMOVE COMMENTS CXX FOR XX<=NN C23 W(IP(23),JP(23)) = 0 C23 I = IP(22) C23 J = JP(22) C23 230 CONTINUE C22 W(IP(22),JP(22)) = 0 C22 I = IP(21) C22 J = JP(21) C22 220 CONTINUE C21 W(IP(21),JP(21)) = 0 C21 I = IP(20) C21 J = JP(20) C21 210 CONTINUE C20 W(IP(20),JP(20)) = 0 C20 I = IP(19) C20 J = JP(19) C20 200 CONTINUE C19 W(IP(19),JP(19)) = 0 C19 I = IP(18) C19 J = JP(18) C19 190 CONTINUE C18 W(IP(18),JP(18)) = 0 C18 I = IP(17) C18 J = JP(17) C18 180 CONTINUE C17 W(IP(17),JP(17)) = 0 C17 I = IP(16) C17 J = JP(16) C17 170 CONTINUE C16 W(IP(16),JP(16)) = 0 C16 I = IP(15) C16 J = JP(15) C16 160 CONTINUE W(IP(15),JP(15)) = 0 I = IP(14) J = JP(14) 150 CONTINUE W(IP(14),JP(14)) = 0 I = IP(13) J = JP(13) 140 CONTINUE W(IP(13),JP(13)) = 0 I = IP(12) J = JP(12) 130 CONTINUE W(IP(12),JP(12)) = 0 I = IP(11) J = JP(11) 120 CONTINUE W(IP(11),JP(11)) = 0 I = IP(10) J = JP(10) 110 CONTINUE W(IP(10),JP(10)) = 0 I = IP(9) J = JP(9) 100 CONTINUE W(IP(9),JP(9)) = 0 I = IP(8) J = JP(8) 90 CONTINUE W(IP(8),JP(8)) = 0 I = IP(7) J = JP(7) 80 CONTINUE W(IP(7),JP(7)) = 0 I = IP(6) J = JP(6) 70 CONTINUE W(IP(6),JP(6)) = 0 I = IP(5) J = JP(5) 60 CONTINUE W(IP(5),JP(5)) = 0 I = IP(4) J = JP(4) 50 CONTINUE W(IP(4),JP(4)) = 0 I = IP(3) J = JP(3) 40 CONTINUE W(IP(3),JP(3)) = 0 I = IP(2) J = JP(2) 30 CONTINUE W(IP(2),JP(2)) = 0 I = IP(1) J = JP(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 A046661) C ITRAP IS THE NUMBER OF TRAPPED PATHS (=2*A077482(N) IN OEIS) WRITE (*,*) 'ITRAP:', ITRAP, ' IPATH:', IPATH IF ( ITRAP .EQ. 0 ) STOP WRITE (*,*) 'PROB = 1/', 1.0D0/PACC C SORT BY INCREASING TRAP PROBABILTIES CSLATEC CALL SSORTI ( IPROB, IPROB, ITRAP, 1) C VSORT8 (N,IA) SORTS AN INTEGER*8 ARRAY IA(N) IN ASCENDING ORDER CALL VSORT8 ( ITRAP, IPROB ) IV = IPROB(1) J = 1 C COUNT PATHS WITH EQUAL TRAP PROBABILTY C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. WRITE (*,1000) 1000 FORMAT (/,' NR 1/PROBABILITY COUNT P3 P2') DO 300 I = 2, ITRAP IF ( IPROB(I) .NE. IV ) THEN IPOT = IPOT + 1 PROB(IPOT) = IV COUN(IPOT) = J WRITE (*,1002) IPOT, IV, J, NPOT8(IV,3), NPOT8(IV,2) 1002 FORMAT ( I4, I15, I8, 2 I4 ) IV = IPROB(I) J = 1 ELSE J = J + 1 ENDIF 300 CONTINUE IPOT = IPOT + 1 PROB(IPOT) = IV COUN(IPOT) = J WRITE (*,1002) IPOT, IV, J, NPOT8(IV,3), NPOT8(IV,2) WRITE (*,*) ' ' CALL PRBSUM ( L, IPOT, PROB, COUN ) C END OF MAIN PROGRAM END C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. SUBROUTINE STEP ( MOVE, * ) PARAMETER ( M=23, MP=M+1) INTEGER INC(4), JNC(4) INTEGER W(-MP:MP,-MP:MP) COMMON /FIELD/ W INTEGER IP(0:M), JP(0:M) COMMON /WALK/ I, J, L, IP, JP DATA INC / 1,-1,0,0 /, JNC / 0,0,1,-1 / C TRY NEXT STEP II = I + INC(MOVE) JJ = J + JNC(MOVE) C CHECK TARGET IF ( W(II,JJ) .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 RETURN 1 ELSE C ADVANCE WALK BY ONE STEP I = II J = JJ C BLOCK VISITED POSITION W(I,J) = 1 C STORE WALK COORDINATES IP(L) = I JP(L) = J C write (11,1000) L, IP(L), JP(L) C1000 FORMAT ( I2, 2I3 ) ENDIF RETURN END C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. SUBROUTINE PRODPR ( P ) C DETERMINE PROBABILITY FOR STORED PATH PARAMETER ( MM=23 ) INTEGER IP(0:MM), JP(0:MM) COMMON /WALK/ I, J, L, IP, JP INTEGER IN(4), JN(4), NF(4) INTEGER*8 P, Q P = 3 DO 10 M = L-1, 2, -1 IN(1) = IP(M)-1 JN(1) = JP(M) IN(2) = IP(M)+1 JN(2) = JP(M) IN(3) = IP(M) JN(3) = JP(M)-1 IN(4) = IP(M) JN(4) = JP(M)+1 NF(1) = 0 NF(2) = 0 NF(3) = 0 NF(4) = 0 DO 20 N = M-1, 0, -1 DO 30 NN = 1, 4 IF ( IP(N) .EQ. IN(NN) .AND. JP(N) .EQ. JN(NN) ) NF(NN) = 1 30 CONTINUE 20 CONTINUE Q = 0 DO 40 NN = 1, 4 IF (NF(NN) .EQ. 0) Q = Q + 1 40 CONTINUE P = P * Q 10 CONTINUE END C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. SUBROUTINE PRBSUM ( LP, IPOT, PROB, COUN ) C C SUM OF PROBABILITIES C C AUTHOR: HUGO PFOERTNER, HTTP://WWW.PFOERTNER.ORG/ C C VERSION HISTORY: C 27.11.2002 OUTPUT WITH NORMALIZED DENOMINATOR C 22.11.2002 16 BYTE REAL FOR NUMERATOR ACCUMULATION C 07.11.2002 INITIAL VERSION C C LP: LENGTH OF PATH C IPOT: NUMBER OF TERMS C COUNT: NUMERATOR C PROB: DENOMINATOR C PRBSUM COMPUTES SUM I=1..IPOT COUN(I)/PROB(I) C THE COMMON DENOMINATOR IS OF THE FORM 3**(LP-1)*2**K C K HAS TO BE SELECTED SUFFICIENTLY LARGE C IMPLICIT INTEGER*8 (A-Z) INTEGER*4 LP, IPOT, I, P2, P3, COUN(IPOT) INTEGER*8 PROB(IPOT) C QUADRUPLE PRECISION REAL AVAILABLE ON VAX COMPUTERS REAL*16 QNUM C C COMMON DENOMINATOR MINIMUM VALUE DEN = 3**(LP-1) C TRY IF DEN IS A MULTIPLE OF PROB DO 10 I = 1, IPOT 20 CONTINUE MUL = DEN / PROB(I) IF ( MUL*PROB(I) .NE. DEN ) THEN DEN = DEN + DEN GOTO 20 ENDIF 10 CONTINUE WRITE (*,*) ' DEN=', DEN C USE QUADRUPLE PRECISION FLOAT TO AVOID OVERFLOW OF PRODUCT COUN(I)*DEN QNUM = 0.0Q0 DO 30 I = 1, IPOT QNUM = QNUM + QFLOAT(COUN(I))*QFLOAT(DEN)/QFLOAT(PROB(I)) 30 CONTINUE C VAX-SPECIFIC EXTENSION OF NINT FUNCTION NUM = KIQNNT ( QNUM ) WRITE (*,*) ' NUM=', NUM C C TRY TO ELIMINATE COMMON FACTORS 2 AND 3 60 CONTINUE IF ( 2*(DEN/2) .EQ. DEN .AND. 2*(NUM/2) .EQ. NUM ) THEN NUM = NUM / 2 DEN = DEN / 2 GOTO 60 ENDIF 70 CONTINUE IF ( 3*(DEN/3) .EQ. DEN .AND. 3*(NUM/3) .EQ. NUM ) THEN NUM = NUM / 3 DEN = DEN / 3 GOTO 70 ENDIF WRITE (*,*) 'NUM/DEN=', NUM, '/', DEN C EXTRACT POWERS OF 2 FROM DENOMINATOR C COUNT EXPONENT OF 2 P2 = 0 80 CONTINUE IF ( 2*(DEN/2) .NE. DEN ) GOTO 90 P2 = P2 + 1 DEN = DEN / 2 GOTO 80 90 CONTINUE C EXTRACT POWERS OF 3 FROM DENOMINATOR C COUNT EXPONENT OF 3 P3 = 0 180 CONTINUE IF ( 3*(DEN/3) .NE. DEN ) GOTO 190 P3 = P3 + 1 DEN = DEN / 3 GOTO 180 190 CONTINUE WRITE (*,*) ' P2=', P2, ' P3=', P3 C ASSUMING THE VALIDITY OF THE CONJECTURE ON THE NUMBER OF C MAXIMALLY 2-CONSTRAINED STWS, C ADJUST OUTPUT TO NORMALIZED FORM (P2S=A076874(LP)-3) P2S = LP - INT ( SQRT (4.0*FLOAT(LP) + 1.0001) ) - 3 P3S = LP - 1 200 CONTINUE IF ( P2S .GT. P2 ) THEN NUM = NUM + NUM P2 = P2 + 1 ELSEIF ( P2S .LT. P2 ) THEN NUM = NUM / 2 P2 = P2 - 1 ELSE GOTO 210 ENDIF GOTO 200 210 CONTINUE 300 CONTINUE IF ( P3S .GT. P3 ) THEN NUM = NUM + NUM + NUM P3 = P3 + 1 ELSEIF ( P3S .LT. P3 ) THEN NUM = NUM / 3 P3 = P3 - 1 ELSE GOTO 310 ENDIF GOTO 300 310 CONTINUE WRITE (*,1003) LP,NUM, P3S, P2S 1003 FORMAT (' P(',I2,') = ', I15, ' / (3^', I2, '*2^', I2 ')' ) RETURN END