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