cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        Program to find difference triangle sets by exhaustive search
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c                  IBM SOFTWARE DISCLAIMER
c
c   dts2.f (version 1.0)
c   Copyright (2000,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: 2000 (partially derived from code written in 1986)
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        This program constructs all difference triangle sets with nt
c   triangles, n marks per triangle and length k (or less) by
c   exhaustive backtrack search.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        This program is based on the Golomb ruler finding programs.  It
c   is most similar to the lightly tuned grs2.f
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      parameter(maxn=50,maxk=1000,nprint=10,nout=1)
c   maxn is maximun allowed number of marks
c   maxk is maximum allowed length
c   nprint is maximum number of solutions to output
c   nout output progress reports when at level nout, =1 for final only
      integer*4 ip(maxn),ie(maxn),m(maxn),mub(maxn),mubv(maxn)
      integer*4 iflag(maxn),ngt(maxn),kgt(maxn)
      integer*4 id(maxk),l(maxn*maxk)
      real*8 t1,t2,time,tit
c get number of marks per ruler, length and number of rulers
c input not checked for validity
   90 read(5,100)n,k,nt
  100 format(3i5)
c n=0 to end program
      if(n.eq.0)go to 300
c start timing
      call cputime(t1,iret)
c compute total number of marks, must not exceed maxn
      nm=n*nt
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, iflag and stuff for loop 55
      do 15 j=1,nm
      mub(j)=k
      mubv(j)=k
      iflag(j)=0
c number of rulers not yet started at depth j
      ntl=nt-1-(j-1)/n
c number of marks placed in current ruler at depth j
      ja=mod(j-1,n)+1
c number of first differences to come at depth j
      ngt(j)=n-ja+ntl*(n-1)
c maximum possible sum of first differences in this and remaining rulers
c at depth j
      kgt(j)=k+ntl*k-(ntl*(ntl+1))/2
   15 continue
c use mub to force middle marks to left of middle position
      do 16 ja=1,nt
c flag used for even number of marks
      if(mod(n,2).eq.0)iflag((ja-1)*n+n/2)=1
      do 16 j=1,(n+1)/2
      jb=(n+1)/2 - j
      mub((ja-1)*n+j)=(k-1)/2 + mod(n,2) - 1 - (jb*(jb+1))/2
   16 continue
c middle positions
      na=n/2
      nb=(n-1)/2
c initialize backtrack search
      ns=0
      tit=0.d0
      m(1)=0
c first differences decrease so must be at least nt in first ruler
      ip(2)=nt-1
      ie(2)=k
      nr=1
      nc=2
      nl=1
      nh=n
c start backtrack search
      go to 40
c backtrack portion of search loop  
c update search depth
   20 nc=nc-1
c release differences, start loop at nl for current ruler only
   25 do 30 j=nl,nc-1
      id(m(nc)-m(j))=1
   30 continue
c check for end of search for current ruler
      if(nc.eq.nl)go to 80
c if nc = nout go output progress report
      if(nc.eq.nout)go to 185
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.nh-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 mark new differences used, start loop at nl for current ruler only
      do 50 j=nl,nc-1
      id(m(nc)-m(j))=0
   50 continue
c do we have a complete ruler?  note this check now after loop 50
c since we may add more rulers.  this means we go to 25 instead of
c 40 when we back up since we have to undo loop 50
      if(nc.eq.nh)go to 70
c check if enough small differences remain to complete the ruler
      j=0
      jc=0
      i=m(nc)
      ng=ngt(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 ng differences yet?
      if(jc.lt.ng)go to 55
c can all remaining rulers 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.kgt(nc))go to 25
c set mubv bound assuming first difference is largest remaining
      mubv(nc+1)=m(nc)+kgt(nc)-i+j
c sum of middle marks must be less than length
c note use l(ie(nc)) which may be less than k
      if(iflag(nc).eq.1)mub(nc+1)=l(ie(nc))-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 current ruler complete, check if it satisfies middle mark property
   70 if(m(nc-na)+m(nc-nb).gt.m(nc))go to 25
c have all rulers been completed
      if(nc.eq.nm)go to 170
c start next ruler
c ntl rulers remain
      ntl=nt-nr
c find ntl smallest unused differences
      j=0
      jc=0
   72 j=j+1
      if(id(j).eq.0)go to 72
      jc=jc+1
      if(jc.lt.ntl)go to 72
c since we are assuming first differences decrease
c j is the mimimum possible for first difference in current ruler
c this must be less than first difference in previous ruler
      if(j.gt.m(nc-n+2))go to 25
c initialize eligible mark position list for current ruler
      ja=j
      l(ie(nc)+1)=0
      do 75 j=ja,k
   75 l(ie(nc)+2+j-ja)=j
      ip(nc+1)=ie(nc)
      ie(nc+1)=ie(nc)+2+k-ja
c upper bound for first difference in current ruler
      mub(nc+2)=m(nc-n+2)
c find sum ntl largest unused differences
      j=k+1
      jc=0
      jsum=0
   76 j=j-1
      if(id(j).eq.0)go to 76
      jsum=jsum+j
      jc=jc+1
      if(jc.lt.ntl)go to 76
c update kgt array
      do 77 j=1,n
      kgt(nr*n+j)=jsum
   77 continue
c find largest unused differences
      j=k+1
   78 j=j-1
      if(id(j).eq.0)go to 78
c update mub array
      ja=(j-1)/2+iand(n,1)-1
      do 79 j=3,(n+1)/2
      jb=(n+1)/2 - j
      mub(nr*n+j)=ja - (jb*(jb+1))/2
   79 continue
c update search status information
      nr=nr+1
      nc=nc+1
      nl=nl+n
      nh=nh+n
c go resume search
      go to 40
c increment number of solutions found
  170 ns=ns+1
c check if we should output solution
      if(ns.gt.nprint)go to 25
      write(1,200)n,m(n),nt
      write(1,210)(m(j),j=1,nm)
      write(6,220)(m(j),j=1,nm)
  200 format(3i10)
  210 format(10i6)
  220 format(1x,10i6)
      go to 25
c search complete at current ruler level
   80 continue
c is this first ruler?  if so we are done
      if(nl.eq.1)go to 180
c else backup to previous ruler and resume search
      nr=nr-1
      nc=nc-1
      nl=nl-n
      nh=nh-n
      go to 25
c stop timing
  180 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=',-6pf16.2,
     x ' nodes=',0pf15.0)
c go get next input
      go to 90
c progress report wanted
  185 call cputime(t2,irc)
c output current search statistics
      write(6,235)n,k,ns,t2-t1,tit,m(nc)
  235 format(1x,'n=',i3,' k=',i4,' solutions=',i4,' time= ',-6pf16.2,
     x ' nodes=',0pf15.0,i5)
c go resume search
      go to 40
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