program PNfinite
! Identification of leading pairwise-nonrecurrant sequences, by Thomas Bewley.
logical :: periodic=.true.  ! Change this flag for the periodic or nonperiodic case.
integer :: nmin=2, nmax=7   ! Change this to specify the values of m to test
! Note: compile with, eg:   g95 PNfinite.f90 -o PNfinite   Run in background with, eg:
! date > periodic.dat; PNfinite >> periodic.dat; date >> periodic.dat &
integer :: a(50),count(7),pairs(7,50),first(7)
integer :: b,e,n,m,flag,s,s1,s2,i,is,k,k2,ks,kb,ks1,t; integer*8 :: iter,ip=100000000
do m=nmin,nmax; iter=0; 
  if (periodic) then; n=m*m; select case (m)
    case (2); a(1:n)=[0,0,1,1]+1;
    case (3); a(1:n)=[0,0,0,1,1,1,2,2,2]+1;
    case (4); a(1:n)=[0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3]+1; 
    case (5); a(1:n)=[0,0,0,0,0,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4]+1; 
    case (6); a(1:n)=[0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,              &
                      4,4,4,4,4,4,5,5,5,5,5,5]+1;
    case (7); a(1:n)=[0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,3,3,      &
                      4,4,4,4,4,4,4,5,5,5,5,5,5,5,6,6,6,6,6,6,6]+1;
  end select; else; select case (m)
    case (2); n=7;  a(1:n)=[0,0,0,0,1,1,1]+1;
    case (3); n=14; a(1:n)=[0,0,0,0,0,1,1,1,1,1,2,2,2,2]+1;
    case (4); n=24; a(1:n)=[0,0,0,0,0,0,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3]+1; 
    case (5); n=35; a(1:n)=[0,0,0,0,0,0,0,1,1,1,1,1,1,1,2,2,2,2,2,2,2,              &
                            3,3,3,3,3,3,3,4,4,4,4,4,4,4]+1; 
    case (6); n=48; a(1:n)=[0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,        &
                            3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5]+1;
  end select; end if
  write(*,'(/A,L1,A,I2,A,I2/)') 'periodicity= ',periodic,', m=',m,', n=',n
  do; iter=iter+1; flag=1;
    ! Sweep to detect & disregard large groups of sequences with repeated pairs
    pairs(1:m,1:n-1)=0;    
    l1: do ks=2,n; do k=1,ks-1; if (a(k)==a(ks)) then; s1=ks-k; s2=k-ks+n
      if (periodic) then
        if (pairs(a(k),s1)==0 .and. s1/=s2) then; pairs(a(k),s1)=k; pairs(a(k),s2)=k
        else; flag=0; exit l1; end if
      else
        if (pairs(a(k),s1)==0) then; pairs(a(k),s1)=k; else; flag=0; exit l1; end if
      end if
    end if; end do; end do l1	
    if (flag==1) exit  ! NOTE: If successful, quit!
    ! Check to see if the first occurances of each symbol appear in order; if not,
    ! adjust ks appropriately, thus again disregarding large groups of invalid sequences
    first(1:m)=0;
    do t=1,m; do k=1,n; if (a(k)==t) then; first(t)=k; exit; end if; end do; end do
    do t=2,m; do s=1,t-1; if (first(t)<first(s)) ks=min(ks,first(t)-1); end do; end do
    ! Print progress on screen
    if (mod(iter,ip)==1) then        
      write(*,'(A,ES9.2,A,I2,A,7I3)') 'iter=',real(iter),', ks=',ks,', first =',first(1:m)
      write(*,'(A,72I2)') 'current: a=',a(1:n)-1; call FLUSH(6);
    end if
    ! Then, reorder elements after ks from largest to smallest, and increment a 
    ! (thus incrementing the ks element) via Alg L below
    count(1:m)=0; do i=ks+1,n; count(a(i))=count(a(i))+1; end do
    t=ks+1; do i=m,1,-1; l=count(i); a(t:t+l-1)=i; t=t+l; end do
    ! Finally, find next permutation (Algorithm L in Knuth 2005, Sec 7.2.1.2)
    b=min(ks,n-1); do while (b>0 .and. a(b)>=a(b+1)); b=b-1; end do;  ! b=begin, e=end
    e=n; do while (a(b)>=a(e)); e=e-1; end do
    if (b==0) exit                         ! NOTE: If exhausted all permutations, quit!
    t=a(b); a(b)=a(e); a(e)=t; b=b+1; e=n; ! (this is just the second half of Algorithm L)
    do while (b<e); t=a(b); a(b)=a(e); a(e)=t; b=b+1; e=e-1; end do
  end do
  write(*,'(A,72I2)') 'done: a=',a(1:n)-1
  write(*,'(A,I2,A,I2,A,I2,A,I32/)') 'success=',flag,', m=',m,', n=',n,', iter=',iter
end do
end program PNfinite