9 #ifdef STANDALONE_LATIN_TEST
10 #define STANDALONE_SOLVER
15 static void assert_f(p
)
20 /* --------------------------------------------------------
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
29 void latin_solver_place(struct latin_solver
*solver
, int x
, int y
, int n
)
34 assert_f(cube(x
,y
,n
));
37 * Rule out all other numbers in this square.
39 for (i
= 1; i
<= o
; i
++)
44 * Rule out this number in all other positions in the row.
46 for (i
= 0; i
< o
; i
++)
51 * Rule out this number in all other positions in the column.
53 for (i
= 0; i
< o
; i
++)
58 * Enter the number in the result grid.
60 solver
->grid
[YUNTRANS(y
)*o
+x
] = n
;
63 * Cross out this number from the list of numbers left to place
64 * in its row, its column and its block.
66 solver
->row
[y
*o
+n
-1] = solver
->col
[x
*o
+n
-1] = TRUE
;
69 int latin_solver_elim(struct latin_solver
*solver
, int start
, int step
70 #ifdef STANDALONE_SOLVER
79 * Count the number of set bits within this section of the
84 for (i
= 0; i
< o
; i
++)
85 if (solver
->cube
[start
+i
*step
]) {
99 if (!solver
->grid
[YUNTRANS(y
)*o
+x
]) {
100 #ifdef STANDALONE_SOLVER
101 if (solver_show_working
) {
103 printf("%*s", solver_recurse_depth
*4, "");
107 printf(":\n%*s placing %d at (%d,%d)\n",
108 solver_recurse_depth
*4, "", n
, x
, YUNTRANS(y
));
111 latin_solver_place(solver
, x
, y
, n
);
115 #ifdef STANDALONE_SOLVER
116 if (solver_show_working
) {
118 printf("%*s", solver_recurse_depth
*4, "");
122 printf(":\n%*s no possibilities available\n",
123 solver_recurse_depth
*4, "");
132 struct latin_solver_scratch
{
133 unsigned char *grid
, *rowidx
, *colidx
, *set
;
134 int *neighbours
, *bfsqueue
;
135 #ifdef STANDALONE_SOLVER
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
150 unsigned char *grid
= scratch
->grid
;
151 unsigned char *rowidx
= scratch
->rowidx
;
152 unsigned char *colidx
= scratch
->colidx
;
153 unsigned char *set
= scratch
->set
;
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.
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
])
169 if (count
== 0) return -1;
171 rowidx
[i
] = colidx
[first
] = FALSE
;
175 * Convert each of rowidx/colidx from a list of 0s and 1s to a
176 * list of the indices of the 1s.
178 for (i
= j
= 0; i
< o
; i
++)
182 for (i
= j
= 0; i
< o
; i
++)
188 * And create the smaller matrix.
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
];
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.
206 * We have a candidate set. If its size is <=1 or >=n-1
207 * then we move on immediately.
209 if (count
> 1 && count
< n
-1) {
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'.
216 for (i
= 0; i
< n
; i
++) {
218 for (j
= 0; j
< n
; j
++)
219 if (set
[j
] && grid
[i
*o
+j
]) {
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
235 if (rows
> n
- count
) {
236 #ifdef STANDALONE_SOLVER
237 if (solver_show_working
) {
239 printf("%*s", solver_recurse_depth
*4,
244 printf(":\n%*s contradiction reached\n",
245 solver_recurse_depth
*4, "");
251 if (rows
>= n
- count
) {
252 int progress
= FALSE
;
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.
261 * This involves referring back through
262 * rowidx/colidx in order to work out which actual
263 * positions in the cube to meddle with.
265 for (i
= 0; i
< n
; i
++) {
267 for (j
= 0; j
< n
; j
++)
268 if (set
[j
] && grid
[i
*o
+j
]) {
273 for (j
= 0; j
< n
; j
++)
274 if (!set
[j
] && grid
[i
*o
+j
]) {
275 int fpos
= (start
+rowidx
[i
]*step1
+
277 #ifdef STANDALONE_SOLVER
278 if (solver_show_working
) {
283 printf("%*s", solver_recurse_depth
*4,
296 printf("%*s ruling out %d at (%d,%d)\n",
297 solver_recurse_depth
*4, "",
298 pn
, px
, YUNTRANS(py
));
302 solver
->cube
[fpos
] = FALSE
;
314 * Binary increment: change the rightmost 0 to a 1, and
315 * change all 1s to the right of it to 0s.
318 while (i
> 0 && set
[i
-1])
319 set
[--i
] = 0, count
--;
321 set
[--i
] = 1, count
++;
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:
335 * (a) Each square on the path has precisely two possible numbers.
337 * (b) Each pair of squares which are adjacent on the path share
338 * at least one possible number in common.
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).
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:
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
353 * (d) The two end squares are both in line with some third
356 * (e) That third square currently has N as a possibility.
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.
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.
365 int latin_solver_forcing(struct latin_solver
*solver
,
366 struct latin_solver_scratch
*scratch
)
369 int *bfsqueue
= scratch
->bfsqueue
;
370 #ifdef STANDALONE_SOLVER
371 int *bfsprev
= scratch
->bfsprev
;
373 unsigned char *number
= scratch
->grid
;
374 int *neighbours
= scratch
->neighbours
;
377 for (y
= 0; y
< o
; y
++)
378 for (x
= 0; x
< o
; x
++) {
382 * If this square doesn't have exactly two candidate
383 * numbers, don't try it.
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
390 for (count
= t
= 0, n
= 1; n
<= o
; n
++)
397 * Now attempt a bfs for each candidate.
399 for (n
= 1; n
<= o
; n
++)
401 int orign
, currn
, head
, tail
;
408 memset(number
, o
+1, o
*o
);
410 bfsqueue
[tail
++] = y
*o
+x
;
411 #ifdef STANDALONE_SOLVER
414 number
[y
*o
+x
] = t
- n
;
416 while (head
< tail
) {
417 int xx
, yy
, nneighbours
, xt
, yt
, i
;
419 xx
= bfsqueue
[head
++];
423 currn
= number
[yy
*o
+xx
];
426 * Find neighbours of yy,xx.
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
;
435 * Try visiting each of those neighbours.
437 for (i
= 0; i
< nneighbours
; i
++) {
440 xt
= neighbours
[i
] % o
;
441 yt
= neighbours
[i
] / o
;
444 * We need this square to not be
445 * already visited, and to include
446 * currn as a possible number.
448 if (number
[yt
*o
+xt
] <= o
)
450 if (!cube(xt
, yt
, currn
))
454 * Don't visit _this_ square a second
457 if (xt
== xx
&& yt
== yy
)
461 * To continue with the bfs, we need
462 * this square to have exactly two
465 for (cc
= tt
= 0, nn
= 1; nn
<= o
; nn
++)
466 if (cube(xt
, yt
, nn
))
469 bfsqueue
[tail
++] = yt
*o
+xt
;
470 #ifdef STANDALONE_SOLVER
471 bfsprev
[yt
*o
+xt
] = yy
*o
+xx
;
473 number
[yt
*o
+xt
] = tt
- currn
;
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.
483 if (currn
== orign
&&
484 (xt
== x
|| yt
== y
)) {
485 #ifdef STANDALONE_SOLVER
486 if (solver_show_working
) {
489 printf("%*sforcing chain, %d at ends of ",
490 solver_recurse_depth
*4, "", orign
);
494 printf("%s(%d,%d)", sep
, xl
,
496 xl
= bfsprev
[yl
*o
+xl
];
503 printf("\n%*s ruling out %d at (%d,%d)\n",
504 solver_recurse_depth
*4, "",
505 orign
, xt
, YUNTRANS(yt
));
508 cube(xt
, yt
, orign
) = FALSE
;
519 struct latin_solver_scratch
*latin_solver_new_scratch(struct latin_solver
*solver
)
521 struct latin_solver_scratch
*scratch
= snew(struct latin_solver_scratch
);
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);
535 void latin_solver_free_scratch(struct latin_solver_scratch
*scratch
)
537 #ifdef STANDALONE_SOLVER
538 sfree(scratch
->bfsprev
);
540 sfree(scratch
->bfsqueue
);
541 sfree(scratch
->neighbours
);
543 sfree(scratch
->colidx
);
544 sfree(scratch
->rowidx
);
545 sfree(scratch
->grid
);
549 void latin_solver_alloc(struct latin_solver
*solver
, digit
*grid
, int 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
);
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
);
563 for (x
= 0; x
< o
; x
++)
564 for (y
= 0; y
< o
; y
++)
566 latin_solver_place(solver
, x
, YTRANS(y
), grid
[y
*o
+x
]);
569 void latin_solver_free(struct latin_solver
*solver
)
576 int latin_solver_diff_simple(struct latin_solver
*solver
)
578 int x
, y
, n
, ret
, o
= solver
->o
;
580 * Row-wise positional elimination.
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
)
591 if (ret
!= 0) return ret
;
594 * Column-wise positional elimination.
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
605 if (ret
!= 0) return ret
;
609 * Numeric elimination.
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
,
620 if (ret
!= 0) return ret
;
625 int latin_solver_diff_set(struct latin_solver
*solver
,
626 struct latin_solver_scratch
*scratch
,
629 int x
, y
, n
, ret
, o
= solver
->o
;
633 * Row-wise set elimination.
635 for (y
= 0; y
< o
; y
++) {
636 ret
= latin_solver_set(solver
, scratch
, cubepos(0,y
,1), o
*o
, 1
637 #ifdef STANDALONE_SOLVER
638 , "set elimination, row %d", YUNTRANS(y
)
641 if (ret
!= 0) return ret
;
644 * Column-wise set elimination.
646 for (x
= 0; x
< o
; x
++) {
647 ret
= latin_solver_set(solver
, scratch
, cubepos(x
,0,1), o
, 1
648 #ifdef STANDALONE_SOLVER
649 , "set elimination, column %d", x
652 if (ret
!= 0) return ret
;
656 * Row-vs-column set elimination on a single number
657 * (much tricker for a human to do!)
659 for (n
= 1; n
<= o
; n
++) {
660 ret
= latin_solver_set(solver
, scratch
, cubepos(0,0,n
), o
*o
, o
661 #ifdef STANDALONE_SOLVER
662 , "positional set elimination, number %d", n
665 if (ret
!= 0) return ret
;
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.
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.
682 * and this function may well assert if given an impossible board.
684 int latin_solver_recurse(struct latin_solver
*solver
, int recdiff
,
685 latin_solver_callback cb
, void *ctx
)
688 int o
= solver
->o
, x
, y
, n
;
693 for (y
= 0; y
< o
; y
++)
694 for (x
= 0; x
< o
; x
++)
695 if (!solver
->grid
[y
*o
+x
]) {
699 * An unfilled square. Count the number of
700 * possible digits in it.
703 for (n
= 1; n
<= o
; n
++)
704 if (cube(x
,YTRANS(y
),n
))
708 * We should have found any impossibilities
709 * already, so this can safely be an assert.
713 if (count
< bestcount
) {
720 /* we were complete already. */
724 digit
*list
, *ingrid
, *outgrid
;
725 int diff
= diff_impossible
; /* no solution found yet */
733 list
= snewn(o
, digit
);
734 ingrid
= snewn(o
*o
, digit
);
735 outgrid
= snewn(o
*o
, digit
);
736 memcpy(ingrid
, solver
->grid
, o
*o
);
738 /* Make a list of the possible digits. */
739 for (j
= 0, n
= 1; n
<= o
; n
++)
740 if (cube(x
,YTRANS(y
),n
))
743 #ifdef STANDALONE_SOLVER
744 if (solver_show_working
) {
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
]);
757 * And step along the list, recursing back into the
758 * main solver at every stage.
760 for (i
= 0; i
< j
; i
++) {
763 memcpy(outgrid
, ingrid
, o
*o
);
764 outgrid
[y
*o
+x
] = list
[i
];
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
++;
773 ret
= cb(outgrid
, o
, recdiff
, ctx
);
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
);
782 /* we recurse as deep as we can, so we should never find
783 * find ourselves giving up on a puzzle without declaring it
785 assert(ret
!= diff_unfinished
);
788 * If we have our first solution, copy it into the
789 * grid we will return.
791 if (diff
== diff_impossible
&& ret
!= diff_impossible
)
792 memcpy(solver
->grid
, outgrid
, o
*o
);
794 if (ret
== diff_ambiguous
)
795 diff
= diff_ambiguous
;
796 else if (ret
== diff_impossible
)
797 /* do not change our return value */;
799 /* the recursion turned up exactly one solution */
800 if (diff
== diff_impossible
)
803 diff
= diff_ambiguous
;
807 * As soon as we've found more than one solution,
808 * give up immediately.
810 if (diff
== diff_ambiguous
)
818 if (diff
== diff_impossible
)
820 else if (diff
== diff_ambiguous
)
823 assert(diff
== recdiff
);
829 enum { diff_simple
= 1, diff_set
, diff_extreme
, diff_recursive
};
831 static int latin_solver_sub(struct latin_solver
*solver
, int maxdiff
, void *ctx
)
833 struct latin_solver_scratch
*scratch
= latin_solver_new_scratch(solver
);
834 int ret
, diff
= diff_simple
;
836 assert(maxdiff
<= diff_recursive
);
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
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.
853 latin_solver_debug(solver
->cube
, solver
->o
);
855 ret
= latin_solver_diff_simple(solver
);
857 diff
= diff_impossible
;
859 } else if (ret
> 0) {
860 diff
= max(diff
, diff_simple
);
864 if (maxdiff
<= diff_simple
)
867 ret
= latin_solver_diff_set(solver
, scratch
, 0);
869 diff
= diff_impossible
;
871 } else if (ret
> 0) {
872 diff
= max(diff
, diff_set
);
876 if (maxdiff
<= diff_set
)
879 ret
= latin_solver_diff_set(solver
, scratch
, 1);
881 diff
= diff_impossible
;
883 } else if (ret
> 0) {
884 diff
= max(diff
, diff_extreme
);
891 if (latin_solver_forcing(solver
, scratch
)) {
892 diff
= max(diff
, diff_extreme
);
897 * If we reach here, we have made no deductions in this
898 * iteration, so the algorithm terminates.
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
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) */
919 * We're forbidden to use recursion, so we just see whether
920 * our grid is fully solved, and return diff_unfinished
923 int x
, y
, o
= solver
->o
;
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
;
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" :
943 latin_solver_free_scratch(scratch
);
948 int latin_solver(digit
*grid
, int o
, int maxdiff
, void *ctx
)
950 struct latin_solver solver
;
953 latin_solver_alloc(&solver
, grid
, o
);
954 diff
= latin_solver_sub(&solver
, maxdiff
, ctx
);
955 latin_solver_free(&solver
);
959 void latin_solver_debug(unsigned char *cube
, int o
)
961 #ifdef STANDALONE_SOLVER
962 if (solver_show_working
) {
963 struct latin_solver ls
, *solver
= &ls
;
967 ls
.cube
= cube
; ls
.o
= o
; /* for cube() to work */
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
++) {
991 void latin_debug(digit
*sq
, int o
)
993 #ifdef STANDALONE_SOLVER
994 if (solver_show_working
) {
997 for (y
= 0; y
< o
; y
++) {
998 for (x
= 0; x
< o
; x
++) {
999 printf("%2d ", sq
[y
*o
+x
]);
1008 /* --------------------------------------------------------
1012 digit
*latin_generate(int o
, random_state
*rs
)
1015 int *edges
, *backedges
, *capacity
, *flow
;
1017 int ne
, scratchsize
;
1019 digit
*row
, *col
, *numinv
, *num
;
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.
1031 * To find a viable row at each stage, we can make use of the
1032 * support functions in maxflow.c.
1035 sq
= snewn(o
*o
, digit
);
1038 * In case this method of generation introduces a really subtle
1039 * top-to-bottom directional bias, we'll generate the rows in
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
++)
1048 shuffle(row
, i
, sizeof(*row
), rs
);
1051 * Set up the infrastructure for the maxflow algorithm.
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. */
1061 for (i
= 0; i
< o
; i
++) {
1062 /* Each LHS vertex is connected to all RHS vertices. */
1063 for (j
= 0; j
< o
; j
++) {
1065 edges
[ne
*2+1] = j
+o
;
1066 /* capacity for this edge is set later on */
1070 for (i
= 0; i
< o
; i
++) {
1071 /* Each RHS vertex is connected to the distinguished sink vertex. */
1073 edges
[ne
*2+1] = o
*2+1;
1077 for (i
= 0; i
< o
; i
++) {
1078 /* And the distinguished source vertex connects to each LHS vertex. */
1084 assert(ne
== o
*o
+ 2*o
);
1085 /* Now set up backedges. */
1086 maxflow_setup_backedges(ne
, edges
, backedges
);
1089 * Now generate each row of the latin square.
1091 for (i
= 0; i
< o
; i
++) {
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.
1097 for (j
= 0; j
< o
; j
++)
1098 col
[j
] = num
[j
] = j
;
1099 shuffle(col
, j
, sizeof(*col
), rs
);
1100 shuffle(num
, j
, sizeof(*num
), rs
);
1101 /* We need the num permutation in both forward and inverse forms. */
1102 for (j
= 0; j
< o
; j
++)
1106 * Set up the capacities for the maxflow run, by examining
1107 * the existing latin square.
1109 for (j
= 0; j
< o
*o
; j
++)
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;
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 */
1125 * And examine the flow array to pick out the new row of
1128 for (j
= 0; j
< o
; j
++) {
1129 for (k
= 0; k
< o
; k
++) {
1134 sq
[row
[i
]*o
+ col
[j
]] = numinv
[k
] + 1;
1139 * Done. Free our internal workspaces...
1152 * ... and return our completed latin square.
1157 /* --------------------------------------------------------
1161 typedef struct lcparams
{
1166 static int latin_check_cmp(void *v1
, void *v2
)
1168 lcparams
*lc1
= (lcparams
*)v1
;
1169 lcparams
*lc2
= (lcparams
*)v2
;
1171 if (lc1
->elt
< lc2
->elt
) return -1;
1172 if (lc1
->elt
> lc2
->elt
) return 1;
1176 #define ELT(sq,x,y) (sq[((y)*order)+(x)])
1178 /* returns non-zero if sq is not a latin square. */
1179 int latin_check(digit
*sq
, int order
)
1181 tree234
*dict
= newtree234(latin_check_cmp
);
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
);
1193 lcp
= snew(lcparams
);
1194 lcp
->elt
= ELT(sq
, c
, r
);
1196 assert_f(add234(dict
, lcp
) == lcp
);
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;
1207 for (c
= 0; (lcp
= index234(dict
, c
)) != NULL
; c
++) {
1208 if (lcp
->count
!= order
) ret
= 1;
1211 for (c
= 0; (lcp
= index234(dict
, c
)) != NULL
; c
++)
1219 /* --------------------------------------------------------
1220 * Testing (and printing).
1223 #ifdef STANDALONE_LATIN_TEST
1230 static void latin_print(digit
*sq
, int order
)
1234 for (y
= 0; y
< order
; y
++) {
1235 for (x
= 0; x
< order
; x
++) {
1236 printf("%2u ", ELT(sq
, x
, y
));
1243 static void gen(int order
, random_state
*rs
, int debug
)
1247 solver_show_working
= debug
;
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!");
1259 void test_soak(int order
, random_state
*rs
)
1263 time_t tt_start
, tt_now
, tt_last
;
1265 solver_show_working
= 0;
1266 tt_now
= tt_start
= time(NULL
);
1269 sq
= latin_generate(order
, rs
);
1273 tt_last
= time(NULL
);
1274 if (tt_last
> tt_now
) {
1276 printf("%d total, %3.1f/s\n", n
,
1277 (double)n
/ (double)(tt_now
- tt_start
));
1282 void usage_exit(const char *msg
)
1285 fprintf(stderr
, "%s: %s\n", quis
, msg
);
1286 fprintf(stderr
, "Usage: %s [--seed SEED] --soak <params> | [game_id [game_id ...]]\n", quis
);
1290 int main(int argc
, char *argv
[])
1294 time_t seed
= time(NULL
);
1297 while (--argc
> 0) {
1298 const char *p
= *++argv
;
1299 if (!strcmp(p
, "--soak"))
1301 else if (!strcmp(p
, "--seed")) {
1303 usage_exit("--seed needs an argument");
1304 seed
= (time_t)atoi(*++argv
);
1306 } else if (*p
== '-')
1307 usage_exit("unrecognised option");
1309 break; /* finished options */
1312 rs
= random_new((void*)&seed
, sizeof(time_t));
1315 if (argc
!= 1) usage_exit("only one argument for --soak");
1316 test_soak(atoi(*argv
), rs
);
1319 for (i
= 0; i
< argc
; i
++) {
1320 gen(atoi(*argv
++), rs
, 1);
1324 i
= random_upto(rs
, 20) + 1;
1335 /* vim: set shiftwidth=4 tabstop=8: */