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 | |
cb0c7d4a |
15 | static void assert_f(p) |
16 | { |
17 | assert(p); |
18 | } |
19 | |
149255d7 |
20 | /* -------------------------------------------------------- |
21 | * Solver. |
22 | */ |
23 | |
24 | /* |
25 | * Function called when we are certain that a particular square has |
26 | * a particular number in it. The y-coordinate passed in here is |
27 | * transformed. |
28 | */ |
29 | void latin_solver_place(struct latin_solver *solver, int x, int y, int n) |
30 | { |
31 | int i, o = solver->o; |
32 | |
33 | assert(n <= o); |
cb0c7d4a |
34 | assert_f(cube(x,y,n)); |
149255d7 |
35 | |
36 | /* |
37 | * Rule out all other numbers in this square. |
38 | */ |
39 | for (i = 1; i <= o; i++) |
40 | if (i != n) |
41 | cube(x,y,i) = FALSE; |
42 | |
43 | /* |
44 | * Rule out this number in all other positions in the row. |
45 | */ |
46 | for (i = 0; i < o; i++) |
47 | if (i != y) |
48 | cube(x,i,n) = FALSE; |
49 | |
50 | /* |
51 | * Rule out this number in all other positions in the column. |
52 | */ |
53 | for (i = 0; i < o; i++) |
54 | if (i != x) |
55 | cube(i,y,n) = FALSE; |
56 | |
57 | /* |
58 | * Enter the number in the result grid. |
59 | */ |
60 | solver->grid[YUNTRANS(y)*o+x] = n; |
61 | |
62 | /* |
63 | * Cross out this number from the list of numbers left to place |
64 | * in its row, its column and its block. |
65 | */ |
66 | solver->row[y*o+n-1] = solver->col[x*o+n-1] = TRUE; |
67 | } |
68 | |
69 | int latin_solver_elim(struct latin_solver *solver, int start, int step |
70 | #ifdef STANDALONE_SOLVER |
71 | , char *fmt, ... |
72 | #endif |
73 | ) |
74 | { |
75 | int o = solver->o; |
76 | int fpos, m, i; |
77 | |
78 | /* |
79 | * Count the number of set bits within this section of the |
80 | * cube. |
81 | */ |
82 | m = 0; |
83 | fpos = -1; |
84 | for (i = 0; i < o; i++) |
85 | if (solver->cube[start+i*step]) { |
86 | fpos = start+i*step; |
87 | m++; |
88 | } |
89 | |
90 | if (m == 1) { |
91 | int x, y, n; |
92 | assert(fpos >= 0); |
93 | |
94 | n = 1 + fpos % o; |
95 | y = fpos / o; |
96 | x = y / o; |
97 | y %= o; |
98 | |
99 | if (!solver->grid[YUNTRANS(y)*o+x]) { |
100 | #ifdef STANDALONE_SOLVER |
101 | if (solver_show_working) { |
102 | va_list ap; |
103 | printf("%*s", solver_recurse_depth*4, ""); |
104 | va_start(ap, fmt); |
105 | vprintf(fmt, ap); |
106 | va_end(ap); |
107 | printf(":\n%*s placing %d at (%d,%d)\n", |
108 | solver_recurse_depth*4, "", n, x, YUNTRANS(y)); |
109 | } |
110 | #endif |
111 | latin_solver_place(solver, x, y, n); |
112 | return +1; |
113 | } |
114 | } else if (m == 0) { |
115 | #ifdef STANDALONE_SOLVER |
116 | if (solver_show_working) { |
117 | va_list ap; |
118 | printf("%*s", solver_recurse_depth*4, ""); |
119 | va_start(ap, fmt); |
120 | vprintf(fmt, ap); |
121 | va_end(ap); |
122 | printf(":\n%*s no possibilities available\n", |
123 | solver_recurse_depth*4, ""); |
124 | } |
125 | #endif |
126 | return -1; |
127 | } |
128 | |
129 | return 0; |
130 | } |
131 | |
132 | struct latin_solver_scratch { |
133 | unsigned char *grid, *rowidx, *colidx, *set; |
134 | int *neighbours, *bfsqueue; |
135 | #ifdef STANDALONE_SOLVER |
136 | int *bfsprev; |
137 | #endif |
138 | }; |
139 | |
140 | int latin_solver_set(struct latin_solver *solver, |
141 | struct latin_solver_scratch *scratch, |
142 | int start, int step1, int step2 |
143 | #ifdef STANDALONE_SOLVER |
144 | , char *fmt, ... |
145 | #endif |
146 | ) |
147 | { |
148 | int o = solver->o; |
149 | int i, j, n, count; |
150 | unsigned char *grid = scratch->grid; |
151 | unsigned char *rowidx = scratch->rowidx; |
152 | unsigned char *colidx = scratch->colidx; |
153 | unsigned char *set = scratch->set; |
154 | |
155 | /* |
156 | * We are passed a o-by-o matrix of booleans. Our first job |
157 | * is to winnow it by finding any definite placements - i.e. |
158 | * any row with a solitary 1 - and discarding that row and the |
159 | * column containing the 1. |
160 | */ |
161 | memset(rowidx, TRUE, o); |
162 | memset(colidx, TRUE, o); |
163 | for (i = 0; i < o; i++) { |
164 | int count = 0, first = -1; |
165 | for (j = 0; j < o; j++) |
166 | if (solver->cube[start+i*step1+j*step2]) |
167 | first = j, count++; |
168 | |
169 | if (count == 0) return -1; |
170 | if (count == 1) |
171 | rowidx[i] = colidx[first] = FALSE; |
172 | } |
173 | |
174 | /* |
175 | * Convert each of rowidx/colidx from a list of 0s and 1s to a |
176 | * list of the indices of the 1s. |
177 | */ |
178 | for (i = j = 0; i < o; i++) |
179 | if (rowidx[i]) |
180 | rowidx[j++] = i; |
181 | n = j; |
182 | for (i = j = 0; i < o; i++) |
183 | if (colidx[i]) |
184 | colidx[j++] = i; |
185 | assert(n == j); |
186 | |
187 | /* |
188 | * And create the smaller matrix. |
189 | */ |
190 | for (i = 0; i < n; i++) |
191 | for (j = 0; j < n; j++) |
192 | grid[i*o+j] = solver->cube[start+rowidx[i]*step1+colidx[j]*step2]; |
193 | |
194 | /* |
195 | * Having done that, we now have a matrix in which every row |
196 | * has at least two 1s in. Now we search to see if we can find |
197 | * a rectangle of zeroes (in the set-theoretic sense of |
198 | * `rectangle', i.e. a subset of rows crossed with a subset of |
199 | * columns) whose width and height add up to n. |
200 | */ |
201 | |
202 | memset(set, 0, n); |
203 | count = 0; |
204 | while (1) { |
205 | /* |
206 | * We have a candidate set. If its size is <=1 or >=n-1 |
207 | * then we move on immediately. |
208 | */ |
209 | if (count > 1 && count < n-1) { |
210 | /* |
211 | * The number of rows we need is n-count. See if we can |
212 | * find that many rows which each have a zero in all |
213 | * the positions listed in `set'. |
214 | */ |
215 | int rows = 0; |
216 | for (i = 0; i < n; i++) { |
217 | int ok = TRUE; |
218 | for (j = 0; j < n; j++) |
219 | if (set[j] && grid[i*o+j]) { |
220 | ok = FALSE; |
221 | break; |
222 | } |
223 | if (ok) |
224 | rows++; |
225 | } |
226 | |
227 | /* |
228 | * We expect never to be able to get _more_ than |
229 | * n-count suitable rows: this would imply that (for |
230 | * example) there are four numbers which between them |
231 | * have at most three possible positions, and hence it |
232 | * indicates a faulty deduction before this point or |
233 | * even a bogus clue. |
234 | */ |
235 | if (rows > n - count) { |
236 | #ifdef STANDALONE_SOLVER |
237 | if (solver_show_working) { |
238 | va_list ap; |
239 | printf("%*s", solver_recurse_depth*4, |
240 | ""); |
241 | va_start(ap, fmt); |
242 | vprintf(fmt, ap); |
243 | va_end(ap); |
244 | printf(":\n%*s contradiction reached\n", |
245 | solver_recurse_depth*4, ""); |
246 | } |
247 | #endif |
248 | return -1; |
249 | } |
250 | |
251 | if (rows >= n - count) { |
252 | int progress = FALSE; |
253 | |
254 | /* |
255 | * We've got one! Now, for each row which _doesn't_ |
256 | * satisfy the criterion, eliminate all its set |
257 | * bits in the positions _not_ listed in `set'. |
258 | * Return +1 (meaning progress has been made) if we |
259 | * successfully eliminated anything at all. |
260 | * |
261 | * This involves referring back through |
262 | * rowidx/colidx in order to work out which actual |
263 | * positions in the cube to meddle with. |
264 | */ |
265 | for (i = 0; i < n; i++) { |
266 | int ok = TRUE; |
267 | for (j = 0; j < n; j++) |
268 | if (set[j] && grid[i*o+j]) { |
269 | ok = FALSE; |
270 | break; |
271 | } |
272 | if (!ok) { |
273 | for (j = 0; j < n; j++) |
274 | if (!set[j] && grid[i*o+j]) { |
275 | int fpos = (start+rowidx[i]*step1+ |
276 | colidx[j]*step2); |
277 | #ifdef STANDALONE_SOLVER |
278 | if (solver_show_working) { |
279 | int px, py, pn; |
280 | |
281 | if (!progress) { |
282 | va_list ap; |
283 | printf("%*s", solver_recurse_depth*4, |
284 | ""); |
285 | va_start(ap, fmt); |
286 | vprintf(fmt, ap); |
287 | va_end(ap); |
288 | printf(":\n"); |
289 | } |
290 | |
291 | pn = 1 + fpos % o; |
292 | py = fpos / o; |
293 | px = py / o; |
294 | py %= o; |
295 | |
296 | printf("%*s ruling out %d at (%d,%d)\n", |
297 | solver_recurse_depth*4, "", |
298 | pn, px, YUNTRANS(py)); |
299 | } |
300 | #endif |
301 | progress = TRUE; |
302 | solver->cube[fpos] = FALSE; |
303 | } |
304 | } |
305 | } |
306 | |
307 | if (progress) { |
308 | return +1; |
309 | } |
310 | } |
311 | } |
312 | |
313 | /* |
314 | * Binary increment: change the rightmost 0 to a 1, and |
315 | * change all 1s to the right of it to 0s. |
316 | */ |
317 | i = n; |
318 | while (i > 0 && set[i-1]) |
319 | set[--i] = 0, count--; |
320 | if (i > 0) |
321 | set[--i] = 1, count++; |
322 | else |
323 | break; /* done */ |
324 | } |
325 | |
326 | return 0; |
327 | } |
328 | |
329 | /* |
330 | * Look for forcing chains. A forcing chain is a path of |
331 | * pairwise-exclusive squares (i.e. each pair of adjacent squares |
332 | * in the path are in the same row, column or block) with the |
333 | * following properties: |
334 | * |
335 | * (a) Each square on the path has precisely two possible numbers. |
336 | * |
337 | * (b) Each pair of squares which are adjacent on the path share |
338 | * at least one possible number in common. |
339 | * |
340 | * (c) Each square in the middle of the path shares _both_ of its |
341 | * numbers with at least one of its neighbours (not the same |
342 | * one with both neighbours). |
343 | * |
344 | * These together imply that at least one of the possible number |
345 | * choices at one end of the path forces _all_ the rest of the |
346 | * numbers along the path. In order to make real use of this, we |
347 | * need further properties: |
348 | * |
349 | * (c) Ruling out some number N from the square at one end |
350 | * of the path forces the square at the other end to |
351 | * take number N. |
352 | * |
353 | * (d) The two end squares are both in line with some third |
354 | * square. |
355 | * |
356 | * (e) That third square currently has N as a possibility. |
357 | * |
358 | * If we can find all of that lot, we can deduce that at least one |
359 | * of the two ends of the forcing chain has number N, and that |
360 | * therefore the mutually adjacent third square does not. |
361 | * |
362 | * To find forcing chains, we're going to start a bfs at each |
363 | * suitable square, once for each of its two possible numbers. |
364 | */ |
365 | int latin_solver_forcing(struct latin_solver *solver, |
366 | struct latin_solver_scratch *scratch) |
367 | { |
368 | int o = solver->o; |
369 | int *bfsqueue = scratch->bfsqueue; |
370 | #ifdef STANDALONE_SOLVER |
371 | int *bfsprev = scratch->bfsprev; |
372 | #endif |
373 | unsigned char *number = scratch->grid; |
374 | int *neighbours = scratch->neighbours; |
375 | int x, y; |
376 | |
377 | for (y = 0; y < o; y++) |
378 | for (x = 0; x < o; x++) { |
379 | int count, t, n; |
380 | |
381 | /* |
382 | * If this square doesn't have exactly two candidate |
383 | * numbers, don't try it. |
384 | * |
385 | * In this loop we also sum the candidate numbers, |
386 | * which is a nasty hack to allow us to quickly find |
387 | * `the other one' (since we will shortly know there |
388 | * are exactly two). |
389 | */ |
390 | for (count = t = 0, n = 1; n <= o; n++) |
391 | if (cube(x, y, n)) |
392 | count++, t += n; |
393 | if (count != 2) |
394 | continue; |
395 | |
396 | /* |
397 | * Now attempt a bfs for each candidate. |
398 | */ |
399 | for (n = 1; n <= o; n++) |
400 | if (cube(x, y, n)) { |
401 | int orign, currn, head, tail; |
402 | |
403 | /* |
404 | * Begin a bfs. |
405 | */ |
406 | orign = n; |
407 | |
408 | memset(number, o+1, o*o); |
409 | head = tail = 0; |
410 | bfsqueue[tail++] = y*o+x; |
411 | #ifdef STANDALONE_SOLVER |
412 | bfsprev[y*o+x] = -1; |
413 | #endif |
414 | number[y*o+x] = t - n; |
415 | |
416 | while (head < tail) { |
417 | int xx, yy, nneighbours, xt, yt, i; |
418 | |
419 | xx = bfsqueue[head++]; |
420 | yy = xx / o; |
421 | xx %= o; |
422 | |
423 | currn = number[yy*o+xx]; |
424 | |
425 | /* |
426 | * Find neighbours of yy,xx. |
427 | */ |
428 | nneighbours = 0; |
429 | for (yt = 0; yt < o; yt++) |
430 | neighbours[nneighbours++] = yt*o+xx; |
431 | for (xt = 0; xt < o; xt++) |
432 | neighbours[nneighbours++] = yy*o+xt; |
433 | |
434 | /* |
435 | * Try visiting each of those neighbours. |
436 | */ |
437 | for (i = 0; i < nneighbours; i++) { |
438 | int cc, tt, nn; |
439 | |
440 | xt = neighbours[i] % o; |
441 | yt = neighbours[i] / o; |
442 | |
443 | /* |
444 | * We need this square to not be |
445 | * already visited, and to include |
446 | * currn as a possible number. |
447 | */ |
448 | if (number[yt*o+xt] <= o) |
449 | continue; |
450 | if (!cube(xt, yt, currn)) |
451 | continue; |
452 | |
453 | /* |
454 | * Don't visit _this_ square a second |
455 | * time! |
456 | */ |
457 | if (xt == xx && yt == yy) |
458 | continue; |
459 | |
460 | /* |
461 | * To continue with the bfs, we need |
462 | * this square to have exactly two |
463 | * possible numbers. |
464 | */ |
465 | for (cc = tt = 0, nn = 1; nn <= o; nn++) |
466 | if (cube(xt, yt, nn)) |
467 | cc++, tt += nn; |
468 | if (cc == 2) { |
469 | bfsqueue[tail++] = yt*o+xt; |
470 | #ifdef STANDALONE_SOLVER |
471 | bfsprev[yt*o+xt] = yy*o+xx; |
472 | #endif |
473 | number[yt*o+xt] = tt - currn; |
474 | } |
475 | |
476 | /* |
477 | * One other possibility is that this |
478 | * might be the square in which we can |
479 | * make a real deduction: if it's |
480 | * adjacent to x,y, and currn is equal |
481 | * to the original number we ruled out. |
482 | */ |
483 | if (currn == orign && |
484 | (xt == x || yt == y)) { |
485 | #ifdef STANDALONE_SOLVER |
486 | if (solver_show_working) { |
487 | char *sep = ""; |
488 | int xl, yl; |
489 | printf("%*sforcing chain, %d at ends of ", |
490 | solver_recurse_depth*4, "", orign); |
491 | xl = xx; |
492 | yl = yy; |
493 | while (1) { |
494 | printf("%s(%d,%d)", sep, xl, |
495 | YUNTRANS(yl)); |
496 | xl = bfsprev[yl*o+xl]; |
497 | if (xl < 0) |
498 | break; |
499 | yl = xl / o; |
500 | xl %= o; |
501 | sep = "-"; |
502 | } |
503 | printf("\n%*s ruling out %d at (%d,%d)\n", |
504 | solver_recurse_depth*4, "", |
505 | orign, xt, YUNTRANS(yt)); |
506 | } |
507 | #endif |
508 | cube(xt, yt, orign) = FALSE; |
509 | return 1; |
510 | } |
511 | } |
512 | } |
513 | } |
514 | } |
515 | |
516 | return 0; |
517 | } |
518 | |
519 | struct latin_solver_scratch *latin_solver_new_scratch(struct latin_solver *solver) |
520 | { |
521 | struct latin_solver_scratch *scratch = snew(struct latin_solver_scratch); |
522 | int o = solver->o; |
523 | scratch->grid = snewn(o*o, unsigned char); |
524 | scratch->rowidx = snewn(o, unsigned char); |
525 | scratch->colidx = snewn(o, unsigned char); |
526 | scratch->set = snewn(o, unsigned char); |
527 | scratch->neighbours = snewn(3*o, int); |
528 | scratch->bfsqueue = snewn(o*o, int); |
529 | #ifdef STANDALONE_SOLVER |
530 | scratch->bfsprev = snewn(o*o, int); |
531 | #endif |
532 | return scratch; |
533 | } |
534 | |
535 | void latin_solver_free_scratch(struct latin_solver_scratch *scratch) |
536 | { |
537 | #ifdef STANDALONE_SOLVER |
538 | sfree(scratch->bfsprev); |
539 | #endif |
540 | sfree(scratch->bfsqueue); |
541 | sfree(scratch->neighbours); |
542 | sfree(scratch->set); |
543 | sfree(scratch->colidx); |
544 | sfree(scratch->rowidx); |
545 | sfree(scratch->grid); |
546 | sfree(scratch); |
547 | } |
548 | |
549 | void latin_solver_alloc(struct latin_solver *solver, digit *grid, int o) |
550 | { |
551 | int x, y; |
552 | |
553 | solver->o = o; |
554 | solver->cube = snewn(o*o*o, unsigned char); |
555 | solver->grid = grid; /* write straight back to the input */ |
556 | memset(solver->cube, TRUE, o*o*o); |
557 | |
558 | solver->row = snewn(o*o, unsigned char); |
559 | solver->col = snewn(o*o, unsigned char); |
560 | memset(solver->row, FALSE, o*o); |
561 | memset(solver->col, FALSE, o*o); |
562 | |
563 | for (x = 0; x < o; x++) |
564 | for (y = 0; y < o; y++) |
565 | if (grid[y*o+x]) |
566 | latin_solver_place(solver, x, YTRANS(y), grid[y*o+x]); |
567 | } |
568 | |
569 | void latin_solver_free(struct latin_solver *solver) |
570 | { |
571 | sfree(solver->cube); |
572 | sfree(solver->row); |
573 | sfree(solver->col); |
574 | } |
575 | |
576 | int latin_solver_diff_simple(struct latin_solver *solver) |
577 | { |
578 | int x, y, n, ret, o = solver->o; |
579 | /* |
580 | * Row-wise positional elimination. |
581 | */ |
582 | for (y = 0; y < o; y++) |
583 | for (n = 1; n <= o; n++) |
584 | if (!solver->row[y*o+n-1]) { |
585 | ret = latin_solver_elim(solver, cubepos(0,y,n), o*o |
586 | #ifdef STANDALONE_SOLVER |
587 | , "positional elimination," |
588 | " %d in row %d", n, YUNTRANS(y) |
589 | #endif |
590 | ); |
591 | if (ret != 0) return ret; |
592 | } |
593 | /* |
594 | * Column-wise positional elimination. |
595 | */ |
596 | for (x = 0; x < o; x++) |
597 | for (n = 1; n <= o; n++) |
598 | if (!solver->col[x*o+n-1]) { |
599 | ret = latin_solver_elim(solver, cubepos(x,0,n), o |
600 | #ifdef STANDALONE_SOLVER |
601 | , "positional elimination," |
602 | " %d in column %d", n, x |
603 | #endif |
604 | ); |
605 | if (ret != 0) return ret; |
606 | } |
607 | |
608 | /* |
609 | * Numeric elimination. |
610 | */ |
611 | for (x = 0; x < o; x++) |
612 | for (y = 0; y < o; y++) |
613 | if (!solver->grid[YUNTRANS(y)*o+x]) { |
614 | ret = latin_solver_elim(solver, cubepos(x,y,1), 1 |
615 | #ifdef STANDALONE_SOLVER |
616 | , "numeric elimination at (%d,%d)", x, |
617 | YUNTRANS(y) |
618 | #endif |
619 | ); |
620 | if (ret != 0) return ret; |
621 | } |
622 | return 0; |
623 | } |
624 | |
625 | int latin_solver_diff_set(struct latin_solver *solver, |
626 | struct latin_solver_scratch *scratch, |
481628b3 |
627 | int extreme) |
149255d7 |
628 | { |
629 | int x, y, n, ret, o = solver->o; |
481628b3 |
630 | |
631 | if (!extreme) { |
632 | /* |
633 | * Row-wise set elimination. |
634 | */ |
635 | for (y = 0; y < o; y++) { |
636 | ret = latin_solver_set(solver, scratch, cubepos(0,y,1), o*o, 1 |
149255d7 |
637 | #ifdef STANDALONE_SOLVER |
481628b3 |
638 | , "set elimination, row %d", YUNTRANS(y) |
149255d7 |
639 | #endif |
481628b3 |
640 | ); |
641 | if (ret != 0) return ret; |
642 | } |
643 | /* |
644 | * Column-wise set elimination. |
645 | */ |
646 | for (x = 0; x < o; x++) { |
647 | ret = latin_solver_set(solver, scratch, cubepos(x,0,1), o, 1 |
149255d7 |
648 | #ifdef STANDALONE_SOLVER |
481628b3 |
649 | , "set elimination, column %d", x |
149255d7 |
650 | #endif |
481628b3 |
651 | ); |
652 | if (ret != 0) return ret; |
653 | } |
654 | } else { |
655 | /* |
656 | * Row-vs-column set elimination on a single number |
657 | * (much tricker for a human to do!) |
658 | */ |
659 | for (n = 1; n <= o; n++) { |
660 | ret = latin_solver_set(solver, scratch, cubepos(0,0,n), o*o, o |
149255d7 |
661 | #ifdef STANDALONE_SOLVER |
481628b3 |
662 | , "positional set elimination, number %d", n |
149255d7 |
663 | #endif |
481628b3 |
664 | ); |
665 | if (ret != 0) return ret; |
666 | } |
149255d7 |
667 | } |
668 | return 0; |
669 | } |
670 | |
671 | /* This uses our own diff_* internally, but doesn't require callers |
672 | * to; this is so it can be used by games that want to rewrite |
673 | * the solver so as to use a different set of difficulties. |
674 | * |
675 | * It returns: |
676 | * 0 for 'didn't do anything' implying it was already solved. |
677 | * -1 for 'impossible' (no solution) |
678 | * 1 for 'single solution' |
679 | * >1 for 'multiple solutions' (you don't get to know how many, and |
680 | * the first such solution found will be set. |
681 | * |
682 | * and this function may well assert if given an impossible board. |
683 | */ |
684 | int latin_solver_recurse(struct latin_solver *solver, int recdiff, |
685 | latin_solver_callback cb, void *ctx) |
686 | { |
687 | int best, bestcount; |
688 | int o = solver->o, x, y, n; |
689 | |
690 | best = -1; |
691 | bestcount = o+1; |
692 | |
693 | for (y = 0; y < o; y++) |
694 | for (x = 0; x < o; x++) |
695 | if (!solver->grid[y*o+x]) { |
696 | int count; |
697 | |
698 | /* |
699 | * An unfilled square. Count the number of |
700 | * possible digits in it. |
701 | */ |
702 | count = 0; |
703 | for (n = 1; n <= o; n++) |
704 | if (cube(x,YTRANS(y),n)) |
705 | count++; |
706 | |
707 | /* |
708 | * We should have found any impossibilities |
709 | * already, so this can safely be an assert. |
710 | */ |
711 | assert(count > 1); |
712 | |
713 | if (count < bestcount) { |
714 | bestcount = count; |
715 | best = y*o+x; |
716 | } |
717 | } |
718 | |
719 | if (best == -1) |
720 | /* we were complete already. */ |
721 | return 0; |
722 | else { |
723 | int i, j; |
724 | digit *list, *ingrid, *outgrid; |
725 | int diff = diff_impossible; /* no solution found yet */ |
726 | |
727 | /* |
728 | * Attempt recursion. |
729 | */ |
730 | y = best / o; |
731 | x = best % o; |
732 | |
733 | list = snewn(o, digit); |
734 | ingrid = snewn(o*o, digit); |
735 | outgrid = snewn(o*o, digit); |
736 | memcpy(ingrid, solver->grid, o*o); |
737 | |
738 | /* Make a list of the possible digits. */ |
739 | for (j = 0, n = 1; n <= o; n++) |
740 | if (cube(x,YTRANS(y),n)) |
741 | list[j++] = n; |
742 | |
743 | #ifdef STANDALONE_SOLVER |
744 | if (solver_show_working) { |
745 | char *sep = ""; |
746 | printf("%*srecursing on (%d,%d) [", |
747 | solver_recurse_depth*4, "", x, y); |
748 | for (i = 0; i < j; i++) { |
749 | printf("%s%d", sep, list[i]); |
750 | sep = " or "; |
751 | } |
752 | printf("]\n"); |
753 | } |
754 | #endif |
755 | |
756 | /* |
757 | * And step along the list, recursing back into the |
758 | * main solver at every stage. |
759 | */ |
760 | for (i = 0; i < j; i++) { |
761 | int ret; |
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", |
769 | solver_recurse_depth*4, "", list[i], x, y); |
770 | solver_recurse_depth++; |
771 | #endif |
772 | |
773 | ret = cb(outgrid, o, recdiff, ctx); |
774 | |
775 | #ifdef STANDALONE_SOLVER |
776 | solver_recurse_depth--; |
777 | if (solver_show_working) { |
778 | printf("%*sretracting %d at (%d,%d)\n", |
779 | solver_recurse_depth*4, "", list[i], x, y); |
780 | } |
781 | #endif |
782 | /* we recurse as deep as we can, so we should never find |
783 | * find ourselves giving up on a puzzle without declaring it |
784 | * impossible. */ |
785 | assert(ret != diff_unfinished); |
786 | |
787 | /* |
788 | * If we have our first solution, copy it into the |
789 | * grid we will return. |
790 | */ |
791 | if (diff == diff_impossible && ret != diff_impossible) |
792 | memcpy(solver->grid, outgrid, o*o); |
793 | |
794 | if (ret == diff_ambiguous) |
795 | diff = diff_ambiguous; |
796 | else if (ret == diff_impossible) |
797 | /* do not change our return value */; |
798 | else { |
799 | /* the recursion turned up exactly one solution */ |
800 | if (diff == diff_impossible) |
801 | diff = recdiff; |
802 | else |
803 | diff = diff_ambiguous; |
804 | } |
805 | |
806 | /* |
807 | * As soon as we've found more than one solution, |
808 | * give up immediately. |
809 | */ |
810 | if (diff == diff_ambiguous) |
811 | break; |
812 | } |
813 | |
814 | sfree(outgrid); |
815 | sfree(ingrid); |
816 | sfree(list); |
817 | |
818 | if (diff == diff_impossible) |
819 | return -1; |
820 | else if (diff == diff_ambiguous) |
821 | return 2; |
822 | else { |
823 | assert(diff == recdiff); |
824 | return 1; |
825 | } |
826 | } |
827 | } |
828 | |
829 | enum { diff_simple = 1, diff_set, diff_extreme, diff_recursive }; |
830 | |
831 | static int latin_solver_sub(struct latin_solver *solver, int maxdiff, void *ctx) |
832 | { |
833 | struct latin_solver_scratch *scratch = latin_solver_new_scratch(solver); |
481628b3 |
834 | int ret, diff = diff_simple; |
149255d7 |
835 | |
836 | assert(maxdiff <= diff_recursive); |
837 | /* |
838 | * Now loop over the grid repeatedly trying all permitted modes |
839 | * of reasoning. The loop terminates if we complete an |
840 | * iteration without making any progress; we then return |
841 | * failure or success depending on whether the grid is full or |
842 | * not. |
843 | */ |
844 | while (1) { |
845 | /* |
846 | * I'd like to write `continue;' inside each of the |
847 | * following loops, so that the solver returns here after |
848 | * making some progress. However, I can't specify that I |
849 | * want to continue an outer loop rather than the innermost |
850 | * one, so I'm apologetically resorting to a goto. |
851 | */ |
852 | cont: |
853 | latin_solver_debug(solver->cube, solver->o); |
854 | |
855 | ret = latin_solver_diff_simple(solver); |
856 | if (ret < 0) { |
857 | diff = diff_impossible; |
858 | goto got_result; |
859 | } else if (ret > 0) { |
860 | diff = max(diff, diff_simple); |
861 | goto cont; |
862 | } |
863 | |
864 | if (maxdiff <= diff_simple) |
865 | break; |
866 | |
481628b3 |
867 | ret = latin_solver_diff_set(solver, scratch, 0); |
149255d7 |
868 | if (ret < 0) { |
869 | diff = diff_impossible; |
870 | goto got_result; |
871 | } else if (ret > 0) { |
481628b3 |
872 | diff = max(diff, diff_set); |
149255d7 |
873 | goto cont; |
874 | } |
875 | |
876 | if (maxdiff <= diff_set) |
877 | break; |
878 | |
481628b3 |
879 | ret = latin_solver_diff_set(solver, scratch, 1); |
880 | if (ret < 0) { |
881 | diff = diff_impossible; |
882 | goto got_result; |
883 | } else if (ret > 0) { |
884 | diff = max(diff, diff_extreme); |
885 | goto cont; |
886 | } |
887 | |
149255d7 |
888 | /* |
889 | * Forcing chains. |
890 | */ |
891 | if (latin_solver_forcing(solver, scratch)) { |
892 | diff = max(diff, diff_extreme); |
893 | goto cont; |
894 | } |
895 | |
896 | /* |
897 | * If we reach here, we have made no deductions in this |
898 | * iteration, so the algorithm terminates. |
899 | */ |
900 | break; |
901 | } |
902 | |
903 | /* |
904 | * Last chance: if we haven't fully solved the puzzle yet, try |
905 | * recursing based on guesses for a particular square. We pick |
906 | * one of the most constrained empty squares we can find, which |
907 | * has the effect of pruning the search tree as much as |
908 | * possible. |
909 | */ |
910 | if (maxdiff == diff_recursive) { |
911 | int nsol = latin_solver_recurse(solver, diff_recursive, latin_solver, ctx); |
912 | if (nsol < 0) diff = diff_impossible; |
913 | else if (nsol == 1) diff = diff_recursive; |
914 | else if (nsol > 1) diff = diff_ambiguous; |
915 | /* if nsol == 0 then we were complete anyway |
916 | * (and thus don't need to change diff) */ |
917 | } else { |
918 | /* |
919 | * We're forbidden to use recursion, so we just see whether |
920 | * our grid is fully solved, and return diff_unfinished |
921 | * otherwise. |
922 | */ |
923 | int x, y, o = solver->o; |
924 | |
925 | for (y = 0; y < o; y++) |
926 | for (x = 0; x < o; x++) |
927 | if (!solver->grid[y*o+x]) |
928 | diff = diff_unfinished; |
929 | } |
930 | |
931 | got_result: |
932 | |
933 | #ifdef STANDALONE_SOLVER |
934 | if (solver_show_working) |
935 | printf("%*s%s found\n", |
936 | solver_recurse_depth*4, "", |
937 | diff == diff_impossible ? "no solution (impossible)" : |
938 | diff == diff_unfinished ? "no solution (unfinished)" : |
939 | diff == diff_ambiguous ? "multiple solutions" : |
940 | "one solution"); |
941 | #endif |
942 | |
943 | latin_solver_free_scratch(scratch); |
944 | |
945 | return diff; |
946 | } |
947 | |
948 | int latin_solver(digit *grid, int o, int maxdiff, void *ctx) |
949 | { |
950 | struct latin_solver solver; |
951 | int diff; |
952 | |
953 | latin_solver_alloc(&solver, grid, o); |
954 | diff = latin_solver_sub(&solver, maxdiff, ctx); |
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 |
962 | if (solver_show_working) { |
963 | struct latin_solver ls, *solver = &ls; |
845a3be0 |
964 | unsigned char *dbg; |
149255d7 |
965 | int x, y, i, c = 0; |
966 | |
967 | ls.cube = cube; ls.o = o; /* for cube() to work */ |
968 | |
969 | dbg = snewn(3*o*o*o, unsigned char); |
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; |
1184 | lcparams *lcp, lc; |
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; |
cb0c7d4a |
1196 | assert_f(add234(dict, lcp) == lcp); |
149255d7 |
1197 | } else { |
1198 | lcp->count++; |
1199 | } |
1200 | } |
1201 | } |
1202 | |
1203 | /* There should be precisely 'order' letters in the alphabet, |
1204 | * each occurring 'order' times (making the OxO tree) */ |
1205 | if (count234(dict) != order) ret = 1; |
1206 | else { |
1207 | for (c = 0; (lcp = index234(dict, c)) != NULL; c++) { |
1208 | if (lcp->count != order) ret = 1; |
1209 | } |
1210 | } |
1211 | for (c = 0; (lcp = index234(dict, c)) != NULL; c++) |
1212 | sfree(lcp); |
1213 | freetree234(dict); |
1214 | |
1215 | return ret; |
1216 | } |
1217 | |
1218 | |
1219 | /* -------------------------------------------------------- |
1220 | * Testing (and printing). |
1221 | */ |
1222 | |
1223 | #ifdef STANDALONE_LATIN_TEST |
1224 | |
1225 | #include <stdio.h> |
1226 | #include <time.h> |
1227 | |
1228 | const char *quis; |
1229 | |
1230 | static void latin_print(digit *sq, int order) |
1231 | { |
1232 | int x, y; |
1233 | |
1234 | for (y = 0; y < order; y++) { |
1235 | for (x = 0; x < order; x++) { |
1236 | printf("%2u ", ELT(sq, x, y)); |
1237 | } |
1238 | printf("\n"); |
1239 | } |
1240 | printf("\n"); |
1241 | } |
1242 | |
1243 | static void gen(int order, random_state *rs, int debug) |
1244 | { |
1245 | digit *sq; |
1246 | |
1247 | solver_show_working = debug; |
1248 | |
1249 | sq = latin_generate(order, rs); |
1250 | latin_print(sq, order); |
1251 | if (latin_check(sq, order)) { |
1252 | fprintf(stderr, "Square is not a latin square!"); |
1253 | exit(1); |
1254 | } |
1255 | |
1256 | sfree(sq); |
1257 | } |
1258 | |
1259 | void test_soak(int order, random_state *rs) |
1260 | { |
1261 | digit *sq; |
1262 | int n = 0; |
1263 | time_t tt_start, tt_now, tt_last; |
1264 | |
1265 | solver_show_working = 0; |
1266 | tt_now = tt_start = time(NULL); |
1267 | |
1268 | while(1) { |
1269 | sq = latin_generate(order, rs); |
1270 | sfree(sq); |
1271 | n++; |
1272 | |
1273 | tt_last = time(NULL); |
1274 | if (tt_last > tt_now) { |
1275 | tt_now = tt_last; |
1276 | printf("%d total, %3.1f/s\n", n, |
1277 | (double)n / (double)(tt_now - tt_start)); |
1278 | } |
1279 | } |
1280 | } |
1281 | |
1282 | void usage_exit(const char *msg) |
1283 | { |
1284 | if (msg) |
1285 | fprintf(stderr, "%s: %s\n", quis, msg); |
1286 | fprintf(stderr, "Usage: %s [--seed SEED] --soak <params> | [game_id [game_id ...]]\n", quis); |
1287 | exit(1); |
1288 | } |
1289 | |
1290 | int main(int argc, char *argv[]) |
1291 | { |
1292 | int i, soak = 0; |
1293 | random_state *rs; |
1294 | time_t seed = time(NULL); |
1295 | |
1296 | quis = argv[0]; |
1297 | while (--argc > 0) { |
1298 | const char *p = *++argv; |
1299 | if (!strcmp(p, "--soak")) |
1300 | soak = 1; |
1301 | else if (!strcmp(p, "--seed")) { |
1302 | if (argc == 0) |
1303 | usage_exit("--seed needs an argument"); |
1304 | seed = (time_t)atoi(*++argv); |
1305 | argc--; |
1306 | } else if (*p == '-') |
1307 | usage_exit("unrecognised option"); |
1308 | else |
1309 | break; /* finished options */ |
1310 | } |
1311 | |
1312 | rs = random_new((void*)&seed, sizeof(time_t)); |
1313 | |
1314 | if (soak == 1) { |
1315 | if (argc != 1) usage_exit("only one argument for --soak"); |
1316 | test_soak(atoi(*argv), rs); |
1317 | } else { |
1318 | if (argc > 0) { |
1319 | for (i = 0; i < argc; i++) { |
1320 | gen(atoi(*argv++), rs, 1); |
1321 | } |
1322 | } else { |
1323 | while (1) { |
1324 | i = random_upto(rs, 20) + 1; |
1325 | gen(i, rs, 0); |
1326 | } |
1327 | } |
1328 | } |
1329 | random_free(rs); |
1330 | return 0; |
1331 | } |
1332 | |
1333 | #endif |
1334 | |
1335 | /* vim: set shiftwidth=4 tabstop=8: */ |