149255d7 |
1 | #include <assert.h> |
2 | #include <string.h> |
3 | #include <stdarg.h> |
4 | |
5 | #include "puzzles.h" |
6 | #include "tree234.h" |
7 | #include "maxflow.h" |
8 | |
9 | #ifdef STANDALONE_LATIN_TEST |
10 | #define STANDALONE_SOLVER |
11 | #endif |
12 | |
13 | #include "latin.h" |
14 | |
15 | /* -------------------------------------------------------- |
16 | * Solver. |
17 | */ |
18 | |
62f8425e |
19 | #ifdef STANDALONE_SOLVER |
20 | int solver_show_working, solver_recurse_depth; |
21 | #endif |
22 | |
149255d7 |
23 | /* |
24 | * Function called when we are certain that a particular square has |
25 | * a particular number in it. The y-coordinate passed in here is |
26 | * transformed. |
27 | */ |
28 | void latin_solver_place(struct latin_solver *solver, int x, int y, int n) |
29 | { |
30 | int i, o = solver->o; |
31 | |
32 | assert(n <= o); |
242a7d91 |
33 | assert(cube(x,y,n)); |
149255d7 |
34 | |
35 | /* |
36 | * Rule out all other numbers in this square. |
37 | */ |
38 | for (i = 1; i <= o; i++) |
39 | if (i != n) |
40 | cube(x,y,i) = FALSE; |
41 | |
42 | /* |
43 | * Rule out this number in all other positions in the row. |
44 | */ |
45 | for (i = 0; i < o; i++) |
46 | if (i != y) |
47 | cube(x,i,n) = FALSE; |
48 | |
49 | /* |
50 | * Rule out this number in all other positions in the column. |
51 | */ |
52 | for (i = 0; i < o; i++) |
53 | if (i != x) |
54 | cube(i,y,n) = FALSE; |
55 | |
56 | /* |
57 | * Enter the number in the result grid. |
58 | */ |
59 | solver->grid[YUNTRANS(y)*o+x] = n; |
60 | |
61 | /* |
62 | * Cross out this number from the list of numbers left to place |
63 | * in its row, its column and its block. |
64 | */ |
65 | solver->row[y*o+n-1] = solver->col[x*o+n-1] = TRUE; |
66 | } |
67 | |
68 | int latin_solver_elim(struct latin_solver *solver, int start, int step |
69 | #ifdef STANDALONE_SOLVER |
70 | , char *fmt, ... |
71 | #endif |
72 | ) |
73 | { |
74 | int o = solver->o; |
75 | int fpos, m, i; |
76 | |
77 | /* |
78 | * Count the number of set bits within this section of the |
79 | * cube. |
80 | */ |
81 | m = 0; |
82 | fpos = -1; |
83 | for (i = 0; i < o; i++) |
84 | if (solver->cube[start+i*step]) { |
85 | fpos = start+i*step; |
86 | m++; |
87 | } |
88 | |
89 | if (m == 1) { |
90 | int x, y, n; |
91 | assert(fpos >= 0); |
92 | |
93 | n = 1 + fpos % o; |
94 | y = fpos / o; |
95 | x = y / o; |
96 | y %= o; |
97 | |
98 | if (!solver->grid[YUNTRANS(y)*o+x]) { |
99 | #ifdef STANDALONE_SOLVER |
100 | if (solver_show_working) { |
101 | va_list ap; |
102 | printf("%*s", solver_recurse_depth*4, ""); |
103 | va_start(ap, fmt); |
104 | vprintf(fmt, ap); |
105 | va_end(ap); |
106 | printf(":\n%*s placing %d at (%d,%d)\n", |
6a7d52aa |
107 | solver_recurse_depth*4, "", n, x+1, YUNTRANS(y)+1); |
149255d7 |
108 | } |
109 | #endif |
110 | latin_solver_place(solver, x, y, n); |
111 | return +1; |
112 | } |
113 | } else if (m == 0) { |
114 | #ifdef STANDALONE_SOLVER |
115 | if (solver_show_working) { |
116 | va_list ap; |
117 | printf("%*s", solver_recurse_depth*4, ""); |
118 | va_start(ap, fmt); |
119 | vprintf(fmt, ap); |
120 | va_end(ap); |
121 | printf(":\n%*s no possibilities available\n", |
122 | solver_recurse_depth*4, ""); |
123 | } |
124 | #endif |
125 | return -1; |
126 | } |
127 | |
128 | return 0; |
129 | } |
130 | |
131 | struct latin_solver_scratch { |
132 | unsigned char *grid, *rowidx, *colidx, *set; |
133 | int *neighbours, *bfsqueue; |
134 | #ifdef STANDALONE_SOLVER |
135 | int *bfsprev; |
136 | #endif |
137 | }; |
138 | |
139 | int latin_solver_set(struct latin_solver *solver, |
140 | struct latin_solver_scratch *scratch, |
141 | int start, int step1, int step2 |
142 | #ifdef STANDALONE_SOLVER |
143 | , char *fmt, ... |
144 | #endif |
145 | ) |
146 | { |
147 | int o = solver->o; |
148 | int i, j, n, count; |
149 | unsigned char *grid = scratch->grid; |
150 | unsigned char *rowidx = scratch->rowidx; |
151 | unsigned char *colidx = scratch->colidx; |
152 | unsigned char *set = scratch->set; |
153 | |
154 | /* |
155 | * We are passed a o-by-o matrix of booleans. Our first job |
156 | * is to winnow it by finding any definite placements - i.e. |
157 | * any row with a solitary 1 - and discarding that row and the |
158 | * column containing the 1. |
159 | */ |
160 | memset(rowidx, TRUE, o); |
161 | memset(colidx, TRUE, o); |
162 | for (i = 0; i < o; i++) { |
163 | int count = 0, first = -1; |
164 | for (j = 0; j < o; j++) |
165 | if (solver->cube[start+i*step1+j*step2]) |
166 | first = j, count++; |
167 | |
168 | if (count == 0) return -1; |
169 | if (count == 1) |
170 | rowidx[i] = colidx[first] = FALSE; |
171 | } |
172 | |
173 | /* |
174 | * Convert each of rowidx/colidx from a list of 0s and 1s to a |
175 | * list of the indices of the 1s. |
176 | */ |
177 | for (i = j = 0; i < o; i++) |
178 | if (rowidx[i]) |
179 | rowidx[j++] = i; |
180 | n = j; |
181 | for (i = j = 0; i < o; i++) |
182 | if (colidx[i]) |
183 | colidx[j++] = i; |
184 | assert(n == j); |
185 | |
186 | /* |
187 | * And create the smaller matrix. |
188 | */ |
189 | for (i = 0; i < n; i++) |
190 | for (j = 0; j < n; j++) |
191 | grid[i*o+j] = solver->cube[start+rowidx[i]*step1+colidx[j]*step2]; |
192 | |
193 | /* |
194 | * Having done that, we now have a matrix in which every row |
195 | * has at least two 1s in. Now we search to see if we can find |
196 | * a rectangle of zeroes (in the set-theoretic sense of |
197 | * `rectangle', i.e. a subset of rows crossed with a subset of |
198 | * columns) whose width and height add up to n. |
199 | */ |
200 | |
201 | memset(set, 0, n); |
202 | count = 0; |
203 | while (1) { |
204 | /* |
205 | * We have a candidate set. If its size is <=1 or >=n-1 |
206 | * then we move on immediately. |
207 | */ |
208 | if (count > 1 && count < n-1) { |
209 | /* |
210 | * The number of rows we need is n-count. See if we can |
211 | * find that many rows which each have a zero in all |
212 | * the positions listed in `set'. |
213 | */ |
214 | int rows = 0; |
215 | for (i = 0; i < n; i++) { |
216 | int ok = TRUE; |
217 | for (j = 0; j < n; j++) |
218 | if (set[j] && grid[i*o+j]) { |
219 | ok = FALSE; |
220 | break; |
221 | } |
222 | if (ok) |
223 | rows++; |
224 | } |
225 | |
226 | /* |
227 | * We expect never to be able to get _more_ than |
228 | * n-count suitable rows: this would imply that (for |
229 | * example) there are four numbers which between them |
230 | * have at most three possible positions, and hence it |
231 | * indicates a faulty deduction before this point or |
232 | * even a bogus clue. |
233 | */ |
234 | if (rows > n - count) { |
235 | #ifdef STANDALONE_SOLVER |
236 | if (solver_show_working) { |
237 | va_list ap; |
238 | printf("%*s", solver_recurse_depth*4, |
239 | ""); |
240 | va_start(ap, fmt); |
241 | vprintf(fmt, ap); |
242 | va_end(ap); |
243 | printf(":\n%*s contradiction reached\n", |
244 | solver_recurse_depth*4, ""); |
245 | } |
246 | #endif |
247 | return -1; |
248 | } |
249 | |
250 | if (rows >= n - count) { |
251 | int progress = FALSE; |
252 | |
253 | /* |
254 | * We've got one! Now, for each row which _doesn't_ |
255 | * satisfy the criterion, eliminate all its set |
256 | * bits in the positions _not_ listed in `set'. |
257 | * Return +1 (meaning progress has been made) if we |
258 | * successfully eliminated anything at all. |
259 | * |
260 | * This involves referring back through |
261 | * rowidx/colidx in order to work out which actual |
262 | * positions in the cube to meddle with. |
263 | */ |
264 | for (i = 0; i < n; i++) { |
265 | int ok = TRUE; |
266 | for (j = 0; j < n; j++) |
267 | if (set[j] && grid[i*o+j]) { |
268 | ok = FALSE; |
269 | break; |
270 | } |
271 | if (!ok) { |
272 | for (j = 0; j < n; j++) |
273 | if (!set[j] && grid[i*o+j]) { |
274 | int fpos = (start+rowidx[i]*step1+ |
275 | colidx[j]*step2); |
276 | #ifdef STANDALONE_SOLVER |
277 | if (solver_show_working) { |
278 | int px, py, pn; |
279 | |
280 | if (!progress) { |
281 | va_list ap; |
282 | printf("%*s", solver_recurse_depth*4, |
283 | ""); |
284 | va_start(ap, fmt); |
285 | vprintf(fmt, ap); |
286 | va_end(ap); |
287 | printf(":\n"); |
288 | } |
289 | |
290 | pn = 1 + fpos % o; |
291 | py = fpos / o; |
292 | px = py / o; |
293 | py %= o; |
294 | |
295 | printf("%*s ruling out %d at (%d,%d)\n", |
296 | solver_recurse_depth*4, "", |
6a7d52aa |
297 | pn, px+1, YUNTRANS(py)+1); |
149255d7 |
298 | } |
299 | #endif |
300 | progress = TRUE; |
301 | solver->cube[fpos] = FALSE; |
302 | } |
303 | } |
304 | } |
305 | |
306 | if (progress) { |
307 | return +1; |
308 | } |
309 | } |
310 | } |
311 | |
312 | /* |
313 | * Binary increment: change the rightmost 0 to a 1, and |
314 | * change all 1s to the right of it to 0s. |
315 | */ |
316 | i = n; |
317 | while (i > 0 && set[i-1]) |
318 | set[--i] = 0, count--; |
319 | if (i > 0) |
320 | set[--i] = 1, count++; |
321 | else |
322 | break; /* done */ |
323 | } |
324 | |
325 | return 0; |
326 | } |
327 | |
328 | /* |
329 | * Look for forcing chains. A forcing chain is a path of |
330 | * pairwise-exclusive squares (i.e. each pair of adjacent squares |
331 | * in the path are in the same row, column or block) with the |
332 | * following properties: |
333 | * |
334 | * (a) Each square on the path has precisely two possible numbers. |
335 | * |
336 | * (b) Each pair of squares which are adjacent on the path share |
337 | * at least one possible number in common. |
338 | * |
339 | * (c) Each square in the middle of the path shares _both_ of its |
340 | * numbers with at least one of its neighbours (not the same |
341 | * one with both neighbours). |
342 | * |
343 | * These together imply that at least one of the possible number |
344 | * choices at one end of the path forces _all_ the rest of the |
345 | * numbers along the path. In order to make real use of this, we |
346 | * need further properties: |
347 | * |
348 | * (c) Ruling out some number N from the square at one end |
349 | * of the path forces the square at the other end to |
350 | * take number N. |
351 | * |
352 | * (d) The two end squares are both in line with some third |
353 | * square. |
354 | * |
355 | * (e) That third square currently has N as a possibility. |
356 | * |
357 | * If we can find all of that lot, we can deduce that at least one |
358 | * of the two ends of the forcing chain has number N, and that |
359 | * therefore the mutually adjacent third square does not. |
360 | * |
361 | * To find forcing chains, we're going to start a bfs at each |
362 | * suitable square, once for each of its two possible numbers. |
363 | */ |
364 | int latin_solver_forcing(struct latin_solver *solver, |
365 | struct latin_solver_scratch *scratch) |
366 | { |
367 | int o = solver->o; |
368 | int *bfsqueue = scratch->bfsqueue; |
369 | #ifdef STANDALONE_SOLVER |
370 | int *bfsprev = scratch->bfsprev; |
371 | #endif |
372 | unsigned char *number = scratch->grid; |
373 | int *neighbours = scratch->neighbours; |
374 | int x, y; |
375 | |
376 | for (y = 0; y < o; y++) |
377 | for (x = 0; x < o; x++) { |
378 | int count, t, n; |
379 | |
380 | /* |
381 | * If this square doesn't have exactly two candidate |
382 | * numbers, don't try it. |
383 | * |
384 | * In this loop we also sum the candidate numbers, |
385 | * which is a nasty hack to allow us to quickly find |
386 | * `the other one' (since we will shortly know there |
387 | * are exactly two). |
388 | */ |
389 | for (count = t = 0, n = 1; n <= o; n++) |
390 | if (cube(x, y, n)) |
391 | count++, t += n; |
392 | if (count != 2) |
393 | continue; |
394 | |
395 | /* |
396 | * Now attempt a bfs for each candidate. |
397 | */ |
398 | for (n = 1; n <= o; n++) |
399 | if (cube(x, y, n)) { |
400 | int orign, currn, head, tail; |
401 | |
402 | /* |
403 | * Begin a bfs. |
404 | */ |
405 | orign = n; |
406 | |
407 | memset(number, o+1, o*o); |
408 | head = tail = 0; |
409 | bfsqueue[tail++] = y*o+x; |
410 | #ifdef STANDALONE_SOLVER |
411 | bfsprev[y*o+x] = -1; |
412 | #endif |
413 | number[y*o+x] = t - n; |
414 | |
415 | while (head < tail) { |
416 | int xx, yy, nneighbours, xt, yt, i; |
417 | |
418 | xx = bfsqueue[head++]; |
419 | yy = xx / o; |
420 | xx %= o; |
421 | |
422 | currn = number[yy*o+xx]; |
423 | |
424 | /* |
425 | * Find neighbours of yy,xx. |
426 | */ |
427 | nneighbours = 0; |
428 | for (yt = 0; yt < o; yt++) |
429 | neighbours[nneighbours++] = yt*o+xx; |
430 | for (xt = 0; xt < o; xt++) |
431 | neighbours[nneighbours++] = yy*o+xt; |
432 | |
433 | /* |
434 | * Try visiting each of those neighbours. |
435 | */ |
436 | for (i = 0; i < nneighbours; i++) { |
437 | int cc, tt, nn; |
438 | |
439 | xt = neighbours[i] % o; |
440 | yt = neighbours[i] / o; |
441 | |
442 | /* |
443 | * We need this square to not be |
444 | * already visited, and to include |
445 | * currn as a possible number. |
446 | */ |
447 | if (number[yt*o+xt] <= o) |
448 | continue; |
449 | if (!cube(xt, yt, currn)) |
450 | continue; |
451 | |
452 | /* |
453 | * Don't visit _this_ square a second |
454 | * time! |
455 | */ |
456 | if (xt == xx && yt == yy) |
457 | continue; |
458 | |
459 | /* |
460 | * To continue with the bfs, we need |
461 | * this square to have exactly two |
462 | * possible numbers. |
463 | */ |
464 | for (cc = tt = 0, nn = 1; nn <= o; nn++) |
465 | if (cube(xt, yt, nn)) |
466 | cc++, tt += nn; |
467 | if (cc == 2) { |
468 | bfsqueue[tail++] = yt*o+xt; |
469 | #ifdef STANDALONE_SOLVER |
470 | bfsprev[yt*o+xt] = yy*o+xx; |
471 | #endif |
472 | number[yt*o+xt] = tt - currn; |
473 | } |
474 | |
475 | /* |
476 | * One other possibility is that this |
477 | * might be the square in which we can |
478 | * make a real deduction: if it's |
479 | * adjacent to x,y, and currn is equal |
480 | * to the original number we ruled out. |
481 | */ |
482 | if (currn == orign && |
483 | (xt == x || yt == y)) { |
484 | #ifdef STANDALONE_SOLVER |
485 | if (solver_show_working) { |
486 | char *sep = ""; |
487 | int xl, yl; |
488 | printf("%*sforcing chain, %d at ends of ", |
489 | solver_recurse_depth*4, "", orign); |
490 | xl = xx; |
491 | yl = yy; |
492 | while (1) { |
6a7d52aa |
493 | printf("%s(%d,%d)", sep, xl+1, |
494 | YUNTRANS(yl)+1); |
149255d7 |
495 | xl = bfsprev[yl*o+xl]; |
496 | if (xl < 0) |
497 | break; |
498 | yl = xl / o; |
499 | xl %= o; |
500 | sep = "-"; |
501 | } |
502 | printf("\n%*s ruling out %d at (%d,%d)\n", |
503 | solver_recurse_depth*4, "", |
6a7d52aa |
504 | orign, xt+1, YUNTRANS(yt)+1); |
149255d7 |
505 | } |
506 | #endif |
507 | cube(xt, yt, orign) = FALSE; |
508 | return 1; |
509 | } |
510 | } |
511 | } |
512 | } |
513 | } |
514 | |
515 | return 0; |
516 | } |
517 | |
518 | struct latin_solver_scratch *latin_solver_new_scratch(struct latin_solver *solver) |
519 | { |
520 | struct latin_solver_scratch *scratch = snew(struct latin_solver_scratch); |
521 | int o = solver->o; |
522 | scratch->grid = snewn(o*o, unsigned char); |
523 | scratch->rowidx = snewn(o, unsigned char); |
524 | scratch->colidx = snewn(o, unsigned char); |
525 | scratch->set = snewn(o, unsigned char); |
526 | scratch->neighbours = snewn(3*o, int); |
527 | scratch->bfsqueue = snewn(o*o, int); |
528 | #ifdef STANDALONE_SOLVER |
529 | scratch->bfsprev = snewn(o*o, int); |
530 | #endif |
531 | return scratch; |
532 | } |
533 | |
534 | void latin_solver_free_scratch(struct latin_solver_scratch *scratch) |
535 | { |
536 | #ifdef STANDALONE_SOLVER |
537 | sfree(scratch->bfsprev); |
538 | #endif |
539 | sfree(scratch->bfsqueue); |
540 | sfree(scratch->neighbours); |
541 | sfree(scratch->set); |
542 | sfree(scratch->colidx); |
543 | sfree(scratch->rowidx); |
544 | sfree(scratch->grid); |
545 | sfree(scratch); |
546 | } |
547 | |
548 | void latin_solver_alloc(struct latin_solver *solver, digit *grid, int o) |
549 | { |
550 | int x, y; |
551 | |
552 | solver->o = o; |
553 | solver->cube = snewn(o*o*o, unsigned char); |
554 | solver->grid = grid; /* write straight back to the input */ |
555 | memset(solver->cube, TRUE, o*o*o); |
556 | |
557 | solver->row = snewn(o*o, unsigned char); |
558 | solver->col = snewn(o*o, unsigned char); |
559 | memset(solver->row, FALSE, o*o); |
560 | memset(solver->col, FALSE, o*o); |
561 | |
562 | for (x = 0; x < o; x++) |
563 | for (y = 0; y < o; y++) |
564 | if (grid[y*o+x]) |
565 | latin_solver_place(solver, x, YTRANS(y), grid[y*o+x]); |
566 | } |
567 | |
568 | void latin_solver_free(struct latin_solver *solver) |
569 | { |
570 | sfree(solver->cube); |
571 | sfree(solver->row); |
572 | sfree(solver->col); |
573 | } |
574 | |
575 | int latin_solver_diff_simple(struct latin_solver *solver) |
576 | { |
577 | int x, y, n, ret, o = solver->o; |
578 | /* |
579 | * Row-wise positional elimination. |
580 | */ |
581 | for (y = 0; y < o; y++) |
582 | for (n = 1; n <= o; n++) |
583 | if (!solver->row[y*o+n-1]) { |
584 | ret = latin_solver_elim(solver, cubepos(0,y,n), o*o |
585 | #ifdef STANDALONE_SOLVER |
586 | , "positional elimination," |
6a7d52aa |
587 | " %d in row %d", n, YUNTRANS(y)+1 |
149255d7 |
588 | #endif |
589 | ); |
590 | if (ret != 0) return ret; |
591 | } |
592 | /* |
593 | * Column-wise positional elimination. |
594 | */ |
595 | for (x = 0; x < o; x++) |
596 | for (n = 1; n <= o; n++) |
597 | if (!solver->col[x*o+n-1]) { |
598 | ret = latin_solver_elim(solver, cubepos(x,0,n), o |
599 | #ifdef STANDALONE_SOLVER |
600 | , "positional elimination," |
6a7d52aa |
601 | " %d in column %d", n, x+1 |
149255d7 |
602 | #endif |
603 | ); |
604 | if (ret != 0) return ret; |
605 | } |
606 | |
607 | /* |
608 | * Numeric elimination. |
609 | */ |
610 | for (x = 0; x < o; x++) |
611 | for (y = 0; y < o; y++) |
612 | if (!solver->grid[YUNTRANS(y)*o+x]) { |
613 | ret = latin_solver_elim(solver, cubepos(x,y,1), 1 |
614 | #ifdef STANDALONE_SOLVER |
6a7d52aa |
615 | , "numeric elimination at (%d,%d)", |
616 | x+1, YUNTRANS(y)+1 |
149255d7 |
617 | #endif |
618 | ); |
619 | if (ret != 0) return ret; |
620 | } |
621 | return 0; |
622 | } |
623 | |
624 | int latin_solver_diff_set(struct latin_solver *solver, |
625 | struct latin_solver_scratch *scratch, |
481628b3 |
626 | int extreme) |
149255d7 |
627 | { |
628 | int x, y, n, ret, o = solver->o; |
481628b3 |
629 | |
630 | if (!extreme) { |
631 | /* |
632 | * Row-wise set elimination. |
633 | */ |
634 | for (y = 0; y < o; y++) { |
635 | ret = latin_solver_set(solver, scratch, cubepos(0,y,1), o*o, 1 |
149255d7 |
636 | #ifdef STANDALONE_SOLVER |
6a7d52aa |
637 | , "set elimination, row %d", YUNTRANS(y)+1 |
149255d7 |
638 | #endif |
481628b3 |
639 | ); |
640 | if (ret != 0) return ret; |
641 | } |
642 | /* |
643 | * Column-wise set elimination. |
644 | */ |
645 | for (x = 0; x < o; x++) { |
646 | ret = latin_solver_set(solver, scratch, cubepos(x,0,1), o, 1 |
149255d7 |
647 | #ifdef STANDALONE_SOLVER |
6a7d52aa |
648 | , "set elimination, column %d", x+1 |
149255d7 |
649 | #endif |
481628b3 |
650 | ); |
651 | if (ret != 0) return ret; |
652 | } |
653 | } else { |
654 | /* |
655 | * Row-vs-column set elimination on a single number |
656 | * (much tricker for a human to do!) |
657 | */ |
658 | for (n = 1; n <= o; n++) { |
659 | ret = latin_solver_set(solver, scratch, cubepos(0,0,n), o*o, o |
149255d7 |
660 | #ifdef STANDALONE_SOLVER |
481628b3 |
661 | , "positional set elimination, number %d", n |
149255d7 |
662 | #endif |
481628b3 |
663 | ); |
664 | if (ret != 0) return ret; |
665 | } |
149255d7 |
666 | } |
667 | return 0; |
668 | } |
669 | |
388b2f00 |
670 | /* |
671 | * Returns: |
149255d7 |
672 | * 0 for 'didn't do anything' implying it was already solved. |
673 | * -1 for 'impossible' (no solution) |
674 | * 1 for 'single solution' |
675 | * >1 for 'multiple solutions' (you don't get to know how many, and |
676 | * the first such solution found will be set. |
677 | * |
678 | * and this function may well assert if given an impossible board. |
679 | */ |
388b2f00 |
680 | static int latin_solver_recurse |
681 | (struct latin_solver *solver, int diff_simple, int diff_set_0, |
682 | int diff_set_1, int diff_forcing, int diff_recursive, |
683 | usersolver_t const *usersolvers, void *ctx, |
684 | ctxnew_t ctxnew, ctxfree_t ctxfree) |
149255d7 |
685 | { |
686 | int best, bestcount; |
687 | int o = solver->o, x, y, n; |
688 | |
689 | best = -1; |
690 | bestcount = o+1; |
691 | |
692 | for (y = 0; y < o; y++) |
693 | for (x = 0; x < o; x++) |
694 | if (!solver->grid[y*o+x]) { |
695 | int count; |
696 | |
697 | /* |
698 | * An unfilled square. Count the number of |
699 | * possible digits in it. |
700 | */ |
701 | count = 0; |
702 | for (n = 1; n <= o; n++) |
703 | if (cube(x,YTRANS(y),n)) |
704 | count++; |
705 | |
706 | /* |
707 | * We should have found any impossibilities |
708 | * already, so this can safely be an assert. |
709 | */ |
710 | assert(count > 1); |
711 | |
712 | if (count < bestcount) { |
713 | bestcount = count; |
714 | best = y*o+x; |
715 | } |
716 | } |
717 | |
718 | if (best == -1) |
719 | /* we were complete already. */ |
720 | return 0; |
721 | else { |
722 | int i, j; |
723 | digit *list, *ingrid, *outgrid; |
724 | int diff = diff_impossible; /* no solution found yet */ |
725 | |
726 | /* |
727 | * Attempt recursion. |
728 | */ |
729 | y = best / o; |
730 | x = best % o; |
731 | |
732 | list = snewn(o, digit); |
733 | ingrid = snewn(o*o, digit); |
734 | outgrid = snewn(o*o, digit); |
735 | memcpy(ingrid, solver->grid, o*o); |
736 | |
737 | /* Make a list of the possible digits. */ |
738 | for (j = 0, n = 1; n <= o; n++) |
739 | if (cube(x,YTRANS(y),n)) |
740 | list[j++] = n; |
741 | |
742 | #ifdef STANDALONE_SOLVER |
743 | if (solver_show_working) { |
744 | char *sep = ""; |
745 | printf("%*srecursing on (%d,%d) [", |
6a7d52aa |
746 | solver_recurse_depth*4, "", x+1, y+1); |
149255d7 |
747 | for (i = 0; i < j; i++) { |
748 | printf("%s%d", sep, list[i]); |
749 | sep = " or "; |
750 | } |
751 | printf("]\n"); |
752 | } |
753 | #endif |
754 | |
755 | /* |
756 | * And step along the list, recursing back into the |
757 | * main solver at every stage. |
758 | */ |
759 | for (i = 0; i < j; i++) { |
760 | int ret; |
388b2f00 |
761 | void *newctx; |
149255d7 |
762 | |
763 | memcpy(outgrid, ingrid, o*o); |
764 | outgrid[y*o+x] = list[i]; |
765 | |
766 | #ifdef STANDALONE_SOLVER |
767 | if (solver_show_working) |
768 | printf("%*sguessing %d at (%d,%d)\n", |
6a7d52aa |
769 | solver_recurse_depth*4, "", list[i], x+1, y+1); |
149255d7 |
770 | solver_recurse_depth++; |
771 | #endif |
772 | |
388b2f00 |
773 | if (ctxnew) { |
774 | newctx = ctxnew(ctx); |
775 | } else { |
776 | newctx = ctx; |
777 | } |
778 | ret = latin_solver(outgrid, o, diff_recursive, |
779 | diff_simple, diff_set_0, diff_set_1, |
780 | diff_forcing, diff_recursive, |
781 | usersolvers, newctx, ctxnew, ctxfree); |
782 | if (ctxnew) |
783 | ctxfree(newctx); |
149255d7 |
784 | |
785 | #ifdef STANDALONE_SOLVER |
786 | solver_recurse_depth--; |
787 | if (solver_show_working) { |
788 | printf("%*sretracting %d at (%d,%d)\n", |
6a7d52aa |
789 | solver_recurse_depth*4, "", list[i], x+1, y+1); |
149255d7 |
790 | } |
791 | #endif |
792 | /* we recurse as deep as we can, so we should never find |
793 | * find ourselves giving up on a puzzle without declaring it |
794 | * impossible. */ |
795 | assert(ret != diff_unfinished); |
796 | |
797 | /* |
798 | * If we have our first solution, copy it into the |
799 | * grid we will return. |
800 | */ |
801 | if (diff == diff_impossible && ret != diff_impossible) |
802 | memcpy(solver->grid, outgrid, o*o); |
803 | |
804 | if (ret == diff_ambiguous) |
805 | diff = diff_ambiguous; |
806 | else if (ret == diff_impossible) |
807 | /* do not change our return value */; |
808 | else { |
809 | /* the recursion turned up exactly one solution */ |
810 | if (diff == diff_impossible) |
388b2f00 |
811 | diff = diff_recursive; |
149255d7 |
812 | else |
813 | diff = diff_ambiguous; |
814 | } |
815 | |
816 | /* |
817 | * As soon as we've found more than one solution, |
818 | * give up immediately. |
819 | */ |
820 | if (diff == diff_ambiguous) |
821 | break; |
822 | } |
823 | |
824 | sfree(outgrid); |
825 | sfree(ingrid); |
826 | sfree(list); |
827 | |
828 | if (diff == diff_impossible) |
829 | return -1; |
830 | else if (diff == diff_ambiguous) |
831 | return 2; |
832 | else { |
388b2f00 |
833 | assert(diff == diff_recursive); |
149255d7 |
834 | return 1; |
835 | } |
836 | } |
837 | } |
838 | |
388b2f00 |
839 | int latin_solver_main(struct latin_solver *solver, int maxdiff, |
840 | int diff_simple, int diff_set_0, int diff_set_1, |
841 | int diff_forcing, int diff_recursive, |
842 | usersolver_t const *usersolvers, void *ctx, |
843 | ctxnew_t ctxnew, ctxfree_t ctxfree) |
149255d7 |
844 | { |
845 | struct latin_solver_scratch *scratch = latin_solver_new_scratch(solver); |
481628b3 |
846 | int ret, diff = diff_simple; |
149255d7 |
847 | |
848 | assert(maxdiff <= diff_recursive); |
849 | /* |
850 | * Now loop over the grid repeatedly trying all permitted modes |
851 | * of reasoning. The loop terminates if we complete an |
852 | * iteration without making any progress; we then return |
853 | * failure or success depending on whether the grid is full or |
854 | * not. |
855 | */ |
856 | while (1) { |
388b2f00 |
857 | int i; |
149255d7 |
858 | |
388b2f00 |
859 | cont: |
149255d7 |
860 | |
388b2f00 |
861 | latin_solver_debug(solver->cube, solver->o); |
481628b3 |
862 | |
388b2f00 |
863 | for (i = 0; i <= maxdiff; i++) { |
864 | if (usersolvers[i]) |
865 | ret = usersolvers[i](solver, ctx); |
866 | else |
867 | ret = 0; |
868 | if (ret == 0 && i == diff_simple) |
869 | ret = latin_solver_diff_simple(solver); |
870 | if (ret == 0 && i == diff_set_0) |
871 | ret = latin_solver_diff_set(solver, scratch, 0); |
872 | if (ret == 0 && i == diff_set_1) |
873 | ret = latin_solver_diff_set(solver, scratch, 1); |
874 | if (ret == 0 && i == diff_forcing) |
875 | ret = latin_solver_forcing(solver, scratch); |
876 | |
877 | if (ret < 0) { |
878 | diff = diff_impossible; |
879 | goto got_result; |
880 | } else if (ret > 0) { |
881 | diff = max(diff, i); |
882 | goto cont; |
883 | } |
884 | } |
149255d7 |
885 | |
886 | /* |
887 | * If we reach here, we have made no deductions in this |
888 | * iteration, so the algorithm terminates. |
889 | */ |
890 | break; |
891 | } |
892 | |
893 | /* |
894 | * Last chance: if we haven't fully solved the puzzle yet, try |
895 | * recursing based on guesses for a particular square. We pick |
896 | * one of the most constrained empty squares we can find, which |
897 | * has the effect of pruning the search tree as much as |
898 | * possible. |
899 | */ |
900 | if (maxdiff == diff_recursive) { |
388b2f00 |
901 | int nsol = latin_solver_recurse(solver, |
902 | diff_simple, diff_set_0, diff_set_1, |
903 | diff_forcing, diff_recursive, |
904 | usersolvers, ctx, ctxnew, ctxfree); |
149255d7 |
905 | if (nsol < 0) diff = diff_impossible; |
906 | else if (nsol == 1) diff = diff_recursive; |
907 | else if (nsol > 1) diff = diff_ambiguous; |
908 | /* if nsol == 0 then we were complete anyway |
909 | * (and thus don't need to change diff) */ |
910 | } else { |
911 | /* |
912 | * We're forbidden to use recursion, so we just see whether |
913 | * our grid is fully solved, and return diff_unfinished |
914 | * otherwise. |
915 | */ |
916 | int x, y, o = solver->o; |
917 | |
918 | for (y = 0; y < o; y++) |
919 | for (x = 0; x < o; x++) |
920 | if (!solver->grid[y*o+x]) |
921 | diff = diff_unfinished; |
922 | } |
923 | |
924 | got_result: |
925 | |
926 | #ifdef STANDALONE_SOLVER |
927 | if (solver_show_working) |
928 | printf("%*s%s found\n", |
929 | solver_recurse_depth*4, "", |
930 | diff == diff_impossible ? "no solution (impossible)" : |
931 | diff == diff_unfinished ? "no solution (unfinished)" : |
932 | diff == diff_ambiguous ? "multiple solutions" : |
933 | "one solution"); |
934 | #endif |
935 | |
936 | latin_solver_free_scratch(scratch); |
937 | |
938 | return diff; |
939 | } |
940 | |
388b2f00 |
941 | int latin_solver(digit *grid, int o, int maxdiff, |
942 | int diff_simple, int diff_set_0, int diff_set_1, |
943 | int diff_forcing, int diff_recursive, |
944 | usersolver_t const *usersolvers, void *ctx, |
945 | ctxnew_t ctxnew, ctxfree_t ctxfree) |
149255d7 |
946 | { |
947 | struct latin_solver solver; |
948 | int diff; |
949 | |
950 | latin_solver_alloc(&solver, grid, o); |
388b2f00 |
951 | diff = latin_solver_main(&solver, maxdiff, |
952 | diff_simple, diff_set_0, diff_set_1, |
953 | diff_forcing, diff_recursive, |
954 | usersolvers, ctx, ctxnew, ctxfree); |
149255d7 |
955 | latin_solver_free(&solver); |
956 | return diff; |
957 | } |
958 | |
959 | void latin_solver_debug(unsigned char *cube, int o) |
960 | { |
961 | #ifdef STANDALONE_SOLVER |
388b2f00 |
962 | if (solver_show_working > 1) { |
149255d7 |
963 | struct latin_solver ls, *solver = &ls; |
242a7d91 |
964 | char *dbg; |
149255d7 |
965 | int x, y, i, c = 0; |
966 | |
967 | ls.cube = cube; ls.o = o; /* for cube() to work */ |
968 | |
3a271e7c |
969 | dbg = snewn(3*o*o*o, char); |
149255d7 |
970 | for (y = 0; y < o; y++) { |
971 | for (x = 0; x < o; x++) { |
972 | for (i = 1; i <= o; i++) { |
973 | if (cube(x,y,i)) |
974 | dbg[c++] = i + '0'; |
975 | else |
976 | dbg[c++] = '.'; |
977 | } |
978 | dbg[c++] = ' '; |
979 | } |
980 | dbg[c++] = '\n'; |
981 | } |
982 | dbg[c++] = '\n'; |
983 | dbg[c++] = '\0'; |
984 | |
985 | printf("%s", dbg); |
986 | sfree(dbg); |
987 | } |
988 | #endif |
989 | } |
990 | |
991 | void latin_debug(digit *sq, int o) |
992 | { |
993 | #ifdef STANDALONE_SOLVER |
994 | if (solver_show_working) { |
995 | int x, y; |
996 | |
997 | for (y = 0; y < o; y++) { |
998 | for (x = 0; x < o; x++) { |
999 | printf("%2d ", sq[y*o+x]); |
1000 | } |
1001 | printf("\n"); |
1002 | } |
1003 | printf("\n"); |
1004 | } |
1005 | #endif |
1006 | } |
1007 | |
1008 | /* -------------------------------------------------------- |
1009 | * Generation. |
1010 | */ |
1011 | |
1012 | digit *latin_generate(int o, random_state *rs) |
1013 | { |
1014 | digit *sq; |
1015 | int *edges, *backedges, *capacity, *flow; |
1016 | void *scratch; |
1017 | int ne, scratchsize; |
1018 | int i, j, k; |
1019 | digit *row, *col, *numinv, *num; |
1020 | |
1021 | /* |
1022 | * To efficiently generate a latin square in such a way that |
1023 | * all possible squares are possible outputs from the function, |
1024 | * we make use of a theorem which states that any r x n latin |
1025 | * rectangle, with r < n, can be extended into an (r+1) x n |
1026 | * latin rectangle. In other words, we can reliably generate a |
1027 | * latin square row by row, by at every stage writing down any |
1028 | * row at all which doesn't conflict with previous rows, and |
1029 | * the theorem guarantees that we will never have to backtrack. |
1030 | * |
1031 | * To find a viable row at each stage, we can make use of the |
1032 | * support functions in maxflow.c. |
1033 | */ |
1034 | |
1035 | sq = snewn(o*o, digit); |
1036 | |
1037 | /* |
1038 | * In case this method of generation introduces a really subtle |
1039 | * top-to-bottom directional bias, we'll generate the rows in |
1040 | * random order. |
1041 | */ |
1042 | row = snewn(o, digit); |
1043 | col = snewn(o, digit); |
1044 | numinv = snewn(o, digit); |
1045 | num = snewn(o, digit); |
1046 | for (i = 0; i < o; i++) |
1047 | row[i] = i; |
1048 | shuffle(row, i, sizeof(*row), rs); |
1049 | |
1050 | /* |
1051 | * Set up the infrastructure for the maxflow algorithm. |
1052 | */ |
1053 | scratchsize = maxflow_scratch_size(o * 2 + 2); |
1054 | scratch = smalloc(scratchsize); |
1055 | backedges = snewn(o*o + 2*o, int); |
1056 | edges = snewn((o*o + 2*o) * 2, int); |
1057 | capacity = snewn(o*o + 2*o, int); |
1058 | flow = snewn(o*o + 2*o, int); |
1059 | /* Set up the edge array, and the initial capacities. */ |
1060 | ne = 0; |
1061 | for (i = 0; i < o; i++) { |
1062 | /* Each LHS vertex is connected to all RHS vertices. */ |
1063 | for (j = 0; j < o; j++) { |
1064 | edges[ne*2] = i; |
1065 | edges[ne*2+1] = j+o; |
1066 | /* capacity for this edge is set later on */ |
1067 | ne++; |
1068 | } |
1069 | } |
1070 | for (i = 0; i < o; i++) { |
1071 | /* Each RHS vertex is connected to the distinguished sink vertex. */ |
1072 | edges[ne*2] = i+o; |
1073 | edges[ne*2+1] = o*2+1; |
1074 | capacity[ne] = 1; |
1075 | ne++; |
1076 | } |
1077 | for (i = 0; i < o; i++) { |
1078 | /* And the distinguished source vertex connects to each LHS vertex. */ |
1079 | edges[ne*2] = o*2; |
1080 | edges[ne*2+1] = i; |
1081 | capacity[ne] = 1; |
1082 | ne++; |
1083 | } |
1084 | assert(ne == o*o + 2*o); |
1085 | /* Now set up backedges. */ |
1086 | maxflow_setup_backedges(ne, edges, backedges); |
1087 | |
1088 | /* |
1089 | * Now generate each row of the latin square. |
1090 | */ |
1091 | for (i = 0; i < o; i++) { |
1092 | /* |
1093 | * To prevent maxflow from behaving deterministically, we |
1094 | * separately permute the columns and the digits for the |
1095 | * purposes of the algorithm, differently for every row. |
1096 | */ |
1097 | for (j = 0; j < o; j++) |
1098 | col[j] = num[j] = j; |
1099 | shuffle(col, j, sizeof(*col), rs); |
3b8fad1a |
1100 | shuffle(num, j, sizeof(*num), rs); |
149255d7 |
1101 | /* We need the num permutation in both forward and inverse forms. */ |
1102 | for (j = 0; j < o; j++) |
1103 | numinv[num[j]] = j; |
1104 | |
1105 | /* |
1106 | * Set up the capacities for the maxflow run, by examining |
1107 | * the existing latin square. |
1108 | */ |
1109 | for (j = 0; j < o*o; j++) |
1110 | capacity[j] = 1; |
1111 | for (j = 0; j < i; j++) |
1112 | for (k = 0; k < o; k++) { |
1113 | int n = num[sq[row[j]*o + col[k]] - 1]; |
1114 | capacity[k*o+n] = 0; |
1115 | } |
1116 | |
1117 | /* |
1118 | * Run maxflow. |
1119 | */ |
1120 | j = maxflow_with_scratch(scratch, o*2+2, 2*o, 2*o+1, ne, |
1121 | edges, backedges, capacity, flow, NULL); |
1122 | assert(j == o); /* by the above theorem, this must have succeeded */ |
1123 | |
1124 | /* |
1125 | * And examine the flow array to pick out the new row of |
1126 | * the latin square. |
1127 | */ |
1128 | for (j = 0; j < o; j++) { |
1129 | for (k = 0; k < o; k++) { |
1130 | if (flow[j*o+k]) |
1131 | break; |
1132 | } |
1133 | assert(k < o); |
1134 | sq[row[i]*o + col[j]] = numinv[k] + 1; |
1135 | } |
1136 | } |
1137 | |
1138 | /* |
1139 | * Done. Free our internal workspaces... |
1140 | */ |
1141 | sfree(flow); |
1142 | sfree(capacity); |
1143 | sfree(edges); |
1144 | sfree(backedges); |
1145 | sfree(scratch); |
1146 | sfree(numinv); |
1147 | sfree(num); |
1148 | sfree(col); |
1149 | sfree(row); |
1150 | |
1151 | /* |
1152 | * ... and return our completed latin square. |
1153 | */ |
1154 | return sq; |
1155 | } |
1156 | |
1157 | /* -------------------------------------------------------- |
1158 | * Checking. |
1159 | */ |
1160 | |
1161 | typedef struct lcparams { |
1162 | digit elt; |
1163 | int count; |
1164 | } lcparams; |
1165 | |
1166 | static int latin_check_cmp(void *v1, void *v2) |
1167 | { |
1168 | lcparams *lc1 = (lcparams *)v1; |
1169 | lcparams *lc2 = (lcparams *)v2; |
1170 | |
1171 | if (lc1->elt < lc2->elt) return -1; |
1172 | if (lc1->elt > lc2->elt) return 1; |
1173 | return 0; |
1174 | } |
1175 | |
1176 | #define ELT(sq,x,y) (sq[((y)*order)+(x)]) |
1177 | |
1178 | /* returns non-zero if sq is not a latin square. */ |
1179 | int latin_check(digit *sq, int order) |
1180 | { |
1181 | tree234 *dict = newtree234(latin_check_cmp); |
1182 | int c, r; |
1183 | int ret = 0; |
242a7d91 |
1184 | lcparams *lcp, lc, *aret; |
149255d7 |
1185 | |
1186 | /* Use a tree234 as a simple hash table, go through the square |
1187 | * adding elements as we go or incrementing their counts. */ |
1188 | for (c = 0; c < order; c++) { |
1189 | for (r = 0; r < order; r++) { |
1190 | lc.elt = ELT(sq, c, r); lc.count = 0; |
1191 | lcp = find234(dict, &lc, NULL); |
1192 | if (!lcp) { |
1193 | lcp = snew(lcparams); |
1194 | lcp->elt = ELT(sq, c, r); |
1195 | lcp->count = 1; |
242a7d91 |
1196 | aret = add234(dict, lcp); |
1197 | assert(aret == lcp); |
149255d7 |
1198 | } else { |
1199 | lcp->count++; |
1200 | } |
1201 | } |
1202 | } |
1203 | |
1204 | /* There should be precisely 'order' letters in the alphabet, |
1205 | * each occurring 'order' times (making the OxO tree) */ |
1206 | if (count234(dict) != order) ret = 1; |
1207 | else { |
1208 | for (c = 0; (lcp = index234(dict, c)) != NULL; c++) { |
1209 | if (lcp->count != order) ret = 1; |
1210 | } |
1211 | } |
1212 | for (c = 0; (lcp = index234(dict, c)) != NULL; c++) |
1213 | sfree(lcp); |
1214 | freetree234(dict); |
1215 | |
1216 | return ret; |
1217 | } |
1218 | |
1219 | |
1220 | /* -------------------------------------------------------- |
1221 | * Testing (and printing). |
1222 | */ |
1223 | |
1224 | #ifdef STANDALONE_LATIN_TEST |
1225 | |
1226 | #include <stdio.h> |
1227 | #include <time.h> |
1228 | |
1229 | const char *quis; |
1230 | |
1231 | static void latin_print(digit *sq, int order) |
1232 | { |
1233 | int x, y; |
1234 | |
1235 | for (y = 0; y < order; y++) { |
1236 | for (x = 0; x < order; x++) { |
1237 | printf("%2u ", ELT(sq, x, y)); |
1238 | } |
1239 | printf("\n"); |
1240 | } |
1241 | printf("\n"); |
1242 | } |
1243 | |
1244 | static void gen(int order, random_state *rs, int debug) |
1245 | { |
1246 | digit *sq; |
1247 | |
1248 | solver_show_working = debug; |
1249 | |
1250 | sq = latin_generate(order, rs); |
1251 | latin_print(sq, order); |
1252 | if (latin_check(sq, order)) { |
1253 | fprintf(stderr, "Square is not a latin square!"); |
1254 | exit(1); |
1255 | } |
1256 | |
1257 | sfree(sq); |
1258 | } |
1259 | |
1260 | void test_soak(int order, random_state *rs) |
1261 | { |
1262 | digit *sq; |
1263 | int n = 0; |
1264 | time_t tt_start, tt_now, tt_last; |
1265 | |
1266 | solver_show_working = 0; |
1267 | tt_now = tt_start = time(NULL); |
1268 | |
1269 | while(1) { |
1270 | sq = latin_generate(order, rs); |
1271 | sfree(sq); |
1272 | n++; |
1273 | |
1274 | tt_last = time(NULL); |
1275 | if (tt_last > tt_now) { |
1276 | tt_now = tt_last; |
1277 | printf("%d total, %3.1f/s\n", n, |
1278 | (double)n / (double)(tt_now - tt_start)); |
1279 | } |
1280 | } |
1281 | } |
1282 | |
1283 | void usage_exit(const char *msg) |
1284 | { |
1285 | if (msg) |
1286 | fprintf(stderr, "%s: %s\n", quis, msg); |
1287 | fprintf(stderr, "Usage: %s [--seed SEED] --soak <params> | [game_id [game_id ...]]\n", quis); |
1288 | exit(1); |
1289 | } |
1290 | |
1291 | int main(int argc, char *argv[]) |
1292 | { |
1293 | int i, soak = 0; |
1294 | random_state *rs; |
1295 | time_t seed = time(NULL); |
1296 | |
1297 | quis = argv[0]; |
1298 | while (--argc > 0) { |
1299 | const char *p = *++argv; |
1300 | if (!strcmp(p, "--soak")) |
1301 | soak = 1; |
1302 | else if (!strcmp(p, "--seed")) { |
1303 | if (argc == 0) |
1304 | usage_exit("--seed needs an argument"); |
1305 | seed = (time_t)atoi(*++argv); |
1306 | argc--; |
1307 | } else if (*p == '-') |
1308 | usage_exit("unrecognised option"); |
1309 | else |
1310 | break; /* finished options */ |
1311 | } |
1312 | |
1313 | rs = random_new((void*)&seed, sizeof(time_t)); |
1314 | |
1315 | if (soak == 1) { |
1316 | if (argc != 1) usage_exit("only one argument for --soak"); |
1317 | test_soak(atoi(*argv), rs); |
1318 | } else { |
1319 | if (argc > 0) { |
1320 | for (i = 0; i < argc; i++) { |
1321 | gen(atoi(*argv++), rs, 1); |
1322 | } |
1323 | } else { |
1324 | while (1) { |
1325 | i = random_upto(rs, 20) + 1; |
1326 | gen(i, rs, 0); |
1327 | } |
1328 | } |
1329 | } |
1330 | random_free(rs); |
1331 | return 0; |
1332 | } |
1333 | |
1334 | #endif |
1335 | |
1336 | /* vim: set shiftwidth=4 tabstop=8: */ |