1## Copyright (C) 2009 Jaroslav Hajek <highegg@gmail.com>
2##
3## This program is free software; you can redistribute it and/or modify it under
4## the terms of the GNU General Public License as published by the Free Software
5## Foundation; either version 3 of the License, or (at your option) any later
6## version.
7##
8## This program is distributed in the hope that it will be useful, but WITHOUT
9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
11## details.
12##
13## You should have received a copy of the GNU General Public License along with
14## this program; if not, see <http://www.gnu.org/licenses/>.
15
16## -*- texinfo -*-
17## @deftypefn{Function File} [@var{x}, @var{ntrial}] = solvesudoku (@var{s})
18## Solves a classical 9x9 sudoku. @var{s} should be a 9x9 array with
19## numbers from 0:9. 0 indicates empty field.
20## Returns the filled table or empty matrix if no solution exists.
21## If requested, @var{ntrial} returns the number of trial-and-error steps needed.
22## @end deftypefn
23
24## This uses a recursive backtracking technique combined with revealing new singleton
25## fields by logic. The beauty of it is that it is completely vectorized.
26
27function [x, ntrial] = solvesudoku (s)
28
29  if (nargin != 1)
30    print_usage ();
31  endif
32
33  if (! (ismatrix (s) && ndims (s) == 2 && all (size (s) == [9, 9])))
34    error ("needs a 9x9 matrix");
35  endif
36
37  if (! ismember (unique (s(:)), 0:9))
38    error ("matrix must contain values from 0:9");
39  endif
40
41  if (! verifysudoku (s))
42    error ("matrix is not a valid sudoku grid");
43  endif
44
45  [x, ntrial] = solvesudoku_rec (s);
46
47endfunction
48
49function ok = verifysudoku (s)
50  [i, j, k] = find (s);
51  b = false (9, 9, 9);
52  b(sub2ind ([9, 9, 9], i, j, k)) = true;
53  okc = sum (b, 1) <= 1;
54  okr = sum (b, 2) <= 1;
55  b = reshape (b, [3, 3, 3, 3, 9]);
56  ok3 = sum (sum (b, 1), 3) <= 1;
57  ok = all (okc(:) & okr(:) & ok3(:));
58endfunction
59
60function [x, ntrial] = solvesudoku_rec (s)
61
62  x = s;
63  ntrial = 0;
64
65  ## Run until the logic is exhausted.
66  do
67    b = getoptions (x);
68    s = x;
69    x = getsingletons (b, x);
70    finished = isempty (x) || all (x(:));
71  until (finished || all ((x == s)(:)));
72
73  if (! finished)
74    x = [];
75    ## Find the field with minimum possibilities.
76    sb = sum (b, 3);
77    sb(s != 0) = 10;
78    [msb, i] = min (sb(:));
79    [i, j] = ind2sub ([9, 9], i);
80    ## Try all guesses.
81    for k = find (b(i,j,:))'
82      s(i,j) = k;
83      [x, ntrial1] = solvesudoku_rec (s);
84      ntrial += 1 + ntrial1;
85      if (! isempty (x))
86        ## Found solutions.
87        break;
88      endif
89      s(i,j) = 0;
90    endfor
91  endif
92
93endfunction
94
95## Given a 9x9x9 logical array of allowed values, get the logical singletons.
96function s = getsingletons (b, s)
97
98  n0 = sum (s(:) != 0);
99
100  ## Check for fields with only one option.
101  sb = sum (b, 3);
102  if (any (sb(:) == 0))
103    s = [];
104    return;
105  else
106    s1 = sb == 1;
107    ## We want to return as soon as some new singletons are found.
108    [s(s1), xx] = find (reshape (b, [], 9)(s1, :).');
109    if (sum (s(:) != 0) > n0)
110      return;
111    endif
112  endif
113
114  ## Check for columns where a number has only one field left.
115  sb = squeeze (sum (b, 1));
116  if (any (sb(:) == 0))
117    s = [];
118    return;
119  else
120    s1 = sb == 1;
121    [j, k] = find (s1);
122    [i, xx] = find (b(:, s1));
123    s(sub2ind ([9, 9], i, j)) = k;
124    if (sum (s(:) != 0) > n0)
125      return;
126    endif
127  endif
128
129  ## Ditto for rows.
130  sb = squeeze (sum (b, 2));
131  if (any (sb(:) == 0))
132    s = [];
133    return;
134  else
135    s1 = sb == 1;
136    [i, k] = find (s1);
137    [j, xx] = find (permute (b, [2, 1, 3])(:, s1));
138    s(sub2ind ([9, 9], i, j)) = k;
139    if (sum (s(:) != 0) > n0)
140      return;
141    endif
142  endif
143
144  ## 3x3 tiles.
145  bb = reshape (b, [3, 3, 3, 3, 9]);
146  sb = squeeze (sum (sum (bb, 1), 3));
147  if (any (sb(:) == 0))
148    s = [];
149    return;
150  else
151    s1 = reshape (sb == 1, 9, 9);
152    [j, k] = find (s1);
153    [i, xx] = find (reshape (permute (bb, [1, 3, 2, 4, 5]), 9, 9*9)(:, s1));
154    [i1, i2] = ind2sub ([3, 3], i);
155    [j1, j2] = ind2sub ([3, 3], j);
156    s(sub2ind ([3, 3, 3, 3], i1, j1, i2, j2)) = k;
157    if (sum (s(:) != 0) > n0)
158      return;
159    endif
160  endif
161
162endfunction
163
164## Given known values (singletons), calculate options.
165function b = getoptions (s)
166
167  ## Find true values.
168  [i, j, s] = find (s);
169  ## Columns.
170  bc = true (9, 9, 9);
171  bc(:, sub2ind ([9, 9], j, s)) = false;
172  ## Rows.
173  br = true (9, 9, 9);
174  br(:, sub2ind ([9, 9], i, s)) = false;
175  ## 3x3 tiles.
176  b3 = true (3, 3, 3, 3, 9);
177  b3(:, :, sub2ind ([3, 3, 9], ceil (i/3), ceil (j/3), s)) = false;
178  ## Permute elements to correct order.
179  br = permute (br, [2, 1, 3]);
180  b3 = reshape (permute (b3, [1, 3, 2, 4, 5]), [9, 9, 9]);
181  ## The singleton fields themselves.
182  bb = true (9*9, 9);
183  bb(sub2ind ([9, 9], i, j), :) = false;
184  bb = reshape (bb, [9, 9, 9]);
185  ## Form result.
186  b = bc & br & b3 & bb;
187  ## Correct singleton fields.
188  b = reshape (b, 9, 9, 9);
189  b(sub2ind ([9, 9, 9], i, j, s)) = true;
190
191endfunction
192
193