program PNfiniteBuneven
! Identification of leading pairwise-nonrecurrant sequences, by Thomas Bewley.
! Note: compile with, eg: gfortran PNfiniteBuneven.f90 -o PNfiniteBuneven
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! user input follows !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
logical :: periodic=.true.  ! Specify the periodic or nonperiodic case.
integer :: mm(2)=(/ 6, 2 /), r(2)=(/ 8, 7 /), show=400
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  end user input !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
integer :: m,n,num(2)=0,location(9)=0,pairs(65)=0,a(9)=0,w(9),answer(65)=0,shown=0,best(65)=0,b,e,t
integer*8 :: group(2000000,2)=0
! Set values for rl, rs (# of times each symbol is repeated) and n (# of total slots)
n=mm(1)*r(1)+mm(2)*r(2); m=mm(1)+mm(2); best(1:n)=m
! Now find those groupings of r locations in the n available slots in which a given
! symbol can be repeated without violating the properties of pairwise nonrecurrance.
write(*,'(A,L2,4(A,I2))') 'periodic=',periodic,', m=',m,', r(1)=',r(1),', r(2)=',r(2),', n=',n
call system('date');
write(*,'(A,3(I2,A))') 'Looking for pairwise nonrecurrent groupings of length n=',n,' with r(1)=',r(1),' and r(2)=',r(2)
do k=1,n; pairs(1:n-1)=100; location(1)=k; call FindGroups(2,r(1),group(:,1),num(1)); end do;
write(*,'(A,I9,3(A,I2))') 'Found ',num(1),' pairwise nonrecurrent groupings with r(1)=',r(1),' (we will use',mm(1),' of these)'
do k=1,n; pairs(1:n-1)=100; location(1)=k; call FindGroups(2,r(2),group(:,2),num(2)); end do;
write(*,'(A,I9,3(A,I2))') 'Found ',num(2),' pairwise nonrecurrent groupings with r(2)=',r(2),' (we will use',mm(2),' of these)'
call system('date');

write(*,'(A,I9,A)') 'Now looking at the possible ways to fit such groupings together'
w(1:mm(1))=1; w(mm(1)+1:m)=2
l7: do; iter=iter+1
  write(*,'(A,9(I2))') 'Now checking with w=',w(1:m)

  l6: do loop=1,num(w(1)); a(1)=loop
    if (.not. BTEST(group(a(1),w(1)),n-1)) then
      call system('date'); write(6,*) 'Exhaustively checked all possible fits with this ordering.'; exit l6
    end if; call Check(2)
  end do l6

  ! Knuth's Algorithm to generate all possible permutations of w (that is, which groups are long, and which are short).
  b=m-1; do while (b>0 .and. w(b)>=w(b+1)); b=b-1; end do;  ! b=begin, e=end
  e=m; do while (w(b)>=w(e)); e=e-1; end do
  if (b==0) exit                         ! NOTE: If exhausted all permutations, quit!
  t=w(b); w(b)=w(e); w(e)=t; b=b+1; e=m; ! (this is just the second half of Algorithm L)
  do while (b<e); t=w(b); w(b)=w(e); w(e)=t; b=b+1; e=e-1; end do
end do l7
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
recursive subroutine FindGroups(level,r,group,num)
  integer :: level,s1,s2,i,j,r
  integer*8 :: group(:)
  logical :: success
  do i=location(level-1)+1,n-(r-level)
  	do j=1,n-1; if (pairs(j)>=level) pairs(j)=100; end do
	location(level)=i; success=.true.
	l1: do j=level-1,1,-1
	  s1=i-location(j); s2=location(j)-i+n
	  if (periodic) then
        if (pairs(s1)>level .and. pairs(s2)>level) then; pairs(s1)=level; pairs(s2)=level
        else; success=.false.; exit l1; end if
      else
        if (pairs(s1)>level) then; pairs(s1)=level
		else; success=.false.; exit l1; end if
      end if
	end do l1
	if (success) then
	  if (level==r) then;
		  num=num+1; do j=1,r; group(num)=IBSET(group(num),n-location(j)); end do
!		  write(6,*) num,location,group(num); call FLUSH(6)
	  else; call FindGroups(level+1,r,group,num); end if
	end if
  end do
end subroutine FindGroups
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
recursive subroutine Check(level)
logical   :: test
integer   :: pos,x,x1,x2
integer*8 :: temp,temp1
temp=group(a(1),w(1)); do j=2,level-1; temp=ior(group(a(j),w(j)),temp); end do
l4: do pos=n-1,0,-1; if (.not. btest(temp,pos)) exit l4; end do l4
temp=0; temp=IBSET(temp,pos);    ! note: pos is the most significant bit of the next group.
temp1=0; temp1=IBSET(temp1,pos+1);  
x1=1; x2=num(w(level)) ! Use bisection algorithm to find starting point in search for next
do while (x2-x1>1)              ! group (find first entry such that pos+1 is the most significant bit)
   x=(x2+x1)/2; if (temp1>group(x,w(level)) .and. temp1 > 0) then; x2=x; else; x1=x; end if
end do                          ! need to check sign above only for n=64, integer*8 case

l3: do i=x1,num(w(level)); a(level)=i; test=.true.
  l2: do j=level-1,1,-1
    test=(test .and. (iand(group(a(j),w(j)),group(a(level),w(level)))==0)); if (.not. test) exit l2
  end do l2                        ! need to check sign below only for n=64, integer*8 case
  if (temp >= group(a(level),w(level)) .and. group(a(level),w(level)) > 0) exit l3  
  if (test) then
    if (level<m) then
      call Check(level+1)
    else
      do l=1,m; do k=1,n; if (BTEST(group(a(l),w(l)),k-1)) answer(n+1-k)=l; end do; end do
      write(6,'(100I2)') answer(1:n)-1; call FLUSH(6); shown=shown+1;
	  l5: do k=1,n;
	    if (best(k)>answer(k)) then
	      write(6,*) 'Above entry is lexographically first of all sequences found thus far'
		  best=answer; exit l5
		else if (best(k)<answer(k)) then
		  exit l5
		end if
	  end do l5
	  if (shown==show) then; call system('date'); stop; end if
    end if
  end if
end do l3
end subroutine check
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end program PNfiniteBuneven