cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        Program to find Golomb rulers by exhaustive backtrack search
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                  IBM SOFTWARE DISCLAIMER
c
c   grs2.f (version 1.0)
c   Copyright (1999,1986)
c   International Business Machines Corporation
c
c   Permission to use, copy, modify and distribute this software for
c any purpose and without fee is hereby granted, provided that this
c copyright and permission notice appear on all copies of the software.
c The name of the IBM corporation may not be used in any advertising or
c publicity pertaining to the use of the software.  IBM makes no
c warranty or representations about the suitability of the software
c for any purpose.  It is provided "AS IS" without any express or
c implied warranty, including the implied warranties of merchantability,
c fitness for a particular purpose and non-infringement.  IBM shall not
c be liable for any direct, indirect, special or consequential damages
c resulting from the loss of use, data or projects, whether in action
c of contract or tort, arising out or in the connection with the use or
c performance of this software.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        Author: James B. Shearer
c        email:  jbs@watson.ibm.com
c        website: http://www.research.ibm.com/people/s/shearer/
c        date: 1999 (partially derived from code written in 1986)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        This program constructs all Golomb rulers with n marks and
c   length k (or less) by exhaustive backtrack search.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c   This is a lightly tuned version of grs1.f.
c
c   Loop 60 was rewritten to eliminate the if test.  This included
c reversing the meaning of 0 and 1 entries in the id array as this
c makes the new form of loop 60 more efficient.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      parameter(maxn=50,maxk=1000,nprint=10)
      integer*4 ip(maxn),ie(maxn),m(maxn),mub(maxn),mubv(maxn)
      integer*4 id(maxk),l(maxn*maxk)
      real*8 t1,t2,time,tit
c get number of marks and length
   90 read(5,100)n,k
  100 format(2i5)
c n=0 to end program
      if(n.eq.0)go to 300
c start timing
      call cputime(t1,iret)
c initialize differences used and mark positions available
      do 10 j=1,k
      id(j)=1
      l(j)=j
   10 continue
c initialize mark position bound arrays
      do 15 j=1,n
      mub(j)=k
      mubv(j)=k
   15 continue
c use mub to force middle mark to left of middle position
      do 16 j=1,(n+1)/2
      jb=(n+1)/2 - j
      mub(j)=(k-1)/2 + mod(n,2) - 1 - jb*(jb+1)/2
   16 continue
c flag used for even number of marks
      iflag=mod(n+1,2)*n/2
c initialize backtrack search
      ns=0
      tit=0.d0
      m(1)=0
      ip(2)=0
      ie(2)=k
      nc=2
c start backtrack search
      go to 40
c backtrack portion of search loop  
c update search depth
   20 nc=nc-1
c release differences
   25 do 30 j=1,nc-1
      id(m(nc)-m(j))=1
   30 continue
c check for end of search
      if(nc.eq.1)go to 80
c find next node in search tree, increment node count
   40 ip(nc)=ip(nc)+1
      tit=tit+1.d0
c tests to prune search
c are there enough valid mark positions remaining to complete ruler?
      if(ie(nc)-ip(nc).lt.n-nc)go to 20
c add mark
      m(nc)=l(ip(nc))
c check mark position against fixed and variable bounds
      if(m(nc).gt.mub(nc))go to 20
      if(m(nc).gt.mubv(nc))go to 20
c do we have a complete ruler
      if(n.eq.nc)go to 70
c mark new differences used
      do 50 j=1,nc-1
      id(m(nc)-m(j))=0
   50 continue
c check if enough small differences remain to complete the ruler
      j=0
      jc=nc
      i=m(nc)
c find next unused difference
   55 j=j+1
      if(id(j).eq.0)go to 55
c add differences to partial ruler length
      i=i+j
      jc=jc+1
c n marks yet?
      if(jc.lt.n)go to 55
c can ruler be completed?
c note failure goes to 25 not 20 as we have not bumped nc yet
c we have to undo 50 loop and we may not be done at this level
      if(i.gt.k)go to 25
c set mubv bound assuming first difference is largest remaining
      mubv(nc+1)=m(nc)+k-i+j
c sum of middle marks must be less than length
      if(iflag.eq.nc)mub(nc+1)=k-m(nc)-1
c update list of eligible mark positions
      ip(nc+1)=ie(nc)
      i=ie(nc)
c1    do 60 j=ip(nc)+1,ie(nc)
c1    if(id(l(j)-m(nc)).eq.0)go to 60
c1    i=i+1
c1    l(i)=l(j)
c1 60 continue
      do 60 j=ip(nc)+1,ie(nc)
      lt=l(j)
      l(i+1)=lt
      i=i+id(lt-m(nc))
   60 continue
      ie(nc+1)=i
c increase depth and continue search
      nc=nc+1
      go to 40
c solution found, check if it satisfies middle mark property
   70 if(m((n+1)/2)+m((n+2)/2).gt.m(n))go to 40
c increment number of solutions found
      ns=ns+1
c check if we should output solution
      if(ns.gt.nprint)go to 40
      write(1,200)n,m(n)
      write(1,210)(m(j),j=1,n)
      write(6,220)(m(j),j=1,n)
  200 format(2i10)
  210 format(10i6)
  220 format(1x,10i6)
      go to 40
c search complete
   80 continue
c stop timing
      call cputime(t2,irc)
c output search statistics
      write(6,230)n,k,ns,t2-t1,tit
  230 format(1x,'n=',i3,' k=',i4,' solutions=',i8,' time=',-6pf15.6,
     x ' nodes=',0pf15.0)
c go get next input
      go to 90
c program done, mark end of output file and stop
  300 continue
      write(1,200)0,0
      stop
      end
c timing subroutine for aix
c cputime is vs fortran timing routine
c t is microseconds
      subroutine cputime(t,irc)
      real*8 t
      t=mclock()*10000.d0
c for g77 on pcs replace above statement with
c     t=second()*1000000.d0
      return
      end