!************************************************************************
!*                                                                      *
!* PROGRAM SUDOKU - a recursive, parallel, Sudoku puzzle solver.        *
!*                                                                      *
!* This Sudoku puzzle solver is able to solve 2x2, 3x3, 4x4, 5x5        *
!* and 6x6 Sudoku puzzles, using the following "digits":                *
!*                                                                      *
!*   2x2 - 1..4                                                         *
!*   3x3 - 1..9                                                         *
!*   4x4 - 0..9, A..F (the "hexadecimal" digits)                        *
!*   5x5 - 0..9, A..O                                                   *
!*   6x6 - 0..9, A..Z                                                   *
!*                                                                      *
!* It also shows ALL solutions to the given clues, so it tells whether  *
!* a puzzle is a genuine Sudoku puzzle (having a unique solution).      *
!*                                                                      *
!* As a consequence, an empty board will give you ALL Sudoku layouts    *
!* possible (if disk space permits).                                    *
!* For 2x2 sudokus this is achievable, as there are 288 layouts;        *
!* an empty board is 4 rows of 8 spaces.                                *
!*                                                                      *
!* Compile:                                                             *
!*   gfortran -g -O3 -fstack-arrays -fopenmp -o sudoku sudoku.f90       *
!*                                                                      *
!* Run:                                                                 *
!*   export OMP_NESTED=true                                             *
!*   export OMP_THREAD_LIMIT=20            # On my machine              *
!*   ./sudoku <arguments>                                               *
!*                                                                      *
!************************************************************************
!
! Integer kind used for "digits" and size of the board.
! Offset to deal with the fact that the 2x2 and 3x3 boards are "special".
!
module settings
   integer, parameter :: ik = selected_int_kind(2) ! -100 .. +100
   integer(kind=ik)   :: ioff = 0 ! Internally, the numbers are 'ioff'.
   integer            :: verbose = 0, solution_number = 0
end module settings
!
! Two functions to convert "digits" to numbers and vice versa
! in the internal representation of the board.
! These functions are "elemental", which means that they will
! work on array arguments, returning a same-shaped array result.
!
module converters
   interface
      pure elemental integer(kind=ik) function alf2num(a)
         use settings, only : ik
         character(len=1), intent(in) :: a
      end function alf2num
   end interface
   interface
      pure elemental character(len=1) function num2alf(i)
         use settings, only : ik
         integer(kind=ik), intent(in) :: i
      end function num2alf
   end interface
end module converters
!
pure elemental integer(kind=ik) function alf2num(a)
   use settings, only : ik, ioff
   implicit none
   character(len=1), intent(in) :: a
   select case (a)
   case (' ')
      alf2num = 0
   case ('0':'9')
      alf2num = ichar(a) - ichar('0') + ioff
   case ('A':'Z')
      alf2num = ichar(a) - ichar('A') + 10 + ioff
   end select
end function alf2num
!
pure elemental character(len=1) function num2alf(i)
   use settings, only : ik, ioff
   implicit none
   integer(kind=ik), intent(in) :: i
   select case (i)
   case ( 0)
      num2alf = ' '
   case ( 1:10)
      num2alf = char(ichar('0') + i - ioff)
   case (11:36)
      num2alf = char(ichar('A') + i - 10 - ioff)
   end select
end function num2alf
!
! Subroutine to read a board (lines of space separated "digits") from stdin.
!
subroutine read_board(board, n)
   use settings, only : ik
   use converters, only : alf2num
   implicit none
   integer(kind=ik), intent(out) :: board(n**2,n**2)
   integer(kind=ik), intent(in)  :: n
   character(len=1)              :: cb(n**2,n**2)
   character(len=10)             :: fmt
   integer(kind=ik)              :: i
   write(fmt, '(a,i2.2,a)') '(',n**2,'(1x,a))'
   do i = 1, n**2
      read fmt, cb(:,i)
   enddo
   board = alf2num(cb)
end subroutine read_board
!
! Subroutine to write a board (not necessarily completed) to stdout.
! Lines of space-separated "digits".
!
subroutine write_board(board, n)
   use settings, only : ik
   use converters, only : num2alf
   implicit none
   integer(kind=ik), intent(in) :: board(n**2,n**2)
   integer(kind=ik), intent(in) :: n
   character(len=1)             :: cb(n**2,n**2)
   character(len=10)            :: fmt
   integer(kind=ik)             :: i
   write(fmt, '(a,i2.2,a)') '(',n**2,'(1x,a))'
   cb = num2alf(board)
   do i = 1, n**2
      print fmt, cb(:,i)
   enddo
end subroutine write_board
!
! Auxiliary routine to update the array of values on the board
! that are made impossible after setting ('column','row') to 'value'.
!
subroutine update_possibles(possibles, n, value, column, row)
   use settings, only : ik
   implicit none
   integer(kind=ik), intent(inout) :: possibles(n**2, n**2, n**2)
   integer(kind=ik), intent(in)    :: n, value, column, row
   integer(kind=ik)                :: iupper, ilower, jupper, jlower
   possibles(:    , column, row) = 0  ! Nix this position ...
   possibles(value, :,      row) = 0  ! ... and the row ...
   possibles(value, column, :  ) = 0  ! ... and the column ...
   jlower = ((column-1) / n) * n + 1  ! ... and the box ...
   jupper = ((column-1) / n) * n + n  ! ...
   ilower = ((row   -1) / n) * n + 1
   iupper = ((row   -1) / n) * n + n
   possibles(value, jlower:jupper, ilower:iupper) = 0
end subroutine update_possibles
!
! The recursive solver subroutine.
! Its input is an incompletely solved board, and its output
! is either a solution, or a still incomplete board.
! Open squares are signified by a 0 content, so the test for
! having a solution is the count of 0's on the board being zero.
!
recursive subroutine solve(board, n)
   use settings, only : ik, solution_number, verbose
   use converters, only : num2alf
   implicit none
   integer(kind=ik), intent(inout) :: board(n**2,n**2)
   integer(kind=ik), intent(in)    :: n
   integer(kind=ik)                :: save_board(n**2,n**2)
   integer(kind=ik)                :: possibles(n**2,n**2,n**2)
   integer(kind=ik)                :: i, j, k, m, ii, jj, alternatives(n**2)
   logical                         :: found
   possibles = 1
   do i = 1, n**2
      do j = 1, n**2
         if (board(j, i) /= 0) call update_possibles(possibles, n, board(j, i), j, i)
      enddo
   enddo
mainloop: &
   do
      do i = 1, n**2
         do j = 1, n**2
            if (board(j, i) == 0 .and. sum(possibles(:, j, i)) == 1) then
               board(j, i) = findloc(possibles(:, j, i), 1, dim=1)
               if (verbose > 3) print*,'Set column',j,' row',i,' to ',num2alf(board(j, i)),' [no choice left]'
               call update_possibles(possibles, n, board(j, i), j, i)
               cycle mainloop
            endif
         enddo
      enddo
      do i = 1, n**2
         do k = 1, n**2
            if (sum(possibles(k, :, i)) == 1) then
               do j = 1, n**2
                  if (board(j, i) == 0 .and. possibles(k, j, i) == 1) then
                     board(j, i) = k
                     if (verbose > 3) print*,'Set column',j,' row',i,' to ',num2alf(k),' [completes row]'
                     call update_possibles(possibles, n, k, j, i)
                     cycle mainloop
                  endif
               enddo
            endif
         enddo
      enddo
      do j = 1, n**2
         do k = 1, n**2
            if (sum(possibles(k, j, :)) == 1) then
               do i = 1, n**2
                  if (board(j, i) == 0 .and. possibles(k, j, i) == 1) then
                     board(j, i) = k
                     if (verbose > 3) print*,'Set column',j,' row',i,' to ',num2alf(k),' [completes column]'
                     call update_possibles(possibles, n, k, j, i)
                     cycle mainloop
                  endif
               enddo
            endif
         enddo
      enddo
      do i = 1, n**2, n
         do j = 1, n**2, n
            do k = 1, n**2
               if (sum(possibles(k, j:j+n-1, i:i+n-1)) == 1) then
                  do ii = i, i+n-1
                     do jj = j, j+n-1
                        if (board(jj, ii) == 0 .and. possibles(k, jj, ii) == 1) then
                           board(jj, ii) = k
                           if (verbose > 3) print*,'Set column',jj,' row',ii,' to ',num2alf(k),' [completes box]'
                           call update_possibles(possibles, n, k, jj, ii)
                           cycle mainloop
                        endif
                     enddo
                  enddo
               endif
            enddo
         enddo
      enddo
      exit mainloop
   enddo mainloop
   if (count(board == 0) == 0) then          ! If zero, solution found - print it.
!$OMP CRITICAL
      solution_number = solution_number + 1
      if (verbose > 0) then
         print*,'SOLUTION number', solution_number
         call write_board(board, n)
      endif
!$OMP END CRITICAL
      return
   endif
   found = .false.                           ! No:  Need to search for alternatives
outer: &                                     ! (remember to place a known elephant in Cairo)
   do k = 2, n**2                            ! Depth first search (prunes more hopeless paths)
!  do k = n**2, 2, -1                        ! Breadth first search (enables more parallellism)
      do i = 1, n**2
         do j = 1, n**2
            if (board(j, i) == 0) then
               found = sum(possibles(:, j, i)) == k
               if (found) exit outer         ! Found some
            endif
         enddo
      enddo
   enddo outer
   if (.not. found) return                   ! Did not find any - backtrack
   k = 0                                     ! Otherwise, enumerate alternatives
   do m = 1, n**2
      if (possibles(m, j, i) == 1) then
         k = k + 1
         alternatives(k) = m
      endif
   enddo                                     ! and search for results there - in parallel.
   if (verbose > 1 .and. k > 3 .or. verbose > 2) &
      print*,'Possible values for column',j,' row',i,' are',(' ',num2alf(alternatives(m)),m=1,k)
!$OMP PARALLEL DO PRIVATE (save_board), NUM_THREADS (k)
   do m = 1, k
      save_board = board
      save_board(j, i) = alternatives(m)
      if (verbose > 3) print*,'Alternative',m,': set column',j,' row',i,' to ',num2alf(alternatives(m))
      call solve(save_board, n)
   enddo
!$OMP END PARALLEL DO
end subroutine solve
!
! The main program (supports various options, see below).
!
program sudoku
   use settings, only : ik, ioff, solution_number, verbose
   implicit none
   integer(kind=ik), allocatable :: board(:,:)
   integer(kind=ik) :: n = 3
   integer          :: i
   character(len=2) :: arg
   do i = 1, command_argument_count(), 2
      call get_command_argument(i, arg)
      select case (arg)
      case ('-n')             ! Size of the board
         call get_command_argument(i+1, arg)
         read(arg, fmt=*) n
         if (n < 2 .or. n > 6) then
            print*,'N=',n,' is out of range 2..6'
            stop 1
         endif
      case ('-v')             ! Verbose; show all decisions
         call get_command_argument(i+1, arg)
         read(arg, fmt=*) verbose
      case default
         print*,'Unknown option: ',arg
         stop 2
      end select
   enddo
   if (n > 3) ioff = 1
   allocate(board(n**2, n**2))
   call read_board(board, n)
   call write_board(board, n)
   call solve(board, n)
   deallocate(board)
   if (solution_number == 0) then
      print*,'NO SOLUTION'
   else
      if (verbose == 0) then
         print*,'NUMBER OF SOLUTIONS:', solution_number
      else
         print*,'NO FURTHER SOLUTIONS'
      endif
   endif
end program sudoku