C Count geometrically connected components when 1,2,...,n positions
C are selected in an n*n square array. Two positions (i,j),(k,l)
C are adjacent if max(abs(i-k),abs(j-l))<=1
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Change history:
C Sep 14 2004 Type DOUBLEPRECISION for occurrence counts 
C Sep 13 2004 modified call to gradfs from GRAFPACK
C Sep 12 2004 Initial version
C
      implicit integer (a-z)
C The parameters m and n have to be adjusted before compilation
C to select an mXm square array an n<=m selectable positions
      parameter ( m=9, m2=m*m, n=7 )
C S: Square array where the selected numbers are marked
C F: S in 1-d arangement
C T: The list of selected numbers
C R: Row positions of selected numbers
C C: Column positions of selected numbers
      dimension S(m,m), F(m2), T(n), R(m2), C(m2)
      doubleprecision TICK, LOOK, COMCNT(n)
C The following equivalence statement is responsible for the
C mapping 1-d lottery numbers 2-d arrangement on ticket.
      equivalence (S,F)
C Comment/uncomment n lines of the following equivalence statement
      equivalence (t(1),z1)
     &           ,(t(2),z2)
     &           ,(t(3),z3)
     &           ,(t(4),z4)
     &           ,(t(5),z5)
     &           ,(t(6),z6)
     &           ,(t(7),z7)
c     &           ,(t(8),z8)
c     &           ,(t(9),z9)
C A: The adjacency matrix
      dimension A(m,m)
C Output of the graph debth first search program
      dimension dad(m), order(m)
C
C Clear counters for numbers of components
      do 4 i = 1, n
4     comcnt(i) = 0.0D0
C Clear square array
      do 5 I = 1, m2
5     f(i) = 0
C Start of loops over all possible choices
C100   continue
C
      tick = 0.0D0
      LOOK = 0.0D0
C
C Comment / uncomment the required number of blocks of 4 lines
C dependent on the value of n
C
      do 110 z1 = 1, m2-n+1
      f(z1) = 1
      c(1) = 1 + mod(t(1)-1,m)
      r(1) = 1 + (t(1)-1)/m
C
      do 120 z2 = z1+1, m2-n+2
      f(z2) = 2
      c(2) = 1 + mod(t(2)-1,m)
      r(2) = 1 + (t(2)-1)/m
C
      do 130 z3 = z2+1, m2-n+3
      f(z3) = 3
      c(3) = 1 + mod(t(3)-1,m)
      r(3) = 1 + (t(3)-1)/m
C
      do 140 z4 = z3+1, m2-n+4
      f(z4) = 4
      c(4) = 1 + mod(t(4)-1,m)
      r(4) = 1 + (t(4)-1)/m
C
      do 150 z5 = z4+1, m2-n+5
      f(z5) = 5
      c(5) = 1 + mod(t(5)-1,m)
      r(5) = 1 + (t(5)-1)/m
C
      do 160 z6 = z5+1, m2-n+6
      f(z6) = 6
      c(6) = 1 + mod(t(6)-1,m)
      r(6) = 1 + (t(6)-1)/m
C
      do 170 z7 = z6+1, m2-n+7
      f(z7) = 7
      c(7) = 1 + mod(t(7)-1,m)
      r(7) = 1 + (t(7)-1)/m
C
c      do 180 z8 = z7+1, m2-n+8
c      f(z7) = 8
c      c(8) = 1 + mod(t(8)-1,m)
c      r(8) = 1 + (t(8)-1)/m
C
c      do 190 z9 = z8+1, m2-n+9
c      f(z7) = 9
c      c(9) = 1 + mod(t(9)-1,m)
c      r(9) = 1 + (t(9)-1)/m
C
C
      tick = tick + 1.0D0
C
C Determine adjacency matrix. First clear
      do 20 i = 1, n
      do 30 j = 1, m
      a(j,i) = 0
30    continue
C
C Check adjacent positions with boundary limitations
      do 40 j = max(1,r(i)-1),min(m,r(i)+1)
      do 50 k = max(1,c(i)-1),min(m,c(i)+1)
      if ( s(k,j) .ne. 0 ) then
        a(s(k,j),i) = 1
      endif
50    continue
40    continue
20    continue
C
C Check corresponding graph using GRAPH_ADJ_DFS
C http://www.csit.fsu.edu/~burkardt/f_src/grafpack/grafpack.html
C
      call gradfs ( a, m, n, dad, order, ncomp )
C
C Count occurrence of ncomp components in the adjacency graph
C
      comcnt(ncomp) = comcnt(ncomp) + 1.0D0
C
C Progress indicator
      if ( tick .GE. look ) then
        write (*,1001) tick, comcnt
        look = look + 1.0D7
1001    format (F14.0,/, 5F14.0,:,/, 5F14.0 )
      endif
C Clear single entry of ticket
c      f(t(9)) = 0
c190   continue
c      f(t(8)) = 0
c180   continue
      f(t(7)) = 0
170   continue
      f(t(6)) = 0
160   continue
      f(t(5)) = 0
150   continue
      f(t(4)) = 0
140   continue
      f(t(3)) = 0
130   continue
      f(t(2)) = 0
120   continue
      f(t(1)) = 0
110   continue
C
C Print result
      write (*,*) m,n
      write (*,1001) tick,  comcnt
      write (*,*) ' Enter 0 to stop'
      read (*,*) istop
      if ( istop .eq. 0 ) stop
      end