c------------------------------------------------------
c
c              Suduko solver 
c
c          "move" algorithm from setn0()[finds unique moves]
c           and  setn() [finds smallest multiple choices if no unique]
c
c        setn0:   take a digit and see how many places it could fit
c                  in each sub-block; row; col.  select if unique else:
c        setn:    enumerate number of possible digits that could go in    
c                  each site: choose from smallest list. try/backtrack
c                     
c          tree exploration from moven()
c
c        enter initial position in data statement -- copied to "in"
c
c    Note: code can find maxsol alternative answers (if they exist)
c   search order can be varied -- see S_O (default kchoice=0, 1,...,movemax)
c
c             (c) Chris Michael 4-2014
c------------------------------------------------------

c *******************************************************
c  Number of solutions,   printout volume (default 1; 2 is most; 0 least)
      parameter (maxsol=10, iprin=1)
c *******************************************************

      integer in(9,9), i0(9,9,82),move(82),movemax(82),nmove(81,82)
      integer  movei(81,82), movej(81,82),in1(9,9),in2(9,9),in3(9,9)
     &   ,in4(9,9) ,in5(9,9) ,in6(9,9) ,in7(9,9)
      integer isol(9,9,maxsol)
      logical isok, isdiff, isnew
       external isdiff

c examples to test
c                        hard 782   2 branches
       data in1 /
     &0,0,8,0,6,0,0,2,9,
     &0,2,0,3,0,0,0,7,8,
     &9,0,1,8,0,0,0,0,0,
     &0,6,7,0,0,5,0,0,0,
     &5,0,0,0,0,0,0,0,3,
     &0,0,0,4,0,0,7,6,0,
     &0,0,0,0,0,9,6,0,2,
     &2,9,0,0,0,3,0,1,0,
     &8,1,0,0,5,0,9,0,0 /

c                           v hard  [ has 2 2* branches ]
       data in2 /
     &0,1,0,0,2,0,0,5,0,
     &9,4,0,0,0,0,0,2,6,
     &0,0,6,0,4,0,3,0,0,
     &0,0,1,0,0,0,8,0,0,
     &6,0,0,3,8,4,0,0,9,
     &0,0,4,0,0,0,6,0,0,
     &0,0,5,0,9,0,2,0,0,
     &1,7,0,0,0,0,0,6,4,
     &0,3,0,0,1,0,0,8,0 /

c  has no valid solution (now tweaked in setn0 to find this)
       data in3 /
     &4,0,0,0,0,0,0,0,0,
     &0,0,0,0,0,0,0,0,0,
     &0,0,0,0,0,0,0,0,0,
     &0,0,0,0,4,0,0,0,0,
     &0,9,0,0,0,0,0,0,0,
     &0,0,0,0,0,0,0,4,0,
     &0,0,0,0,0,0,0,0,0,
     &0,0,4,0,0,0,0,0,0,
     &0,0,0,0,0,0,0,0,0 /
c  very hard 1019 [ seems unique path ]
       data in4 /
     &0,0,0,0,0,0,0,0,0,
     &0,8,0,0,0,0,0,3,4,
     &0,3,0,8,5,0,7,0,9,
     &8,0,0,0,0,7,6,0,0,
     &7,0,0,4,0,5,0,0,8,
     &0,0,9,1,0,0,0,0,5,
     &2,0,3,0,6,9,0,4,0,
     &1,6,0,0,0,0,0,2,0,
     &0,0,0,0,0,0,0,0,0 /

c  extremely hard  [ 10 steps ahead  ] 2073 moves in default order
       data in7 /
     &8,0,0,0,0,0,0,0,0,
     &0,0,3,6,0,0,0,0,0,
     &0,7,0,0,9,0,2,0,0,
     &0,5,0,0,0,7,0,0,0,
     &0,0,0,0,4,5,7,0,0,
     &0,0,0,1,0,0,0,3,0,
     &0,0,1,0,0,0,0,6,8,
     &0,0,8,5,0,0,0,1,0,
     &0,9,0,0,0,0,4,0,0 /

c                            null  lots of possible solutions
       data in5 /81*0/
c                            over 10000 possible solutions
       data in6 /1,2,3,4,5,6,7,8,9,9,8*0,8,8*0,7,8*0,6,8*0,5,8*0,
     $4,8*0,3,8*0,2,8*0/

c   choose case to evaluate
c  **********************************
      call copyn(in2,in)
c  **********************************

       nsol=0

       call check(in,isok)

      if( isok) write(*,*)' initial: (valid)'
       call iprint(in)
      if( .not. isok) write(*,*)' initial: (invalid) '
      if( .not. isok) STOP ' initial: (invalid) '



      nct=0
      do i=1,9
       do j=1,9
        if( in(i,j) .eq. 0) nct=nct+1
       enddo
      enddo
      write(*,*)' number to fill ',nct
      write(*,*)

      call copyn(in,i0(1,1,1))
 
c                                      find possible locations for "n"
c                                      first move
      call setn0(i0(1,1,1),nmove(1,1),ict,movei(1,1),movej(1,1))
      if( ict .eq. -1) STOP ' no possible move'
      if( ict .ne. 1) 
     & call setn(i0(1,1,1),nmove(1,1),ict,movei(1,1),movej(1,1))
       if( iprin .gt. 1) 
     &   write(*,901)ict,(nmove(j,1),j=1,ict)
 901   format(' number initial moves= ',i2,' digit= ',9i2,/)
      move(1)=1
      movemax(1)=ict

       il=1
       im=0
c                       endless loop (effectively recursive)
100     continue
         m=move(il)
        if( ((iprin .gt. 0 .and. movemax(il) .gt. 1)
     &   .or.(iprin .gt. 1)) .and. (move(il) .le. movemax(il))  )
     &    write(*,94) il,nmove(m,il),move(il),movemax(il)
 94       format(' lev=',i2,' n=',i1,' move=',i2,' of ',i2)

        if( move(il) .le. movemax(il) ) im=im+1

        call moven(i0,il,move,movemax,nmove,movei,movej)
c                                            care: il can be 82 here 
c                                            care: moven can backtrack
        if( movemax(il) .eq. 99 ) then 
c                                         all places filled
           if( nsol .eq. 0) then 
c                                                first solution
             write(*,*)' no. moves=',im,'  solution no.=',nsol+1
             call iprint(i0(1,1,il))
             call copyn(i0(1,1,il),isol(1,1,1))
             nsol=1
           else
c                                  are further solutions different?
              isnew=.true.
              do ipr=1,nsol
               isnew=isnew .and.  isdiff(i0(1,1,il), isol(1,1,ipr) ) 
              enddo
             if( isnew ) then
              call copyn(i0(1,1,il),isol(1,1,nsol+1))
              write(*,*)' no. moves=',im,'  solution no.=',nsol+1
              call iprint(i0(1,1,il))
              nsol=nsol+1
             endif
           endif
c option to look for further solutions --  back track
           il=il-1
c                         moven backtracks if move(il) > movemax(il)
           move(il)=move(il)+1

           if(nsol .ge. maxsol) write(*,*) maxsol,'  solutions found'
           if(nsol .ge. maxsol) STOP ' maxsol solutions found'
        endif

        goto 100

      end

c--------------------------------------------
c
c check valid (if has zeroes treat as irrelevant)
c
c--------------------------------------------
      subroutine check(in,isok)
      integer in(9,9)
      logical isok


      isok=.true.
c                         rows
      do i=1,9

       do n=1,9
        nn=0
         do j=1,9
          if(in(i,j) .eq. n) nn=nn+1
         enddo
        if ( nn .gt. 1) then
         isok=.false.
         return
        endif
       enddo
      enddo

c                       columns
      do j=1,9

       do n=1,9
        nn=0
         do i=1,9
          if(in(i,j) .eq. n) nn=nn+1
         enddo
       if ( nn .gt. 1) then 
          isok=.false.
          return
       endif
       enddo
      enddo

c                  check 3*3 blocks
      do i3=1,3
       do j3=1,3

       do n=1,9
        nn=0
        do i=1,3
         do j=1,3
          if(in(3*(i3-1)+i,3*(j3-1)+j) .eq. n) nn=nn+1
         enddo
        enddo
       if ( nn .gt. 1) then
        isok=.false.
        return
       endif
         enddo
 
       Enddo
      Enddo


      return
      end


c-----------------------------------------------------
c
c       find optimum locations for a  "move" in i0()
c
c       place digit nmove() in location movei(),movej()
c
c       in movemax ways;   movemax=0 if no options available
c                          movemax=99 if all done
c
c     count and store possible equivalent moves with least options
c
c-----------------------------------------------------
       subroutine setn(i0,nmove,movemax,movei,movej)
       parameter (iprin=1)
       integer i0(9,9)
       integer  movei(81), movej(81),nmove(81),movemax
       integer mti(789,10),mtj(789,10),nt(789,10),nhit(81)
       integer ipgrid(9,9,9)
       logical isok
       character cns(9)*9, cns0(9)*5, cns1(9)*4

c test if all done
        izero=0
        do it=1,9
         do jt=1,9
          if(i0(it,jt) .eq. 0)  izero=izero+1
         enddo
        enddo
        if( izero .eq. 0) then 
          movemax=99
          return   
        endif

c             find locations with least choice (nct is number of choices)
c     try 1,..9 and see how many are possible at each location

        nct=10
        k=1
        do it=1,9
         do jt=1,9
          do n=1,9
             ipgrid(it,jt,n)=0
          enddo
          if(i0(it,jt) .eq. 0)then

          icount =0
          do n=1,9
             ipgrid(it,jt,n)=0
           i0(it,jt)=n
           call check(i0,isok)   
           if(isok) then
             icount =icount +1
             nhit(icount)=n
             ipgrid(it,jt,n)=1
           endif
          enddo

          if(icount  .eq. 0) then
c                            no option for one site -- so backtrack
c           write(*,*)it,jt,nct,k,i0(it,jt)
            movemax=0
            RETURN
          endif

c                                       reset count k for a better case
          if(icount  .lt. nct) k=1
c                                       save possible moves (only some used)

          nct=min(nct,icount )
          do kv=0,nct-1
           if(k+kv .gt. 729 ) STOP 'k+kv too big'
           mti(k+kv,nct)=it
           mtj(k+kv,nct)=jt
           nt(k+kv,nct)= nhit(kv+1)
c           write(*,*)nct,it,jt,icount ,k+kv,nhit(kv+1)
          enddo
          if(icount  .eq. nct) k=k+nct
           i0(it,jt)=0
c -------------------end of empty space loop -----------
          else
           cycle
          endif      

         enddo
        enddo

         if(iprin .gt. 0) then

          call iprint(i0(1,1))

         do it=1,9
          do jt=1,9
c                                collect numbers
            in=0
           cns0(jt)="....."
           cns1(jt)="    "
           do n=1,9
            if(ipgrid(jt,it,n) .eq. 1) Then 
              cns(jt)(n:n)=char(n+48)
               in=in+1
               if( in .lt. 6) cns0(jt)(in:in)=cns(jt)(n:n)
               if( in .gt. 5) cns1(jt)(in:in)=cns(jt)(n:n)
             else
              cns(jt)(n:n)="."
             endif
           enddo
c either output 5 plus 4; or just give first 5 matches     
c            cns0(jt)=cns(jt)(1:5)
c            cns1(jt)=cns(jt)(6:9)

          enddo
    
c          write(*,*)
           write(*,'(3(3(a5,2x),2x))') (cns0(jt),jt=1,9)
           write(*,'(3(3(a4,3x),2x))') (cns1(jt),jt=1,9)
c          write(*,'(3(3(3i1,x),x) )') 
c     & (((ipgrid(jt,it,nc,nr),nc=1,3),jt=1,9),nr=1,3)
         enddo
         write(*,*)
         write(*,*)

         endif


c select among equivalent best moves: nct possible digits at nk sites

c could consider forced moves (nct=1) in any order

        if ( nct .gt. 0) then
           movemax=nct
            nk= (k-1)/nct 
c   S_O     could choose any of equivalent choices kchoice=0,..nk-1
            kchoice=0
            kchoice=mod(kchoice,nk)
            koff= nct* kchoice
           do kv=1,nct
c    could select nct possible digits in any order (code uses order 1,,,nct)
            kvs=mod(kv-1,nct)+1          
            nmove(kv)=nt(kvs+koff,nct)
            movei(kv)=mti(kvs+koff,nct)
            movej(kv)=mtj(kvs+koff,nct)
           enddo
        endif

       return
       end
c-----------------------------------------------------
c
c       find optimum locations for a  "move" in i0()
c
c       place digit nmove() in location movei(),movej() if unique
c        [also test for no solution cases ]
c
c       in movemax ways;   movemax=-1 if no options available
c                          movemax=1  if unique move exists
c                          movemax=0  if no unique move
c                          movemax=99 if all done
c
c  search sub-blocks; rows and columns first
c
c-----------------------------------------------------
       subroutine setn0(i0,nmove,movemax,movei,movej)
       integer i0(9,9)
       integer  movei(81), movej(81),nmove(81),movemax
       integer mti(789,10),mtj(789,10),nt(789,10),nhit(81)
       logical isok

c test if all done
        izero=0
        do it=1,9
         do jt=1,9
          if(i0(it,jt) .eq. 0)  izero=izero+1
         enddo
        enddo
        if( izero .eq. 0) then 
          movemax=99
          return   
        endif

c   take a digit and see how many places it could fit
c    in each sub-block; row; col


        do n=1,9

c                                  search block
        do ib=1,3
         do jb=1,3

          icount =0
          isn=0
          isz=0
         do ir=1,3
          do jr=1,3
           it=3*(ib-1)+ir
           jt=3*(jb-1)+jr
           if( i0(it,jt) .eq. n) isn=1 
           if( i0(it,jt) .eq. 0) then
            isz=1
            i0(it,jt)=n
            call check(i0,isok)   
            if(isok) then
             icount =icount +1
             nmove(icount)=n
             movei(icount)=it
             movej(icount)=jt
            endif
            i0(it,jt)=0
           endif
          enddo
          enddo
c                     if no possible solution return to back-track
          if( icount .eq. 0 .and. isn .eq. 0 .and. isz .eq. 1) then

           nmove(1)=n
           movemax=-1
           return
          endif 
c                                 if unique move -- do it
          if( icount .eq. 1) then
           movemax=1
           RETURN
          endif

          enddo
          enddo

c                             search column

         do it=1,9
          icount =0
          isn=0
          isz=0

          do jt=1,9
           if( i0(it,jt) .eq. n)  isn=1 
           if(i0(it,jt) .eq. 0) then
            isz=1
            i0(it,jt)=n
            call check(i0,isok)   
            if(isok) then
             icount =icount +1
             nmove(icount)=n
             movei(icount)=it
             movej(icount)=jt
            endif
            i0(it,jt)=0
           endif
          enddo
c                     if no possible solution
          if( icount .eq. 0 .and. isn .eq. 0 .and. isz .eq. 1) then

           movemax=-1
           nmove(1)=n
           return
          endif 
c                                 if unique move -- do it
          if( icount .eq. 1) then
           movemax=1
c           i0(movei(1),movej(1))=n
           RETURN
          endif

          enddo

c                             search row

         do jt=1,9
          icount =0
          isz=0
          isn=0
 
          do it=1,9
           if( i0(it,jt) .eq. n)  isn=1 
           if(i0(it,jt) .eq. 0) then
            isz=0
            i0(it,jt)=n
            call check(i0,isok)   
            if(isok) then
             icount =icount +1
             nmove(icount)=n
             movei(icount)=it
             movej(icount)=jt
            endif
            i0(it,jt)=0
           endif
          enddo
c                     if no possible solution
          if( icount .eq. 0 .and. isn .eq. 0 .and. isz .eq. 1) then

           movemax=-1
           nmove(1)=n
           return
          endif 
c                                 if unique move -- do it
          if( icount .eq. 1) then
           movemax=1
c           i0(movei(1),movej(1))=n
           RETURN
          endif

          enddo


          enddo
c signal no unique move
          movemax=0
       return
       end
c---------------------------------------------------
copy
c---------------------------------------------------
       subroutine copyn(in,i0)
       integer in(9,9),i0(9,9)
       logical isok
       do i=1,9
       do j=1,9
        i0(i,j)=in(i,j)
       enddo
       enddo

       return
       end

c---------------------------------------------------
c different
c---------------------------------------------------
       logical function isdiff(in,i0)
       integer in(9,9),i0(9,9)
        isdiff=.false.
       do i=1,9
       do j=1,9
        if(i0(i,j) .ne. in(i,j)) isdiff=.true.
        if( isdiff) return
       enddo
       enddo

       return
       end

c----------------------------------------------------------------
c
c   "recursive" selection of a digit
c   "move" means place a digit in 1 location and increemnt level (il)
c
c       move // select digit to set at next level 
c
c----------------------------------------------------------------

      subroutine moven(i0,il,move,movemax,nmove,movei,movej)
       integer i0(9,9,82),move(82),movemax(82),nmove(81,82)
       integer  movei(81,82), movej(81,82)
c                                  current proposed move:   move(il)


       m=move(il)

c                                  test if move will be valid
       if( m .gt. movemax(il)) then
c                                                     back-track
        il=il-1
         if( il .lt. 1) write(*,*) ' all possible moves explored '
         if( il .lt. 1) STOP

c                                                     next move
        move(il)=move(il)+1

        RETURN
       endif

c                                                 make  move
         call copyn(i0(1,1,il),i0(1,1,il+1))
         i0(movei(m,il),movej(m,il),il+1)=nmove(m,il)

c                                                 explore one level further

c                           find and  count possible locations for next move
c                   look for unique moves:
       call setn0(i0(1,1,il+1),nmove(1,il+1),ict,
     &            movei(1,il+1),movej(1,il+1))
          if( ict .eq. 99) then
c                                                  solution found  
           movemax(il+1)=99
           il=il+1           
           return
          endif
c                             no valid move -- try next previous move
          if( ict .eq. -1) then
           write(*,'(a,i2,a,i1)') ' lev=',il+1,
     &       ' no valid solution for n=', nmove(1,il+1)
           move(il)=move(il)+1
           return
          endif
c                         no unique move - so collect possible moves
       if( ict .ne. 1) 
     &  call setn(i0(1,1,il+1),nmove(1,il+1),ict,
     &            movei(1,il+1),movej(1,il+1))


c       write(*,*)' number ways ',ict

c                           if no ways then dead-end and try next move
       if( ict .lt. 1) then
        move(il)=move(il)+1
c                                  check <=  movemax at top of code
        RETURN        
       endif

c set up to move at deeper level [ nmove movei movej set already ]

       movemax(il+1)=ict
       move(il+1)=1
       il=il+1
       return
       end
c------------------------------------------------------
c
c  print out
c
c------------------------------------------------------
      subroutine iprint(i0)
      integer i0(9,9)
        do i=1,9
          write(*,90) (i0(j,i),j=1,9)
 90       format(9i2)
        enddo
        return
        end