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
