Solve Sudokus with R and PicoSAT

At least on number per cell

For all \(i,j\) we build a clause \(\bigvee \limits_{v=1}^{9} x_{i,j,v}\). So in order for this clause to evaluate to TRUE, one of its variables needs to be TRUE and thus at least one number is present per cell.

# we nest the data per cell i,j
per_cell <- nest(grid_helper, v, index)$data

# then the clauses are straight forward.
at_least_one <- map(per_cell, ~.x$index)
head(at_least_one, 3)

## [[1]]
## [1] 1 2 3 4 5 6 7 8 9
## 
## [[2]]
## [1] 10 11 12 13 14 15 16 17 18
## 
## [[3]]
## [1] 19 20 21 22 23 24 25 26 27

In total we have 81 such clauses.

At most one number per cell

In order to have exactly one number per cell, we create a clause for each cell and all combinations of two numbers in that cell. This clause forbids both numbers to be present at the same time.

So for each cell \(i,j\) and each pair of values \(1\leq v_1 < v_2 \leq 9\): \(\neg(x_{i,j,v_1} \wedge x_{i,j,v_2})\).

The code looks a bit strange, but just goes through our data.frame and selects the correct pairs of variables.

at_most_one  <- flatten(map(per_cell, function(ij_values) {
  flatten(map(1:8, function(v1) {
    map((v1 + 1):9, function(v2) {
      -1 * c(
        ij_values[ij_values$v == v1, ]$index,
        ij_values[ij_values$v == v2, ]$index)
    })
  }))
}))
head(at_most_one, 3)

## [[1]]
## [1] -1 -2
## 
## [[2]]
## [1] -1 -3
## 
## [[3]]
## [1] -1 -4

This gives us 2916 clauses.

No two equal numbers per row

Again for each row \(i\) and each value \(v\) we find all tuples for columns and negate them: \(\forall1\leq j_1 < j_2 \leq 9\): \(\neg(x_{i,j_1,v} \wedge x_{i,j_2,v})\)

no_two_per_line <- flatten(map(1:9, function(i) {
  line_df <- grid_helper[grid_helper$i == i, ]
  flatten(map(1:9, function(v) {
    line_df_v <- line_df[line_df$v == v, ]
    flatten(map(1:8, function(j1) {
      map((j1 + 1):9, function(j2) {
          -1 * c(
            line_df_v[line_df_v$j == j1, ]$index,
            line_df_v[line_df_v$j == j2, ]$index
          )
      })
    }))  
  }))
}))
head(no_two_per_line, 3)

## [[1]]
## [1]  -1 -10
## 
## [[2]]
## [1]  -1 -19
## 
## [[3]]
## [1]  -1 -28

In total we have 2916 clauses.

No two equal numbers per column

The same can be done for each column:

no_two_per_col <- flatten(map(1:9, function(j) {
  col_df <- grid_helper[grid_helper$j == j, ]
  flatten(map(1:9, function(v) {
    col_df_v <- col_df[col_df$v == v, ]
    flatten(map(1:8, function(i1) {
      map((i1 + 1):9, function(i2) {
          -1 * c(
            col_df_v[col_df_v$i == i1, ]$index,
            col_df_v[col_df_v$i == i2, ]$index
          )
      })
    }))  
  }))
}))
head(no_two_per_col, 3)

## [[1]]
## [1]  -1 -82
## 
## [[2]]
## [1]   -1 -163
## 
## [[3]]
## [1]   -1 -244

In total we have 2916 clauses.

No two equal numbers per block

The last constraint forbids the same number being in a block twice. Again we loop through each block, value and pairs of variables. Here the block column in the grid_helper comes in handy.

no_two_per_block <- flatten(map(1:9, function(b) {
  line_df <- grid_helper[grid_helper$block == b, ]
  flatten(map(1:9, function(v) {
    indexes <- line_df[line_df$v == v, ]$index
    flatten(map(1:8, function(i1) {
      map((i1 + 1):9, function(i2) {
          -1 * c(
            indexes[i1],
            indexes[i2]
          )
      })
    }))  
  }))
}))
head(no_two_per_block, 3)

## [[1]]
## [1]  -1 -10
## 
## [[2]]
## [1]  -1 -19
## 
## [[3]]
## [1]  -1 -82

Again 2916 clauses were added.

Solve it

Having prepared all the clauses, we can combine them to one giant CNF.

formula <- c(
  at_least_one,
  at_most_one,
  no_two_per_line,
  no_two_per_col,
  no_two_per_block
)
length(formula)

## [1] 11745

Even though we create a fairly large number of clauses, PicoSAT is able to find a solution in less than a second.

This CNF can then be passed to PicoSAT:

library(rpicosat)
result <- picosat_sat(formula)
result

## Variables: 729
## Clauses: 11745
## Solver status: PICOSAT_SATISFIABLE

The result is also a data.frame and we use tidyr to spread the \(j\) as columns. This then roughly looks like a Sudoku ;)

format_result <- function(x) {
  df <- x %>%
    right_join(grid_helper, by = c("variable" = "index")) %>% 
    filter(value == TRUE) %>% 
    select(i, j, v) %>% 
    spread(j, v)
  select(df, -i) %>% 
    knitr::kable(row.names = 1:9)
}

format_result(result)

Of course a Sudoku is only fun if some numbers are set beforehand. We can fix this by setting the respective variables to TRUE. This can be done by passing an integer vector to picosat_sat. Within picosat these are called assumptions.

So let us assume the first row is fixed ascending from 1 to 9 and the last col is set to c(9, 1:8).

fixed_vars1 <- filter(grid_helper, i == 1, j == v)$index
fixed_vars2 <- filter(grid_helper, j == 9, i == v + 1, i > 1)$index
result <- picosat_sat(formula, assumptions = c(fixed_vars1, fixed_vars2))

With these assumed starting values, the resulting Sudoku looks like this:

format_result(result)

Finding all possible solutions

Now the last question one could ask is: is this solution unique? To check this we can enumerate all possible solutions to a Sudoku with these assumptions. In order to do that, whenever we find a solution, we build another formula disallowing the found solution until the formula becomes unsatisfiable. This can easily be done by adding the negation of the solution as another clause (read=“ and not this solution”).

As an example we will enumerate all solutions to the Sudoku example on Wikipedia by by Tim Stellmach (license CC-0):

However as this Sudoku is unique, we leave out the numbers in the center block:

sample_sudoku <- with(grid_helper, index[
  i == 1 & j == 1 & v == 5 | 
  i == 1 & j == 2 & v == 3 | 
  i == 2 & j == 1 & v == 6 | 
  i == 3 & j == 2 & v == 9 | 
  i == 3 & j == 3 & v == 8 | 
                                                      
  i == 1 & j == 5 & v == 7 | 
  i == 2 & j == 4 & v == 1 | 
  i == 2 & j == 5 & v == 9 | 
  i == 2 & j == 6 & v == 5 |
  
  i == 3 & j == 8 & v == 6 | 
                                                      
  i == 4 & j == 1 & v == 8 | 
  i == 5 & j == 1 & v == 4 | 
  i == 6 & j == 1 & v == 7 |
                                                      
  #i == 4 & j == 5 & v == 6 | 
  #i == 5 & j == 4 & v == 8 | 
  #i == 5 & j == 6 & v == 3 | 
  #i == 6 & j == 5 & v == 2 |
  
  i == 4 & j == 9 & v == 3 | 
  i == 5 & j == 9 & v == 1 | 
  i == 6 & j == 9 & v == 6 | 
                                                      
  i == 7 & j == 2 & v == 6 |
  
  i == 8 & j == 4 & v == 4 | 
  i == 8 & j == 5 & v == 1 | 
  i == 8 & j == 6 & v == 9 | 
  i == 9 & j == 5 & v == 8 |
                                                      
  i == 7 & j == 7 & v == 2 | 
  i == 7 & j == 8 & v == 8 | 
  i == 8 & j == 9 & v == 5 | 
  i == 9 & j == 8 & v == 7 | 
  i == 9 & j == 9 & v == 9])

current_formula <- formula
results <- list()
no_results <- 0
system.time(while(TRUE) {
  result <- picosat_sat(current_formula, assumptions = sample_sudoku)
  if (picosat_solution_status(result) == "PICOSAT_SATISFIABLE") {
    no_results <- no_results + 1
    results[[no_results]] <- result
    new_clause <- -1 * (result$variable * ifelse(result$value, 1L, -1L))
    current_formula <- c(list(new_clause), current_formula)
  } else {
    break
  }
})

##    user  system elapsed 
##   1.603   0.033   1.672

The resulting Sudoku has 63 solutions and is therefore not unique.