C DETERMINE DISTRIBUTION OF TERMINATION POINTS OF C PATHS WITH GIVEN NUMBER OF STEPS FOR SELF-AVOIDING RANDOM WALKS C ON THE SQUARE LATTICE C C AUTHOR: HUGO PFOERTNER, HTTP://WWW.PFOERTNER.ORG/ C C VERSION HISTORY: C 05.02.2002 OUTPUT OF SUM OF SQUARES C 04.12.2002 AVERAGES DOUBLE PRECISION, RUN UP TO N=24 C 03.12.2002 INITIAL VERSION C MAXIMUM NUMBER OF STEPS OF 2D WALK PARAMETER ( M=25,MP=M+1,M2=M*M ) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. INTEGER W(-MP:MP,-MP:MP), IP(0:M), JP(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, L, IP, JP C TO EXCLUDE U-TURNS (REVERSE DIRECTION OF PREVIOUS MOVE) INTEGER EX(4) 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 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) = 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 DO 160 I16 = 1, 4 L = 16 IF ( I16 .EQ. EX(I15) ) GOTO 160 CALL STEP ( I16, *160 ) DO 170 I17 = 1, 4 L = 17 IF ( I17 .EQ. EX(I16) ) GOTO 170 CALL STEP ( I17, *170 ) DO 180 I18 = 1, 4 L = 18 IF ( I18 .EQ. EX(I17) ) GOTO 180 CALL STEP ( I18, *180 ) DO 190 I19 = 1, 4 L = 19 IF ( I19 .EQ. EX(I18) ) GOTO 190 CALL STEP ( I19, *190 ) DO 200 I20 = 1, 4 L = 20 IF ( I20 .EQ. EX(I19) ) GOTO 200 CALL STEP ( I20, *200 ) DO 210 I21 = 1, 4 L = 21 IF ( I21 .EQ. EX(I20) ) GOTO 210 CALL STEP ( I21, *210 ) DO 220 I22 = 1, 4 L = 22 IF ( I22 .EQ. EX(I21) ) GOTO 220 CALL STEP ( I22, *220 ) DO 230 I23 = 1, 4 L = 23 IF ( I23 .EQ. EX(I22) ) GOTO 230 CALL STEP ( I23, *230 ) DO 240 I24 = 1, 4 L = 24 IF ( I24 .EQ. EX(I23) ) GOTO 240 CALL STEP ( I24, *240 ) DO 250 I25 = 1, 4 L = 25 IF ( I25 .EQ. EX(I24) ) GOTO 250 CALL STEP ( I25, *250 ) C...+....1....+....2....+....3....+....4....+....5....+....6....+....7.. IPATH = IPATH + 1 MAND = ABS(I) + ABS(J) MANALL(MAND) = MANALL(MAND) + 1 NEUC = I*I + J*J EUDALL(NEUC) = EUDALL(NEUC) + 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 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 W(IP(25),JP(25)) = 0 I = IP(24) J = JP(24) 250 CONTINUE W(IP(24),JP(24)) = 0 I = IP(23) J = JP(23) 240 CONTINUE W(IP(23),JP(23)) = 0 I = IP(22) J = JP(22) 230 CONTINUE W(IP(22),JP(22)) = 0 I = IP(21) J = JP(21) 220 CONTINUE W(IP(21),JP(21)) = 0 I = IP(20) J = JP(20) 210 CONTINUE W(IP(20),JP(20)) = 0 I = IP(19) J = JP(19) 200 CONTINUE W(IP(19),JP(19)) = 0 I = IP(18) J = JP(18) 190 CONTINUE W(IP(18),JP(18)) = 0 I = IP(17) J = JP(17) 180 CONTINUE W(IP(17),JP(17)) = 0 I = IP(16) J = JP(16) 170 CONTINUE W(IP(16),JP(16)) = 0 I = IP(15) J = JP(15) 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 (*,*) '#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=25, 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