Just to add another computer-based solution, then using the MiniZinc modelling language you can write the following program:
int: n;
array[1..n, 1..n] of 0..n: initial_grid;
int: reg;
array[1..n, 1..n] of 1..reg: regions;
array[1..n, 1..n] of var 1..n: final_grid;
include "alldifferent.mzn";
constraint forall(r, c in 1..n)(initial_grid[r, c] = 0 \/ initial_grid[r, c] = final_grid[r, c]);
constraint forall(r in 1..n)(alldifferent([ final_grid[r, c] | c in 1..n ]));
constraint forall(c in 1..n)(alldifferent([ final_grid[r, c] | r in 1..n ]));
constraint forall(region in 1..reg)(alldifferent([ final_grid[r, c] | r, c in 1..n where regions[r, c] = region ]));
solve satisfy;
output [ show_int(1, final_grid[r, c]) ++
if c = n then
("\n"
++ if (r mod 3 = 0 /\ r < n) then "---------------------\n" else "" endif
)
elseif c mod 3 = 0 then " | "
else " "
endif
| r, c in 1..n ];
Along with the appropriate data file:
n = 9;
reg = 9;
regions = array2d(1..9, 1..9, [ 3 * (row div 3) + col div 3 + 1 | row, col in 0..8 ]);
initial_grid =
[| 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 |]
;
And using the default solver on a fairly standard laptop the solution comes out in 100ms, which does beat PM Lee's C++ implementation by a considerable margin.