C Two-colorings of an nXn grid such that no 4 points
C of any color form a square
C Additionally: Count numbers of different sizes of squares
C
C Author: Hugo Pfoertner http://www.pfoertner.org/
C
C Change history:
C June 6 2005 Diagnostic output for 1st square of new size
C June 5 2005 Initial version
C
      implicit integer (a-z)
      parameter (n=6,m=n*n, ld=2*m)
      integer*8 c, limit
      dimension xx(n,n), yy(n,n), x(m), y(m), d(6), dia(ld),
     &           p(n,n), q1(10000), q2(10000), q3(10000), q4(10000)
      logical ss(n,n), s(m)
      equivalence (x,xx),(y,yy),
     &      (d(1),d12), (d(2),d23), (d(3),d31),
     &      (d(4),d14), (d(5),d24), (d(6),d34),
     &      (s,ss)
C
      pos = -1
      do 5 ix = 1, n
      do 5 iy = 1, n
      pos = pos + 1
      xx(ix,iy) = ix
      yy(ix,iy) = iy
      p(ix,iy) = pos
5     continue
C preset occurrence counts of diagonal lengths
      do 10 i = 1, ld
10    dia(i) = 0
      nsq = 0
C Output file
      open ( unit=11, file='sqout.txt', form='formatted',
     &        status='unknown' )
C Loop over all combinations of 4 points
      do 100 k1 = 1, m-3
      x1 = x(k1)
      y1 = y(k1)
      do 200 k2 = k1+1, m-2
      x2 = x(k2)
      y2 = y(k2)
      d12 = (x2-x1)**2 + (y2-y1)**2
      do 300 k3 = k2+1, m-1
      x3 = x(k3)
      y3 = y(k3)
      d23 = (x2-x3)**2 + (y2-y3)**2
      d31 = (x1-x3)**2 + (y1-y3)**2
      do 400 k4 = k3+1, m
      x4 = x(k4)
      y4 = y(k4)
      d14 = (x1-x4)**2 + (y1-y4)**2
      d24 = (x2-x4)**2 + (y2-y4)**2
      d34 = (x3-x4)**2 + (y3-y4)**2
      dmi = min ( d12, d23, d31, d14, d24, d34 ) 
      dma = max ( d12, d23, d31, d14, d24, d34 )
C Check
C Pythagoras
      if ( dmi + dmi .ne. dma ) goto 400
      nmi = 0
      nma = 0
C 4 sides and 2 diagonals
      do 500 i = 1, 6
      if ( d(i) .eq. dmi ) nmi = nmi + 1
      if ( d(i) .eq. dma ) nma = nma + 1
500   continue
C      write (*,*) nmi, nma
      if ( nmi .ne. 4 ) goto 400
      if ( nma .ne. 2 ) goto 400
      nsq = nsq + 1
C Diagnostic output
C      write (*,1000) nsq,
C     &  x1, y1,  x2, y2,  x3, y3,  x4, y4
C Code of square
      Q1(nsq) = p(x1,y1)
      Q2(nsq) = p(x2,y2)
      Q3(nsq) = p(x3,y3)
      Q4(nsq) = p(x4,y4)
C count occurrences of diagonal lengths
      dia(dma) = dia(dma) + 1 
C Write first example of a particular size to output
      if ( dia(dma) .eq. 1 ) then
        write (11,*) ' Diag^2=', dma
        write (11,1000) nsq,
     &  x1, y1,  x2, y2,  x3, y3,  x4, y4
1000    format ( i5, 4 ( 2x, 2I2 ) )
      endif
400   continue
300   continue
200   continue
100   continue
      write (*,*) ' n=', n, ' Squares:', nsq
      write (11,*) ' n=', n, ' Squares:', nsq
      ndial = 0
      do 600 i = 1, ld
      if ( dia(i) .gt. 0 ) then
        ndial = ndial + 1
        write (11,*) i, dia(i)
      endif
600   continue
      write (*,*) ' ndial:', ndial
      write (11,*) ' ndial:', ndial
      close ( unit=11 )
      write (*,*) ' Color=1, stop=0:'
      read (*,*) icolor
      if ( icolor .ne. 1 ) stop
C Generate all colorings
      nosqua = 0
      mm = m
      if ( mm .lt. 63 ) then
        limit = 2**mm - 1
      else
        limit = 2**62 + (2**62-1)
      endif
C Loop over all colorings (written as 64bit binary)
      do 1100 c = 0, limit
C Progress indicator
      if ( mod(c,2**30) .eq. 0) write (*,*) c
C Test all squares
      do 1200 i = 1, nsq
      bit = 0
      if ( btest(c,q1(i)) ) bit = bit + 1
      if ( btest(c,q2(i)) ) bit = bit + 1
      if ( btest(c,q3(i)) ) bit = bit + 1
      if ( btest(c,q4(i)) ) bit = bit + 1
      if ( bit .eq. 0  .or.  bit .eq. 4 ) goto 1210
1200  continue
C No square found, write binary equivalent to file 10
      nosqua = nosqua + 1
      write (*,*) nosqua, c 
      write (10,*) nosqua, c 
1210  continue
1100  continue
C write number of square avoiding configs
      write (*,*) ' nosqua:', nosqua
C To inhibit closing of text window at program stop
      write (*,*) ' Stop:?'
      read (*,*) nsq
      end
