Recursion is a powerful technique wherein a function makes calls to itself to solve part of the problem together with a base case that yields an answer for the simplest case. Many problems can be solved in a reasonable amount of time using recursion occasionally when they are coupled with better heuristics. The famous "NQueens problem", "Knight's Tour" and "Boggle words" are just to name a few.

Here is a quick and easy recursive way to solve a Sudoku puzzle

suppose we have a grid as shown below

grid = {{9, 0, 0, 2, 3, 7, 6, 8, 0},
   {0, 2, 0, 8, 4, 0, 0, 7, 3},
   {8, 0, 7, 1, 0, 5, 0, 2, 9},
   {0, 0, 4, 5, 9, 8, 3, 0, 0},
   {2, 0, 0, 0, 0, 1, 0, 0, 6},
   {5, 1, 0, 0, 0, 0, 0, 4, 7},
   {4, 0, 1, 3, 0, 6, 2, 9, 5},
   {0, 5, 0, 9, 1, 0, 7, 3, 8},
   {3, 0, 8, 0, 5, 0, 0, 0, 0}};

now the idea is to replace all zeros by numbers ranging from 1 to 9 in a manner that the row and the column as well as the 3x3 submatrix the element belongs to do not have any repeating integers.

this check can be conveniently implemented by using getSafeNumbers

getSafeNumbers[grid_, {row_, col_}] := 
 Block[{rownum, colnum, submatrix},
  rownum = Cases[grid[[row]], Except[0, _Integer]];
  colnum = Cases[grid[[All, col]], Except[0, _Integer]];
  submatrix = 
   Cases[#, Except[0, _Integer]] &@
    Flatten[Which[row <= 3 && col <= 3, Take[grid, 3, 3],
      row <= 3 && col <= 6, Take[grid, 3, 4 ;; 6], 
      row <= 3 && col <= 9, Take[grid, 3, 7 ;; 9],
      row <= 6 && col <= 3, Take[grid, 4 ;; 6, 3], 
      row <= 6 && col <= 6, Take[grid, 4 ;; 6, 4 ;; 6], 
      row <= 6 && col <= 9, Take[grid, 4 ;; 6, 7 ;; 9],
      row <= 9 && col <= 3, Take[grid, 7 ;; 9, 3], 
      row <= 9 && col <= 6, Take[grid, 7 ;; 9, 4 ;; 6], 
      row <= 9 && col <= 9, Take[grid, 7 ;; 9, 7 ;; 9]]];
  Complement[Range@9, Union[rownum, colnum, submatrix]]
  ] 

All we now require is the recursive calls:

sudoku[grid_] := Module[{zeroposition, safenumbers},
   zeroposition = FirstPosition[grid, 0];
   safenumbers = getSafeNumbers[grid, zeroposition];
   If[Head[zeroposition] =!= Missing ,
    sudoku[ReplacePart[grid, zeroposition -> #]] & /@ safenumbers, 
    Sow[grid]]
   ];

the initialization is done as follows:

(sol = Reap[sudoku[grid]][[2]] // Flatten[#, 2] &;) // AbsoluteTiming
(* {0.00504693, Null} *)

now we can add some styling and print our matrices (the original one with blue entries and the filled matrix with green entries)

rule = Map[# ->  Framed[Style[sol[[Sequence @@ #]], Bold], RoundingRadius -> 50, 
Background -> LightGreen ] &, Position[grid, 0]];

GraphicsRow[
{Graphics[ Map[# /. 0 :>  Framed[Style[0, Bold], RoundingRadius -> 50, Background -> LightBlue] &, grid] //
MatrixForm // Rasterize[#, "Graphics", RasterSize -> 1000] &, ImageSize -> Full], 
Graphics[ReplacePart[sol, Dispatch[rule]] // MatrixForm // Rasterize[#, "Graphics", RasterSize -> 1000] &, 
   ImageSize -> Full]}]

Using Integer Linear Programming

we can deviate from recursion and notice that our system is just a set of equations with integer solutions that satisfy some constraints. Therefore, FindInstance can be used to find a unique solution to our Sudoku (Disclaimer: in contrast to the recursive solution, our initial grid now needs to be undetermined to a lesser extent).

Remaining true to our words, we take the above grid with fewer unknowns than before

grid = {{9, 0, 0, 2, 3, 7, 6, 8, 0}, {0, 2, 0, 8, 4, 0, 0, 7, 3}, {8, 
    0, 7, 1, 0, 5, 0, 2, 9}, {0, 0, 4, 5, 9, 8, 3, 0, 0}, {2, 0, 0, 4,
     0, 1, 0, 0, 6}, {5, 1, 0, 0, 0, 3, 0, 4, 7}, {4, 0, 1, 3, 0, 6, 
    2, 9, 5}, {0, 5, 0, 9, 1, 0, 7, 3, 8}, {3, 0, 8, 0, 5, 0, 1, 0, 
    4}};

symgrid = grid /. 0 :> Unique@x;

symbols = Cases[symgrid, _Symbol, Infinity]; (* symbols in the grid *)

After replacing the unknowns with symbols we proceed to setting up the solution. The necessary conditions require the individual rows and columns to add to 45 and the product of their constituent members to be equal to 9!. Moreover, we need to ensure that all 3x3 submatrices also add to 45.

the equations for submatrices can be generated using the getSubMatrices function

Clear@getSubMatrices;
getSubMatrices[grid_, x_Symbol] := 
 Block[{row, col, rownum, colnum, submatrix, safemem},
  {row, col} = First@Position[grid, x];
  submatrix = (Plus @@ # == 45) &@
    Flatten[Which[row <= 3 && col <= 3, Take[grid, 3, 3],
      row <= 3 && col <= 6, Take[grid, 3, 4 ;; 6], 
      row <= 3 && col <= 9, Take[grid, 3, 7 ;; 9],
      row <= 6 && col <= 3, Take[grid, 4 ;; 6, 3], 
      row <= 6 && col <= 6, Take[grid, 4 ;; 6, 4 ;; 6], 
      row <= 6 && col <= 9, Take[grid, 4 ;; 6, 7 ;; 9],
      row <= 9 && col <= 3, Take[grid, 7 ;; 9, 3], 
      row <= 9 && col <= 6, Take[grid, 7 ;; 9, 4 ;; 6], 
      row <= 9 && col <= 9, Take[grid, 7 ;; 9, 7 ;; 9]]]
] 

the conditions are set up as shown below:

sol = FindInstance[And @@ Join[Table[Plus @@ symgrid[[i]] == 45, {i, 9}], 
Table[Plus @@ symgrid[[All, i]] == 45, {i, 9}], 
Table[Times @@ symgrid[[i]] == 9!, {i, 9}], 
Table[Times @@ symgrid[[All, i]] == 9!, {i, 9}], 
DeleteDuplicates[getSubMatrices[symgrid, #] & /@ symbols]], symbols]

(*  {{x$298 -> 4, x$299 -> 5, x$300 -> 1, x$301 -> 1, x$302 -> 6, x$303 -> 9, x$304 -> 5, x$305 -> 3, x$306 -> 6, x$307 -> 4, 
  x$308 -> 7, x$309 -> 6, x$310 -> 1, x$311 -> 2, x$312 -> 8, x$313 -> 3, x$314 -> 7, x$315 -> 9, x$316 -> 5, x$317 -> 9, 
  x$318 -> 6, x$319 -> 2, x$320 -> 8, x$321 -> 7, x$322 -> 8, x$323 -> 6, x$324 -> 2, x$325 -> 4, x$326 -> 9, x$327 -> 7, 
  x$328 -> 2, x$329 -> 6}}*)

now we can display our solution as a matrix

sympos = Flatten[Position[symgrid, #] & /@ symbols, 1];
symfilled = First[symgrid /. sol];
rule = Map[# ->  Framed[Style[symfilled[[Sequence @@ #]], Bold],  RoundingRadius -> 50, Background -> LightGreen ] &, sympos];
ReplacePart[symgrid, Dispatch[rule]] // MatrixForm // Rasterize[#, "Graphics", RasterSize -> 1000] &

Note: since i had to switch some "unknowns" to "known" for the case where i implemented ILP, recursion is not that bad at solving some problems

Attachments: