Efficient Sudoku solver in 13 lines of Caml

open Array

let bits = init 9 (fun b -> b, 1 lsl b)
let group = init 81 (fun i -> [| i/9; 9 + i mod 9; 18 + (i/27)*3 + (i mod 9) / 3 |])
let card = init 512 (fun i -> fold_left (fun a (b,s) -> if i land s > 0 then a+1 else a) 0 bits)

let rec solve grid = 
   let free = make 27 511 in
   iteri (fun i v -> iter (fun g -> free.(g) <- free.(g) - (1 lsl (v-1))) group.(i)) grid;
   let poss = mapi (fun i v -> (i,v), fold_left (fun a g -> a land free.(g)) 511 group.(i)) grid in
   let m,c = fold_left (fun (m,c) ((i,v),p) -> if v = 0 && card.(p) < m then card.(p),i else m,c) (10,-1) poss in
   let subsolve value = let grid2 = copy grid in grid2.(c) <- value; solve grid2 in
   if m = 10 then (iteri (fun i v -> print_int v; if i mod 9 = 8 then print_char '\n') grid; print_char '\n')
   else if m > 0 then iter (fun (b,s) -> if (snd poss.(c)) land s > 0 then (subsolve (b+1))) bits

let _ = 
   let inp() = int_of_char (input_char stdin) - 48 in
   solve (init 81 (fun i -> (if i > 0 && i mod 9 = 0 then ignore (inp())); inp()))

a rather efficient Sudoku solver
written in 13 lines of OCaml (some 1024 bytes to be precise)
using nothing else than modules Pervasive and Array (no advanced data structures)
it displays all solutions for the input grid (easy modification if you only want one, read further)

Warning: this is just an experiment to study how compact can the implementation of an algorithm can be.
It is not a good idea in general to write source code that looks like this one!

How does it work?

The algorithm looks for the cell which has the more constraint on it, then there is 3 possible cases:

if a cell has no possibilities left, backtrack
if all cells are labeled, print grid, and return
otherwise try all legal values for this cell, backtracking

About the program

Some comments about the code:

it mixes a lot functionnal and imperative programming.
it uses only arrays for data structures.
it indexes linearly all cells and all groups of constraints.
it uses integer division and modulo to work out relationship between cells and groups of constraints.
it uses logical operations on bits to perform the insersection of the constraints.

About these comments:

Point 1) => this is usual when programming algorithm in Caml
Point 2) => it means that only basic algorithmics is involved
Point 3), 4) and 5) => these are standard techniques when implementing algorithms

Some statistics about the code: upon the 8 arrays used, there are plenty of operations performed:

4 Array.init, 1 Array.make
1 Array.mapi, 1 Array.copy
4 Array.fold_left, 2 Array.iteri

1) The array bits contains the sequence of values (i,2^i) for i from 0 to 9.
   let bits = init 9 (fun b -> b, 1 lsl b)

2) The array group contains for each cell the indices of the three groups of constraints containing that cell.
   let group = init 81 (fun i -> [| i/9; 9 + i mod 9; 18 + (i/27)*3 + (i mod 9) / 3 |])

3) The card array stores for each integer the number of bits set to 1 in its binary representation.   
   let card = init 512 (fun i -> fold_left (fun a (b,s) -> if i land s > 0 then a+1 else a) 0 bits)

4) The recursive solve function is the backtrack function   
   let rec solve grid = 
   
5) The array free contains for each group of constraints a field of bit (an int) 
the k-ith bit is set if and only if the value (k+1) is not yet used in the group.
   let free = make 27 511 in
   
6) For each cell, for each group containing that cell, we unset the bit corresponding to the value of the cell.
   iteri (fun i v -> iter (fun g -> free.(g) <- free.(g) - (1 lsl (v-1))) group.(i)) grid;
   
7) For each cell, the possibilities left are the intersection of the values unused 
in the three groups of contraints containing the cell
   let poss = mapi (fun i v -> (i,v), fold_left (fun a g -> a land free.(g)) 511 group.(i)) grid in

8) We extract the index of the cell c which has not yet a value and which minimises the number of possibilities.
   let m,c = fold_left (fun (m,c) ((i,v),p) -> if v = 0 && card.(p) < m then card.(p),i else m,c) (10,-1) poss in

9) In case all cells have a value, we have a solution, we print it.
   if m = 10 then (iteri (fun i v -> print_int v; if i mod 9 = 8 then print_char '\n') grid; print_char '\n')

10) In case m = 0, at least one cell has no possibility left, there is no solution, we simply return.
Otherwise, m > 0, we try recursively all legal values for the cell c selected previously.
   else if m > 0 then iter (fun (b,s) -> if (snd poss.(c)) land s > 0 then (subsolve (b+1))) bits

11) The subsolve value copy the grid, place a value, and call recursively the function solve 
   let subsolve value = let grid2 = copy grid in grid2.(c) <- value; solve grid2 in
   
12) To read one character and obtain the corresponding digit:
   let inp() = int_of_char (input_char stdin) - 48 in
   
13) To call function solve on the input grid:   
   solve (init 81 (fun i -> (if i > 0 && i mod 9 = 0 then ignore (inp())); inp()))

A) If you want to see only one solution, you can add an exit 0. To do that, replace line 9:
   if m = 10 then (iteri (fun i v -> print_int v; if i mod 9 = 8 then print_char '\n') grid; print_char '\n')
   with:
   if m = 10 then (iteri (fun i v -> print_int v; if i mod 9 = 8 then print_char '\n') grid; print_char '\n'; exit 0)

B) If you only want to count the number of solutions, do the following:

   1) at the top of your source, add the line:
      let nb_solutions = ref 0

   2) in the middle of the source, replace line 9:
      if m = 10 then (iteri (fun i v -> print_int v; if i mod 9 = 8 then print_char '\n') grid; print_char '\n')
      with:
      if m = 10 then incr nb_solutions

   3) at the bottom of your source, add the line:
      print_int !nb_solutions

1) Compile the program:

ocamlopt -o sudoku code.ml

2) Run the program:

./sudoku

3) Type in your input problem. Here is an example:

4) The solution is displayed. For the example above:

This example is copied from ffconsultancy website, and it has a unique solution.

Note that lines of input must be ended by \n. If your lines end with \r\n, you need to change the last line to:

   solve (init 81 (fun i -> (if i > 0 && i mod 9 = 0 then ignore (inp(),inp())); inp()))

If you have any question, please mail me on: arthur@france-ioi.org

This code is free for anyone to play with it, simply mention my name when using it.