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