+ if (m == 1) {
+ int x, y, n;
+ assert(fpos >= 0);
+
+ n = 1 + fpos % cr;
+ x = fpos / cr;
+ y = x / cr;
+ x %= cr;
+
+ if (!usage->grid[y*cr+x]) {
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ va_list ap;
+ printf("%*s", solver_recurse_depth*4, "");
+ va_start(ap, fmt);
+ vprintf(fmt, ap);
+ va_end(ap);
+ printf(":\n%*s placing %d at (%d,%d)\n",
+ solver_recurse_depth*4, "", n, 1+x, 1+y);
+ }
+#endif
+ solver_place(usage, x, y, n);
+ return +1;
+ }
+ } else if (m == 0) {
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ va_list ap;
+ printf("%*s", solver_recurse_depth*4, "");
+ va_start(ap, fmt);
+ vprintf(fmt, ap);
+ va_end(ap);
+ printf(":\n%*s no possibilities available\n",
+ solver_recurse_depth*4, "");
+ }
+#endif
+ return -1;
+ }
+
+ return 0;
+}
+
+static int solver_intersect(struct solver_usage *usage,
+ int *indices1, int *indices2
+#ifdef STANDALONE_SOLVER
+ , char *fmt, ...
+#endif
+ )
+{
+ int cr = usage->cr;
+ int ret, i, j;
+
+ /*
+ * Loop over the first domain and see if there's any set bit
+ * not also in the second.
+ */
+ for (i = j = 0; i < cr; i++) {
+ int p = indices1[i];
+ while (j < cr && indices2[j] < p)
+ j++;
+ if (usage->cube[p]) {
+ if (j < cr && indices2[j] == p)
+ continue; /* both domains contain this index */
+ else
+ return 0; /* there is, so we can't deduce */
+ }
+ }
+
+ /*
+ * We have determined that all set bits in the first domain are
+ * within its overlap with the second. So loop over the second
+ * domain and remove all set bits that aren't also in that
+ * overlap; return +1 iff we actually _did_ anything.
+ */
+ ret = 0;
+ for (i = j = 0; i < cr; i++) {
+ int p = indices2[i];
+ while (j < cr && indices1[j] < p)
+ j++;
+ if (usage->cube[p] && (j >= cr || indices1[j] != p)) {
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ int px, py, pn;
+
+ if (!ret) {
+ va_list ap;
+ printf("%*s", solver_recurse_depth*4, "");
+ va_start(ap, fmt);
+ vprintf(fmt, ap);
+ va_end(ap);
+ printf(":\n");
+ }
+
+ pn = 1 + p % cr;
+ px = p / cr;
+ py = px / cr;
+ px %= cr;
+
+ printf("%*s ruling out %d at (%d,%d)\n",
+ solver_recurse_depth*4, "", pn, 1+px, 1+py);
+ }
+#endif
+ ret = +1; /* we did something */
+ usage->cube[p] = 0;
+ }
+ }
+
+ return ret;
+}
+
+struct solver_scratch {
+ unsigned char *grid, *rowidx, *colidx, *set;
+ int *neighbours, *bfsqueue;
+ int *indexlist, *indexlist2;
+#ifdef STANDALONE_SOLVER
+ int *bfsprev;
+#endif
+};
+
+static int solver_set(struct solver_usage *usage,
+ struct solver_scratch *scratch,
+ int *indices
+#ifdef STANDALONE_SOLVER
+ , char *fmt, ...
+#endif
+ )
+{
+ int cr = usage->cr;
+ int i, j, n, count;
+ unsigned char *grid = scratch->grid;
+ unsigned char *rowidx = scratch->rowidx;
+ unsigned char *colidx = scratch->colidx;
+ unsigned char *set = scratch->set;
+
+ /*
+ * We are passed a cr-by-cr matrix of booleans. Our first job
+ * is to winnow it by finding any definite placements - i.e.
+ * any row with a solitary 1 - and discarding that row and the
+ * column containing the 1.
+ */
+ memset(rowidx, TRUE, cr);
+ memset(colidx, TRUE, cr);
+ for (i = 0; i < cr; i++) {
+ int count = 0, first = -1;
+ for (j = 0; j < cr; j++)
+ if (usage->cube[indices[i*cr+j]])
+ first = j, count++;
+
+ /*
+ * If count == 0, then there's a row with no 1s at all and
+ * the puzzle is internally inconsistent. However, we ought
+ * to have caught this already during the simpler reasoning
+ * methods, so we can safely fail an assertion if we reach
+ * this point here.
+ */
+ assert(count > 0);
+ if (count == 1)
+ rowidx[i] = colidx[first] = FALSE;
+ }
+
+ /*
+ * Convert each of rowidx/colidx from a list of 0s and 1s to a
+ * list of the indices of the 1s.
+ */
+ for (i = j = 0; i < cr; i++)
+ if (rowidx[i])
+ rowidx[j++] = i;
+ n = j;
+ for (i = j = 0; i < cr; i++)
+ if (colidx[i])
+ colidx[j++] = i;
+ assert(n == j);
+
+ /*
+ * And create the smaller matrix.
+ */
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ grid[i*cr+j] = usage->cube[indices[rowidx[i]*cr+colidx[j]]];
+
+ /*
+ * Having done that, we now have a matrix in which every row
+ * has at least two 1s in. Now we search to see if we can find
+ * a rectangle of zeroes (in the set-theoretic sense of
+ * `rectangle', i.e. a subset of rows crossed with a subset of
+ * columns) whose width and height add up to n.
+ */
+
+ memset(set, 0, n);
+ count = 0;
+ while (1) {
+ /*
+ * We have a candidate set. If its size is <=1 or >=n-1
+ * then we move on immediately.
+ */
+ if (count > 1 && count < n-1) {
+ /*
+ * The number of rows we need is n-count. See if we can
+ * find that many rows which each have a zero in all
+ * the positions listed in `set'.
+ */
+ int rows = 0;
+ for (i = 0; i < n; i++) {
+ int ok = TRUE;
+ for (j = 0; j < n; j++)
+ if (set[j] && grid[i*cr+j]) {
+ ok = FALSE;
+ break;
+ }
+ if (ok)
+ rows++;
+ }
+
+ /*
+ * We expect never to be able to get _more_ than
+ * n-count suitable rows: this would imply that (for
+ * example) there are four numbers which between them
+ * have at most three possible positions, and hence it
+ * indicates a faulty deduction before this point or
+ * even a bogus clue.
+ */
+ if (rows > n - count) {
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ va_list ap;
+ printf("%*s", solver_recurse_depth*4,
+ "");
+ va_start(ap, fmt);
+ vprintf(fmt, ap);
+ va_end(ap);
+ printf(":\n%*s contradiction reached\n",
+ solver_recurse_depth*4, "");
+ }
+#endif
+ return -1;
+ }
+
+ if (rows >= n - count) {
+ int progress = FALSE;
+
+ /*
+ * We've got one! Now, for each row which _doesn't_
+ * satisfy the criterion, eliminate all its set
+ * bits in the positions _not_ listed in `set'.
+ * Return +1 (meaning progress has been made) if we
+ * successfully eliminated anything at all.
+ *
+ * This involves referring back through
+ * rowidx/colidx in order to work out which actual
+ * positions in the cube to meddle with.
+ */
+ for (i = 0; i < n; i++) {
+ int ok = TRUE;
+ for (j = 0; j < n; j++)
+ if (set[j] && grid[i*cr+j]) {
+ ok = FALSE;
+ break;
+ }
+ if (!ok) {
+ for (j = 0; j < n; j++)
+ if (!set[j] && grid[i*cr+j]) {
+ int fpos = indices[rowidx[i]*cr+colidx[j]];
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ int px, py, pn;
+
+ if (!progress) {
+ va_list ap;
+ printf("%*s", solver_recurse_depth*4,
+ "");
+ va_start(ap, fmt);
+ vprintf(fmt, ap);
+ va_end(ap);
+ printf(":\n");
+ }
+
+ pn = 1 + fpos % cr;
+ px = fpos / cr;
+ py = px / cr;
+ px %= cr;
+
+ printf("%*s ruling out %d at (%d,%d)\n",
+ solver_recurse_depth*4, "",
+ pn, 1+px, 1+py);
+ }
+#endif
+ progress = TRUE;
+ usage->cube[fpos] = FALSE;
+ }
+ }
+ }
+
+ if (progress) {
+ return +1;
+ }
+ }
+ }
+
+ /*
+ * Binary increment: change the rightmost 0 to a 1, and
+ * change all 1s to the right of it to 0s.
+ */
+ i = n;
+ while (i > 0 && set[i-1])
+ set[--i] = 0, count--;
+ if (i > 0)
+ set[--i] = 1, count++;
+ else
+ break; /* done */
+ }
+
+ return 0;
+}
+
+/*
+ * Look for forcing chains. A forcing chain is a path of
+ * pairwise-exclusive squares (i.e. each pair of adjacent squares
+ * in the path are in the same row, column or block) with the
+ * following properties:
+ *
+ * (a) Each square on the path has precisely two possible numbers.
+ *
+ * (b) Each pair of squares which are adjacent on the path share
+ * at least one possible number in common.
+ *
+ * (c) Each square in the middle of the path shares _both_ of its
+ * numbers with at least one of its neighbours (not the same
+ * one with both neighbours).
+ *
+ * These together imply that at least one of the possible number
+ * choices at one end of the path forces _all_ the rest of the
+ * numbers along the path. In order to make real use of this, we
+ * need further properties:
+ *
+ * (c) Ruling out some number N from the square at one end of the
+ * path forces the square at the other end to take the same
+ * number N.
+ *
+ * (d) The two end squares are both in line with some third
+ * square.
+ *
+ * (e) That third square currently has N as a possibility.
+ *
+ * If we can find all of that lot, we can deduce that at least one
+ * of the two ends of the forcing chain has number N, and that
+ * therefore the mutually adjacent third square does not.
+ *
+ * To find forcing chains, we're going to start a bfs at each
+ * suitable square, once for each of its two possible numbers.
+ */
+static int solver_forcing(struct solver_usage *usage,
+ struct solver_scratch *scratch)
+{
+ int cr = usage->cr;
+ int *bfsqueue = scratch->bfsqueue;
+#ifdef STANDALONE_SOLVER
+ int *bfsprev = scratch->bfsprev;
+#endif
+ unsigned char *number = scratch->grid;
+ int *neighbours = scratch->neighbours;
+ int x, y;
+
+ for (y = 0; y < cr; y++)
+ for (x = 0; x < cr; x++) {
+ int count, t, n;
+
+ /*
+ * If this square doesn't have exactly two candidate
+ * numbers, don't try it.
+ *
+ * In this loop we also sum the candidate numbers,
+ * which is a nasty hack to allow us to quickly find
+ * `the other one' (since we will shortly know there
+ * are exactly two).
+ */
+ for (count = t = 0, n = 1; n <= cr; n++)
+ if (cube(x, y, n))
+ count++, t += n;
+ if (count != 2)
+ continue;
+
+ /*
+ * Now attempt a bfs for each candidate.
+ */
+ for (n = 1; n <= cr; n++)
+ if (cube(x, y, n)) {
+ int orign, currn, head, tail;
+
+ /*
+ * Begin a bfs.
+ */
+ orign = n;
+
+ memset(number, cr+1, cr*cr);
+ head = tail = 0;
+ bfsqueue[tail++] = y*cr+x;
+#ifdef STANDALONE_SOLVER
+ bfsprev[y*cr+x] = -1;
+#endif
+ number[y*cr+x] = t - n;
+
+ while (head < tail) {
+ int xx, yy, nneighbours, xt, yt, i;
+
+ xx = bfsqueue[head++];
+ yy = xx / cr;
+ xx %= cr;
+
+ currn = number[yy*cr+xx];
+
+ /*
+ * Find neighbours of yy,xx.
+ */
+ nneighbours = 0;
+ for (yt = 0; yt < cr; yt++)
+ neighbours[nneighbours++] = yt*cr+xx;
+ for (xt = 0; xt < cr; xt++)
+ neighbours[nneighbours++] = yy*cr+xt;
+ xt = usage->blocks->whichblock[yy*cr+xx];
+ for (yt = 0; yt < cr; yt++)
+ neighbours[nneighbours++] = usage->blocks->blocks[xt][yt];
+ if (usage->diag) {
+ int sqindex = yy*cr+xx;
+ if (ondiag0(sqindex)) {
+ for (i = 0; i < cr; i++)
+ neighbours[nneighbours++] = diag0(i);
+ }
+ if (ondiag1(sqindex)) {
+ for (i = 0; i < cr; i++)
+ neighbours[nneighbours++] = diag1(i);
+ }
+ }
+
+ /*
+ * Try visiting each of those neighbours.
+ */
+ for (i = 0; i < nneighbours; i++) {
+ int cc, tt, nn;
+
+ xt = neighbours[i] % cr;
+ yt = neighbours[i] / cr;
+
+ /*
+ * We need this square to not be
+ * already visited, and to include
+ * currn as a possible number.
+ */
+ if (number[yt*cr+xt] <= cr)
+ continue;
+ if (!cube(xt, yt, currn))
+ continue;
+
+ /*
+ * Don't visit _this_ square a second
+ * time!
+ */
+ if (xt == xx && yt == yy)
+ continue;
+
+ /*
+ * To continue with the bfs, we need
+ * this square to have exactly two
+ * possible numbers.
+ */
+ for (cc = tt = 0, nn = 1; nn <= cr; nn++)
+ if (cube(xt, yt, nn))
+ cc++, tt += nn;
+ if (cc == 2) {
+ bfsqueue[tail++] = yt*cr+xt;
+#ifdef STANDALONE_SOLVER
+ bfsprev[yt*cr+xt] = yy*cr+xx;
+#endif
+ number[yt*cr+xt] = tt - currn;
+ }
+
+ /*
+ * One other possibility is that this
+ * might be the square in which we can
+ * make a real deduction: if it's
+ * adjacent to x,y, and currn is equal
+ * to the original number we ruled out.
+ */
+ if (currn == orign &&
+ (xt == x || yt == y ||
+ (usage->blocks->whichblock[yt*cr+xt] == usage->blocks->whichblock[y*cr+x]) ||
+ (usage->diag && ((ondiag0(yt*cr+xt) && ondiag0(y*cr+x)) ||
+ (ondiag1(yt*cr+xt) && ondiag1(y*cr+x)))))) {
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ char *sep = "";
+ int xl, yl;
+ printf("%*sforcing chain, %d at ends of ",
+ solver_recurse_depth*4, "", orign);
+ xl = xx;
+ yl = yy;
+ while (1) {
+ printf("%s(%d,%d)", sep, 1+xl,
+ 1+yl);
+ xl = bfsprev[yl*cr+xl];
+ if (xl < 0)
+ break;
+ yl = xl / cr;
+ xl %= cr;
+ sep = "-";
+ }
+ printf("\n%*s ruling out %d at (%d,%d)\n",
+ solver_recurse_depth*4, "",
+ orign, 1+xt, 1+yt);
+ }
+#endif
+ cube(xt, yt, orign) = FALSE;
+ return 1;
+ }
+ }
+ }
+ }
+ }
+
+ return 0;
+}
+
+static int solver_killer_minmax(struct solver_usage *usage,
+ struct block_structure *cages, digit *clues,
+ int b
+#ifdef STANDALONE_SOLVER
+ , const char *extra
+#endif
+ )
+{
+ int cr = usage->cr;
+ int i;
+ int ret = 0;
+ int nsquares = cages->nr_squares[b];
+
+ if (clues[b] == 0)
+ return 0;
+
+ for (i = 0; i < nsquares; i++) {
+ int n, x = cages->blocks[b][i];
+
+ for (n = 1; n <= cr; n++)
+ if (cube2(x, n)) {
+ int maxval = 0, minval = 0;
+ int j;
+ for (j = 0; j < nsquares; j++) {
+ int m;
+ int y = cages->blocks[b][j];
+ if (i == j)
+ continue;
+ for (m = 1; m <= cr; m++)
+ if (cube2(y, m)) {
+ minval += m;
+ break;
+ }
+ for (m = cr; m > 0; m--)
+ if (cube2(y, m)) {
+ maxval += m;
+ break;
+ }
+ }
+ if (maxval + n < clues[b]) {
+ cube2(x, n) = FALSE;
+ ret = 1;
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working)
+ printf("%*s ruling out %d at (%d,%d) as too low %s\n",
+ solver_recurse_depth*4, "killer minmax analysis",
+ n, 1 + x%cr, 1 + x/cr, extra);
+#endif
+ }
+ if (minval + n > clues[b]) {
+ cube2(x, n) = FALSE;
+ ret = 1;
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working)
+ printf("%*s ruling out %d at (%d,%d) as too high %s\n",
+ solver_recurse_depth*4, "killer minmax analysis",
+ n, 1 + x%cr, 1 + x/cr, extra);
+#endif
+ }
+ }
+ }
+ return ret;
+}
+
+static int solver_killer_sums(struct solver_usage *usage, int b,
+ struct block_structure *cages, int clue,
+ int cage_is_region
+#ifdef STANDALONE_SOLVER
+ , const char *cage_type
+#endif
+ )
+{
+ int cr = usage->cr;
+ int i, ret, max_sums;
+ int nsquares = cages->nr_squares[b];
+ unsigned long *sumbits, possible_addends;
+
+ if (clue == 0) {
+ assert(nsquares == 0);
+ return 0;
+ }
+ assert(nsquares > 0);
+
+ if (nsquares > 4)
+ return 0;
+
+ if (!cage_is_region) {
+ int known_row = -1, known_col = -1, known_block = -1;
+ /*
+ * Verify that the cage lies entirely within one region,
+ * so that using the precomputed sums is valid.
+ */
+ for (i = 0; i < nsquares; i++) {
+ int x = cages->blocks[b][i];
+
+ assert(usage->grid[x] == 0);
+
+ if (i == 0) {
+ known_row = x/cr;
+ known_col = x%cr;
+ known_block = usage->blocks->whichblock[x];
+ } else {
+ if (known_row != x/cr)
+ known_row = -1;
+ if (known_col != x%cr)
+ known_col = -1;
+ if (known_block != usage->blocks->whichblock[x])
+ known_block = -1;
+ }
+ }
+ if (known_block == -1 && known_col == -1 && known_row == -1)
+ return 0;
+ }
+ if (nsquares == 2) {
+ if (clue < 3 || clue > 17)
+ return -1;
+
+ sumbits = sum_bits2[clue];
+ max_sums = MAX_2SUMS;
+ } else if (nsquares == 3) {
+ if (clue < 6 || clue > 24)
+ return -1;
+
+ sumbits = sum_bits3[clue];
+ max_sums = MAX_3SUMS;
+ } else {
+ if (clue < 10 || clue > 30)
+ return -1;
+
+ sumbits = sum_bits4[clue];
+ max_sums = MAX_4SUMS;
+ }
+ /*
+ * For every possible way to get the sum, see if there is
+ * one square in the cage that disallows all the required
+ * addends. If we find one such square, this way to compute
+ * the sum is impossible.
+ */
+ possible_addends = 0;
+ for (i = 0; i < max_sums; i++) {
+ int j;
+ unsigned long bits = sumbits[i];
+
+ if (bits == 0)
+ break;
+
+ for (j = 0; j < nsquares; j++) {
+ int n;
+ unsigned long square_bits = bits;
+ int x = cages->blocks[b][j];
+ for (n = 1; n <= cr; n++)
+ if (!cube2(x, n))
+ square_bits &= ~(1L << n);
+ if (square_bits == 0) {
+ break;
+ }
+ }
+ if (j == nsquares)
+ possible_addends |= bits;
+ }
+ /*
+ * Now we know which addends can possibly be used to
+ * compute the sum. Remove all other digits from the
+ * set of possibilities.
+ */
+ if (possible_addends == 0)
+ return -1;
+
+ ret = 0;
+ for (i = 0; i < nsquares; i++) {
+ int n;
+ int x = cages->blocks[b][i];
+ for (n = 1; n <= cr; n++) {
+ if (!cube2(x, n))
+ continue;
+ if ((possible_addends & (1 << n)) == 0) {
+ cube2(x, n) = FALSE;
+ ret = 1;
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ printf("%*s using %s\n",
+ solver_recurse_depth*4, "killer sums analysis",
+ cage_type);
+ printf("%*s ruling out %d at (%d,%d) due to impossible %d-sum\n",
+ solver_recurse_depth*4, "",
+ n, 1 + x%cr, 1 + x/cr, nsquares);
+ }
+#endif
+ }
+ }
+ }
+ return ret;
+}
+
+static int filter_whole_cages(struct solver_usage *usage, int *squares, int n,
+ int *filtered_sum)
+{
+ int b, i, j, off;
+ *filtered_sum = 0;
+
+ /* First, filter squares with a clue. */
+ for (i = j = 0; i < n; i++)
+ if (usage->grid[squares[i]])
+ *filtered_sum += usage->grid[squares[i]];
+ else
+ squares[j++] = squares[i];
+ n = j;
+
+ /*
+ * Filter all cages that are covered entirely by the list of
+ * squares.
+ */
+ off = 0;
+ for (b = 0; b < usage->kblocks->nr_blocks && off < n; b++) {
+ int b_squares = usage->kblocks->nr_squares[b];
+ int matched = 0;
+
+ if (b_squares == 0)
+ continue;
+
+ /*
+ * Find all squares of block b that lie in our list,
+ * and make them contiguous at off, which is the current position
+ * in the output list.
+ */
+ for (i = 0; i < b_squares; i++) {
+ for (j = off; j < n; j++)
+ if (squares[j] == usage->kblocks->blocks[b][i]) {
+ int t = squares[off + matched];
+ squares[off + matched] = squares[j];
+ squares[j] = t;
+ matched++;
+ break;
+ }
+ }
+ /* If so, filter out all squares of b from the list. */
+ if (matched != usage->kblocks->nr_squares[b]) {
+ off += matched;
+ continue;
+ }
+ memmove(squares + off, squares + off + matched,
+ (n - off - matched) * sizeof *squares);
+ n -= matched;
+
+ *filtered_sum += usage->kclues[b];
+ }
+ assert(off == n);
+ return off;
+}
+
+static struct solver_scratch *solver_new_scratch(struct solver_usage *usage)
+{
+ struct solver_scratch *scratch = snew(struct solver_scratch);
+ int cr = usage->cr;
+ scratch->grid = snewn(cr*cr, unsigned char);
+ scratch->rowidx = snewn(cr, unsigned char);
+ scratch->colidx = snewn(cr, unsigned char);
+ scratch->set = snewn(cr, unsigned char);
+ scratch->neighbours = snewn(5*cr, int);
+ scratch->bfsqueue = snewn(cr*cr, int);
+#ifdef STANDALONE_SOLVER
+ scratch->bfsprev = snewn(cr*cr, int);
+#endif
+ scratch->indexlist = snewn(cr*cr, int); /* used for set elimination */
+ scratch->indexlist2 = snewn(cr, int); /* only used for intersect() */
+ return scratch;
+}
+
+static void solver_free_scratch(struct solver_scratch *scratch)
+{
+#ifdef STANDALONE_SOLVER
+ sfree(scratch->bfsprev);
+#endif
+ sfree(scratch->bfsqueue);
+ sfree(scratch->neighbours);
+ sfree(scratch->set);
+ sfree(scratch->colidx);
+ sfree(scratch->rowidx);
+ sfree(scratch->grid);
+ sfree(scratch->indexlist);
+ sfree(scratch->indexlist2);
+ sfree(scratch);
+}
+
+/*
+ * Used for passing information about difficulty levels between the solver
+ * and its callers.
+ */
+struct difficulty {
+ /* Maximum levels allowed. */
+ int maxdiff, maxkdiff;
+ /* Levels reached by the solver. */
+ int diff, kdiff;
+};
+
+static void solver(int cr, struct block_structure *blocks,
+ struct block_structure *kblocks, int xtype,
+ digit *grid, digit *kgrid, struct difficulty *dlev)
+{
+ struct solver_usage *usage;
+ struct solver_scratch *scratch;
+ int x, y, b, i, n, ret;
+ int diff = DIFF_BLOCK;
+ int kdiff = DIFF_KSINGLE;
+
+ /*
+ * Set up a usage structure as a clean slate (everything
+ * possible).
+ */
+ usage = snew(struct solver_usage);
+ usage->cr = cr;
+ usage->blocks = blocks;
+ if (kblocks) {
+ usage->kblocks = dup_block_structure(kblocks);
+ usage->extra_cages = alloc_block_structure (kblocks->c, kblocks->r,
+ cr * cr, cr, cr * cr);
+ usage->extra_clues = snewn(cr*cr, digit);
+ } else {
+ usage->kblocks = usage->extra_cages = NULL;
+ usage->extra_clues = NULL;
+ }
+ usage->cube = snewn(cr*cr*cr, unsigned char);
+ usage->grid = grid; /* write straight back to the input */
+ if (kgrid) {
+ int nclues;
+
+ assert(kblocks);
+ nclues = kblocks->nr_blocks;
+ /*
+ * Allow for expansion of the killer regions, the absolute
+ * limit is obviously one region per square.
+ */
+ usage->kclues = snewn(cr*cr, digit);
+ for (i = 0; i < nclues; i++) {
+ for (n = 0; n < kblocks->nr_squares[i]; n++)
+ if (kgrid[kblocks->blocks[i][n]] != 0)
+ usage->kclues[i] = kgrid[kblocks->blocks[i][n]];
+ assert(usage->kclues[i] > 0);
+ }
+ memset(usage->kclues + nclues, 0, cr*cr - nclues);
+ } else {
+ usage->kclues = NULL;
+ }
+
+ memset(usage->cube, TRUE, cr*cr*cr);
+
+ usage->row = snewn(cr * cr, unsigned char);
+ usage->col = snewn(cr * cr, unsigned char);
+ usage->blk = snewn(cr * cr, unsigned char);
+ memset(usage->row, FALSE, cr * cr);
+ memset(usage->col, FALSE, cr * cr);
+ memset(usage->blk, FALSE, cr * cr);
+
+ if (xtype) {
+ usage->diag = snewn(cr * 2, unsigned char);
+ memset(usage->diag, FALSE, cr * 2);
+ } else
+ usage->diag = NULL;
+
+ usage->nr_regions = cr * 3 + (xtype ? 2 : 0);
+ usage->regions = snewn(cr * usage->nr_regions, int);
+ usage->sq2region = snewn(cr * cr * 3, int *);
+
+ for (n = 0; n < cr; n++) {
+ for (i = 0; i < cr; i++) {
+ x = n*cr+i;
+ y = i*cr+n;
+ b = usage->blocks->blocks[n][i];
+ usage->regions[cr*n*3 + i] = x;
+ usage->regions[cr*n*3 + cr + i] = y;
+ usage->regions[cr*n*3 + 2*cr + i] = b;
+ usage->sq2region[x*3] = usage->regions + cr*n*3;
+ usage->sq2region[y*3 + 1] = usage->regions + cr*n*3 + cr;
+ usage->sq2region[b*3 + 2] = usage->regions + cr*n*3 + 2*cr;
+ }
+ }
+
+ scratch = solver_new_scratch(usage);
+
+ /*
+ * Place all the clue numbers we are given.
+ */
+ for (x = 0; x < cr; x++)
+ for (y = 0; y < cr; y++)
+ if (grid[y*cr+x])
+ solver_place(usage, x, y, grid[y*cr+x]);
+
+ /*
+ * Now loop over the grid repeatedly trying all permitted modes
+ * of reasoning. The loop terminates if we complete an
+ * iteration without making any progress; we then return
+ * failure or success depending on whether the grid is full or
+ * not.
+ */
+ while (1) {
+ /*
+ * I'd like to write `continue;' inside each of the
+ * following loops, so that the solver returns here after
+ * making some progress. However, I can't specify that I
+ * want to continue an outer loop rather than the innermost
+ * one, so I'm apologetically resorting to a goto.
+ */
+ cont:
+
+ /*
+ * Blockwise positional elimination.
+ */
+ for (b = 0; b < cr; b++)
+ for (n = 1; n <= cr; n++)
+ if (!usage->blk[b*cr+n-1]) {
+ for (i = 0; i < cr; i++)
+ scratch->indexlist[i] = cubepos2(usage->blocks->blocks[b][i],n);
+ ret = solver_elim(usage, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "positional elimination,"
+ " %d in block %s", n,
+ usage->blocks->blocknames[b]
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_BLOCK);
+ goto cont;
+ }
+ }
+
+ if (usage->kclues != NULL) {
+ int changed = FALSE;
+
+ /*
+ * First, bring the kblocks into a more useful form: remove
+ * all filled-in squares, and reduce the sum by their values.
+ * Walk in reverse order, since otherwise remove_from_block
+ * can move element past our loop counter.
+ */
+ for (b = 0; b < usage->kblocks->nr_blocks; b++)
+ for (i = usage->kblocks->nr_squares[b] -1; i >= 0; i--) {
+ int x = usage->kblocks->blocks[b][i];
+ int t = usage->grid[x];
+
+ if (t == 0)
+ continue;
+ remove_from_block(usage->kblocks, b, x);
+ if (t > usage->kclues[b]) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ usage->kclues[b] -= t;
+ /*
+ * Since cages are regions, this tells us something
+ * about the other squares in the cage.
+ */
+ for (n = 0; n < usage->kblocks->nr_squares[b]; n++) {
+ cube2(usage->kblocks->blocks[b][n], t) = FALSE;
+ }
+ }
+
+ /*
+ * The most trivial kind of solver for killer puzzles: fill
+ * single-square cages.
+ */
+ for (b = 0; b < usage->kblocks->nr_blocks; b++) {
+ int squares = usage->kblocks->nr_squares[b];
+ if (squares == 1) {
+ int v = usage->kclues[b];
+ if (v < 1 || v > cr) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ x = usage->kblocks->blocks[b][0] % cr;
+ y = usage->kblocks->blocks[b][0] / cr;
+ if (!cube(x, y, v)) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ solver_place(usage, x, y, v);
+
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ printf("%*s placing %d at (%d,%d)\n",
+ solver_recurse_depth*4, "killer single-square cage",
+ v, 1 + x%cr, 1 + x/cr);
+ }
+#endif
+ changed = TRUE;
+ }
+ }
+
+ if (changed) {
+ kdiff = max(kdiff, DIFF_KSINGLE);
+ goto cont;
+ }
+ }
+ if (dlev->maxkdiff >= DIFF_KINTERSECT && usage->kclues != NULL) {
+ int changed = FALSE;
+ /*
+ * Now, create the extra_cages information. Every full region
+ * (row, column, or block) has the same sum total (45 for 3x3
+ * puzzles. After we try to cover these regions with cages that
+ * lie entirely within them, any squares that remain must bring
+ * the total to this known value, and so they form additional
+ * cages which aren't immediately evident in the displayed form
+ * of the puzzle.
+ */
+ usage->extra_cages->nr_blocks = 0;
+ for (i = 0; i < 3; i++) {
+ for (n = 0; n < cr; n++) {
+ int *region = usage->regions + cr*n*3 + i*cr;
+ int sum = cr * (cr + 1) / 2;
+ int nsquares = cr;
+ int filtered;
+ int n_extra = usage->extra_cages->nr_blocks;
+ int *extra_list = usage->extra_cages->blocks[n_extra];
+ memcpy(extra_list, region, cr * sizeof *extra_list);
+
+ nsquares = filter_whole_cages(usage, extra_list, nsquares, &filtered);
+ sum -= filtered;
+ if (nsquares == cr || nsquares == 0)
+ continue;
+ if (dlev->maxdiff >= DIFF_RECURSIVE) {
+ if (sum <= 0) {
+ dlev->diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ }
+ assert(sum > 0);
+
+ if (nsquares == 1) {
+ if (sum > cr) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ x = extra_list[0] % cr;
+ y = extra_list[0] / cr;
+ if (!cube(x, y, sum)) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ solver_place(usage, x, y, sum);
+ changed = TRUE;
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ printf("%*s placing %d at (%d,%d)\n",
+ solver_recurse_depth*4, "killer single-square deduced cage",
+ sum, 1 + x, 1 + y);
+ }
+#endif
+ }
+
+ b = usage->kblocks->whichblock[extra_list[0]];
+ for (x = 1; x < nsquares; x++)
+ if (usage->kblocks->whichblock[extra_list[x]] != b)
+ break;
+ if (x == nsquares) {
+ assert(usage->kblocks->nr_squares[b] > nsquares);
+ split_block(usage->kblocks, extra_list, nsquares);
+ assert(usage->kblocks->nr_squares[usage->kblocks->nr_blocks - 1] == nsquares);
+ usage->kclues[usage->kblocks->nr_blocks - 1] = sum;
+ usage->kclues[b] -= sum;
+ } else {
+ usage->extra_cages->nr_squares[n_extra] = nsquares;
+ usage->extra_cages->nr_blocks++;
+ usage->extra_clues[n_extra] = sum;
+ }
+ }
+ }
+ if (changed) {
+ kdiff = max(kdiff, DIFF_KINTERSECT);
+ goto cont;
+ }
+ }
+
+ /*
+ * Another simple killer-type elimination. For every square in a
+ * cage, find the minimum and maximum possible sums of all the
+ * other squares in the same cage, and rule out possibilities
+ * for the given square based on whether they are guaranteed to
+ * cause the sum to be either too high or too low.
+ * This is a special case of trying all possible sums across a
+ * region, which is a recursive algorithm. We should probably
+ * implement it for a higher difficulty level.
+ */
+ if (dlev->maxkdiff >= DIFF_KMINMAX && usage->kclues != NULL) {
+ int changed = FALSE;
+ for (b = 0; b < usage->kblocks->nr_blocks; b++) {
+ int ret = solver_killer_minmax(usage, usage->kblocks,
+ usage->kclues, b
+#ifdef STANDALONE_SOLVER
+ , ""
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0)
+ changed = TRUE;
+ }
+ for (b = 0; b < usage->extra_cages->nr_blocks; b++) {
+ int ret = solver_killer_minmax(usage, usage->extra_cages,
+ usage->extra_clues, b
+#ifdef STANDALONE_SOLVER
+ , "using deduced cages"
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0)
+ changed = TRUE;
+ }
+ if (changed) {
+ kdiff = max(kdiff, DIFF_KMINMAX);
+ goto cont;
+ }
+ }
+
+ /*
+ * Try to use knowledge of which numbers can be used to generate
+ * a given sum.
+ * This can only be used if a cage lies entirely within a region.
+ */
+ if (dlev->maxkdiff >= DIFF_KSUMS && usage->kclues != NULL) {
+ int changed = FALSE;
+
+ for (b = 0; b < usage->kblocks->nr_blocks; b++) {
+ int ret = solver_killer_sums(usage, b, usage->kblocks,
+ usage->kclues[b], TRUE
+#ifdef STANDALONE_SOLVER
+ , "regular clues"
+#endif
+ );
+ if (ret > 0) {
+ changed = TRUE;
+ kdiff = max(kdiff, DIFF_KSUMS);
+ } else if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ }
+
+ for (b = 0; b < usage->extra_cages->nr_blocks; b++) {
+ int ret = solver_killer_sums(usage, b, usage->extra_cages,
+ usage->extra_clues[b], FALSE
+#ifdef STANDALONE_SOLVER
+ , "deduced clues"
+#endif
+ );
+ if (ret > 0) {
+ changed = TRUE;
+ kdiff = max(kdiff, DIFF_KINTERSECT);
+ } else if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ }
+ }
+
+ if (changed)
+ goto cont;
+ }
+
+ if (dlev->maxdiff <= DIFF_BLOCK)
+ break;
+
+ /*
+ * Row-wise positional elimination.
+ */
+ for (y = 0; y < cr; y++)
+ for (n = 1; n <= cr; n++)
+ if (!usage->row[y*cr+n-1]) {
+ for (x = 0; x < cr; x++)
+ scratch->indexlist[x] = cubepos(x, y, n);
+ ret = solver_elim(usage, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "positional elimination,"
+ " %d in row %d", n, 1+y
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SIMPLE);
+ goto cont;
+ }
+ }
+ /*
+ * Column-wise positional elimination.
+ */
+ for (x = 0; x < cr; x++)
+ for (n = 1; n <= cr; n++)
+ if (!usage->col[x*cr+n-1]) {
+ for (y = 0; y < cr; y++)
+ scratch->indexlist[y] = cubepos(x, y, n);
+ ret = solver_elim(usage, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "positional elimination,"
+ " %d in column %d", n, 1+x
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SIMPLE);
+ goto cont;
+ }
+ }
+
+ /*
+ * X-diagonal positional elimination.
+ */
+ if (usage->diag) {
+ for (n = 1; n <= cr; n++)
+ if (!usage->diag[n-1]) {
+ for (i = 0; i < cr; i++)
+ scratch->indexlist[i] = cubepos2(diag0(i), n);
+ ret = solver_elim(usage, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "positional elimination,"
+ " %d in \\-diagonal", n
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SIMPLE);
+ goto cont;
+ }
+ }
+ for (n = 1; n <= cr; n++)
+ if (!usage->diag[cr+n-1]) {
+ for (i = 0; i < cr; i++)
+ scratch->indexlist[i] = cubepos2(diag1(i), n);
+ ret = solver_elim(usage, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "positional elimination,"
+ " %d in /-diagonal", n
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SIMPLE);
+ goto cont;
+ }
+ }
+ }
+
+ /*
+ * Numeric elimination.
+ */
+ for (x = 0; x < cr; x++)
+ for (y = 0; y < cr; y++)
+ if (!usage->grid[y*cr+x]) {
+ for (n = 1; n <= cr; n++)
+ scratch->indexlist[n-1] = cubepos(x, y, n);
+ ret = solver_elim(usage, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "numeric elimination at (%d,%d)",
+ 1+x, 1+y
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SIMPLE);
+ goto cont;
+ }
+ }
+
+ if (dlev->maxdiff <= DIFF_SIMPLE)
+ break;
+
+ /*
+ * Intersectional analysis, rows vs blocks.
+ */
+ for (y = 0; y < cr; y++)
+ for (b = 0; b < cr; b++)
+ for (n = 1; n <= cr; n++) {
+ if (usage->row[y*cr+n-1] ||
+ usage->blk[b*cr+n-1])
+ continue;
+ for (i = 0; i < cr; i++) {
+ scratch->indexlist[i] = cubepos(i, y, n);
+ scratch->indexlist2[i] = cubepos2(usage->blocks->blocks[b][i], n);
+ }
+ /*
+ * solver_intersect() never returns -1.
+ */
+ if (solver_intersect(usage, scratch->indexlist,
+ scratch->indexlist2
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in row %d vs block %s",
+ n, 1+y, usage->blocks->blocknames[b]
+#endif
+ ) ||
+ solver_intersect(usage, scratch->indexlist2,
+ scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in block %s vs row %d",
+ n, usage->blocks->blocknames[b], 1+y
+#endif
+ )) {
+ diff = max(diff, DIFF_INTERSECT);
+ goto cont;
+ }
+ }
+
+ /*
+ * Intersectional analysis, columns vs blocks.
+ */
+ for (x = 0; x < cr; x++)
+ for (b = 0; b < cr; b++)
+ for (n = 1; n <= cr; n++) {
+ if (usage->col[x*cr+n-1] ||
+ usage->blk[b*cr+n-1])
+ continue;
+ for (i = 0; i < cr; i++) {
+ scratch->indexlist[i] = cubepos(x, i, n);
+ scratch->indexlist2[i] = cubepos2(usage->blocks->blocks[b][i], n);
+ }
+ if (solver_intersect(usage, scratch->indexlist,
+ scratch->indexlist2
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in column %d vs block %s",
+ n, 1+x, usage->blocks->blocknames[b]
+#endif
+ ) ||
+ solver_intersect(usage, scratch->indexlist2,
+ scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in block %s vs column %d",
+ n, usage->blocks->blocknames[b], 1+x
+#endif
+ )) {
+ diff = max(diff, DIFF_INTERSECT);
+ goto cont;
+ }
+ }
+
+ if (usage->diag) {
+ /*
+ * Intersectional analysis, \-diagonal vs blocks.
+ */
+ for (b = 0; b < cr; b++)
+ for (n = 1; n <= cr; n++) {
+ if (usage->diag[n-1] ||
+ usage->blk[b*cr+n-1])
+ continue;
+ for (i = 0; i < cr; i++) {
+ scratch->indexlist[i] = cubepos2(diag0(i), n);
+ scratch->indexlist2[i] = cubepos2(usage->blocks->blocks[b][i], n);
+ }
+ if (solver_intersect(usage, scratch->indexlist,
+ scratch->indexlist2
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in \\-diagonal vs block %s",
+ n, usage->blocks->blocknames[b]
+#endif
+ ) ||
+ solver_intersect(usage, scratch->indexlist2,
+ scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in block %s vs \\-diagonal",
+ n, usage->blocks->blocknames[b]
+#endif
+ )) {
+ diff = max(diff, DIFF_INTERSECT);
+ goto cont;
+ }
+ }
+
+ /*
+ * Intersectional analysis, /-diagonal vs blocks.
+ */
+ for (b = 0; b < cr; b++)
+ for (n = 1; n <= cr; n++) {
+ if (usage->diag[cr+n-1] ||
+ usage->blk[b*cr+n-1])
+ continue;
+ for (i = 0; i < cr; i++) {
+ scratch->indexlist[i] = cubepos2(diag1(i), n);
+ scratch->indexlist2[i] = cubepos2(usage->blocks->blocks[b][i], n);
+ }
+ if (solver_intersect(usage, scratch->indexlist,
+ scratch->indexlist2
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in /-diagonal vs block %s",
+ n, usage->blocks->blocknames[b]
+#endif
+ ) ||
+ solver_intersect(usage, scratch->indexlist2,
+ scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "intersectional analysis,"
+ " %d in block %s vs /-diagonal",
+ n, usage->blocks->blocknames[b]
+#endif
+ )) {
+ diff = max(diff, DIFF_INTERSECT);
+ goto cont;
+ }
+ }
+ }
+
+ if (dlev->maxdiff <= DIFF_INTERSECT)
+ break;
+
+ /*
+ * Blockwise set elimination.
+ */
+ for (b = 0; b < cr; b++) {
+ for (i = 0; i < cr; i++)
+ for (n = 1; n <= cr; n++)
+ scratch->indexlist[i*cr+n-1] = cubepos2(usage->blocks->blocks[b][i], n);
+ ret = solver_set(usage, scratch, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "set elimination, block %s",
+ usage->blocks->blocknames[b]
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SET);
+ goto cont;
+ }
+ }
+
+ /*
+ * Row-wise set elimination.
+ */
+ for (y = 0; y < cr; y++) {
+ for (x = 0; x < cr; x++)
+ for (n = 1; n <= cr; n++)
+ scratch->indexlist[x*cr+n-1] = cubepos(x, y, n);
+ ret = solver_set(usage, scratch, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "set elimination, row %d", 1+y
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SET);
+ goto cont;
+ }
+ }
+
+ /*
+ * Column-wise set elimination.
+ */
+ for (x = 0; x < cr; x++) {
+ for (y = 0; y < cr; y++)
+ for (n = 1; n <= cr; n++)
+ scratch->indexlist[y*cr+n-1] = cubepos(x, y, n);
+ ret = solver_set(usage, scratch, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "set elimination, column %d", 1+x
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SET);
+ goto cont;
+ }
+ }
+
+ if (usage->diag) {
+ /*
+ * \-diagonal set elimination.
+ */
+ for (i = 0; i < cr; i++)
+ for (n = 1; n <= cr; n++)
+ scratch->indexlist[i*cr+n-1] = cubepos2(diag0(i), n);
+ ret = solver_set(usage, scratch, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "set elimination, \\-diagonal"
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SET);
+ goto cont;
+ }
+
+ /*
+ * /-diagonal set elimination.
+ */
+ for (i = 0; i < cr; i++)
+ for (n = 1; n <= cr; n++)
+ scratch->indexlist[i*cr+n-1] = cubepos2(diag1(i), n);
+ ret = solver_set(usage, scratch, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "set elimination, \\-diagonal"
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_SET);
+ goto cont;
+ }
+ }
+
+ if (dlev->maxdiff <= DIFF_SET)
+ break;
+
+ /*
+ * Row-vs-column set elimination on a single number.
+ */
+ for (n = 1; n <= cr; n++) {
+ for (y = 0; y < cr; y++)
+ for (x = 0; x < cr; x++)
+ scratch->indexlist[y*cr+x] = cubepos(x, y, n);
+ ret = solver_set(usage, scratch, scratch->indexlist
+#ifdef STANDALONE_SOLVER
+ , "positional set elimination, number %d", n
+#endif
+ );
+ if (ret < 0) {
+ diff = DIFF_IMPOSSIBLE;
+ goto got_result;
+ } else if (ret > 0) {
+ diff = max(diff, DIFF_EXTREME);
+ goto cont;
+ }
+ }
+
+ /*
+ * Forcing chains.
+ */
+ if (solver_forcing(usage, scratch)) {
+ diff = max(diff, DIFF_EXTREME);
+ goto cont;
+ }
+
+ /*
+ * If we reach here, we have made no deductions in this
+ * iteration, so the algorithm terminates.
+ */
+ break;
+ }
+
+ /*
+ * Last chance: if we haven't fully solved the puzzle yet, try
+ * recursing based on guesses for a particular square. We pick
+ * one of the most constrained empty squares we can find, which
+ * has the effect of pruning the search tree as much as
+ * possible.
+ */
+ if (dlev->maxdiff >= DIFF_RECURSIVE) {
+ int best, bestcount;
+
+ best = -1;
+ bestcount = cr+1;
+
+ for (y = 0; y < cr; y++)
+ for (x = 0; x < cr; x++)
+ if (!grid[y*cr+x]) {
+ int count;
+
+ /*
+ * An unfilled square. Count the number of
+ * possible digits in it.
+ */
+ count = 0;
+ for (n = 1; n <= cr; n++)
+ if (cube(x,y,n))
+ count++;
+
+ /*
+ * We should have found any impossibilities
+ * already, so this can safely be an assert.
+ */
+ assert(count > 1);
+
+ if (count < bestcount) {
+ bestcount = count;
+ best = y*cr+x;
+ }
+ }
+
+ if (best != -1) {
+ int i, j;
+ digit *list, *ingrid, *outgrid;
+
+ diff = DIFF_IMPOSSIBLE; /* no solution found yet */
+
+ /*
+ * Attempt recursion.
+ */
+ y = best / cr;
+ x = best % cr;
+
+ list = snewn(cr, digit);
+ ingrid = snewn(cr * cr, digit);
+ outgrid = snewn(cr * cr, digit);
+ memcpy(ingrid, grid, cr * cr);
+
+ /* Make a list of the possible digits. */
+ for (j = 0, n = 1; n <= cr; n++)
+ if (cube(x,y,n))
+ list[j++] = n;
+
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ char *sep = "";
+ printf("%*srecursing on (%d,%d) [",
+ solver_recurse_depth*4, "", x + 1, y + 1);
+ for (i = 0; i < j; i++) {
+ printf("%s%d", sep, list[i]);
+ sep = " or ";
+ }
+ printf("]\n");
+ }
+#endif
+
+ /*
+ * And step along the list, recursing back into the
+ * main solver at every stage.
+ */
+ for (i = 0; i < j; i++) {
+ memcpy(outgrid, ingrid, cr * cr);
+ outgrid[y*cr+x] = list[i];
+
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working)
+ printf("%*sguessing %d at (%d,%d)\n",
+ solver_recurse_depth*4, "", list[i], x + 1, y + 1);
+ solver_recurse_depth++;
+#endif
+
+ solver(cr, blocks, kblocks, xtype, outgrid, kgrid, dlev);
+
+#ifdef STANDALONE_SOLVER
+ solver_recurse_depth--;
+ if (solver_show_working) {
+ printf("%*sretracting %d at (%d,%d)\n",
+ solver_recurse_depth*4, "", list[i], x + 1, y + 1);
+ }
+#endif
+
+ /*
+ * If we have our first solution, copy it into the
+ * grid we will return.
+ */
+ if (diff == DIFF_IMPOSSIBLE && dlev->diff != DIFF_IMPOSSIBLE)
+ memcpy(grid, outgrid, cr*cr);
+
+ if (dlev->diff == DIFF_AMBIGUOUS)
+ diff = DIFF_AMBIGUOUS;
+ else if (dlev->diff == DIFF_IMPOSSIBLE)
+ /* do not change our return value */;
+ else {
+ /* the recursion turned up exactly one solution */
+ if (diff == DIFF_IMPOSSIBLE)
+ diff = DIFF_RECURSIVE;
+ else
+ diff = DIFF_AMBIGUOUS;
+ }
+
+ /*
+ * As soon as we've found more than one solution,
+ * give up immediately.
+ */
+ if (diff == DIFF_AMBIGUOUS)
+ break;
+ }
+
+ sfree(outgrid);
+ sfree(ingrid);
+ sfree(list);
+ }
+
+ } else {
+ /*
+ * We're forbidden to use recursion, so we just see whether
+ * our grid is fully solved, and return DIFF_IMPOSSIBLE
+ * otherwise.
+ */
+ for (y = 0; y < cr; y++)
+ for (x = 0; x < cr; x++)
+ if (!grid[y*cr+x])
+ diff = DIFF_IMPOSSIBLE;
+ }
+
+ got_result:
+ dlev->diff = diff;
+ dlev->kdiff = kdiff;
+
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working)
+ printf("%*s%s found\n",
+ solver_recurse_depth*4, "",
+ diff == DIFF_IMPOSSIBLE ? "no solution" :
+ diff == DIFF_AMBIGUOUS ? "multiple solutions" :
+ "one solution");
+#endif
+
+ sfree(usage->sq2region);
+ sfree(usage->regions);
+ sfree(usage->cube);
+ sfree(usage->row);
+ sfree(usage->col);
+ sfree(usage->blk);
+ if (usage->kblocks) {
+ free_block_structure(usage->kblocks);
+ free_block_structure(usage->extra_cages);
+ sfree(usage->extra_clues);
+ }
+ if (usage->kclues) sfree(usage->kclues);
+ sfree(usage);
+
+ solver_free_scratch(scratch);
+}
+
+/* ----------------------------------------------------------------------
+ * End of solver code.
+ */
+
+/* ----------------------------------------------------------------------
+ * Killer set generator.
+ */
+
+/* ----------------------------------------------------------------------
+ * Solo filled-grid generator.
+ *
+ * This grid generator works by essentially trying to solve a grid
+ * starting from no clues, and not worrying that there's more than
+ * one possible solution. Unfortunately, it isn't computationally
+ * feasible to do this by calling the above solver with an empty
+ * grid, because that one needs to allocate a lot of scratch space
+ * at every recursion level. Instead, I have a much simpler
+ * algorithm which I shamelessly copied from a Python solver
+ * written by Andrew Wilkinson (which is GPLed, but I've reused
+ * only ideas and no code). It mostly just does the obvious
+ * recursive thing: pick an empty square, put one of the possible
+ * digits in it, recurse until all squares are filled, backtrack
+ * and change some choices if necessary.
+ *
+ * The clever bit is that every time it chooses which square to
+ * fill in next, it does so by counting the number of _possible_
+ * numbers that can go in each square, and it prioritises so that
+ * it picks a square with the _lowest_ number of possibilities. The
+ * idea is that filling in lots of the obvious bits (particularly
+ * any squares with only one possibility) will cut down on the list
+ * of possibilities for other squares and hence reduce the enormous
+ * search space as much as possible as early as possible.
+ *
+ * The use of bit sets implies that we support puzzles up to a size of
+ * 32x32 (less if anyone finds a 16-bit machine to compile this on).
+ */
+
+/*
+ * Internal data structure used in gridgen to keep track of
+ * progress.
+ */
+struct gridgen_coord { int x, y, r; };
+struct gridgen_usage {
+ int cr;
+ struct block_structure *blocks, *kblocks;
+ /* grid is a copy of the input grid, modified as we go along */
+ digit *grid;
+ /*
+ * Bitsets. In each of them, bit n is set if digit n has been placed
+ * in the corresponding region. row, col and blk are used for all
+ * puzzles. cge is used only for killer puzzles, and diag is used
+ * only for x-type puzzles.
+ * All of these have cr entries, except diag which only has 2,
+ * and cge, which has as many entries as kblocks.
+ */
+ unsigned int *row, *col, *blk, *cge, *diag;
+ /* This lists all the empty spaces remaining in the grid. */
+ struct gridgen_coord *spaces;
+ int nspaces;
+ /* If we need randomisation in the solve, this is our random state. */
+ random_state *rs;
+};
+
+static void gridgen_place(struct gridgen_usage *usage, int x, int y, digit n)
+{
+ unsigned int bit = 1 << n;
+ int cr = usage->cr;
+ usage->row[y] |= bit;
+ usage->col[x] |= bit;
+ usage->blk[usage->blocks->whichblock[y*cr+x]] |= bit;
+ if (usage->cge)
+ usage->cge[usage->kblocks->whichblock[y*cr+x]] |= bit;
+ if (usage->diag) {
+ if (ondiag0(y*cr+x))
+ usage->diag[0] |= bit;
+ if (ondiag1(y*cr+x))
+ usage->diag[1] |= bit;
+ }
+ usage->grid[y*cr+x] = n;
+}
+
+static void gridgen_remove(struct gridgen_usage *usage, int x, int y, digit n)
+{
+ unsigned int mask = ~(1 << n);
+ int cr = usage->cr;
+ usage->row[y] &= mask;
+ usage->col[x] &= mask;
+ usage->blk[usage->blocks->whichblock[y*cr+x]] &= mask;
+ if (usage->cge)
+ usage->cge[usage->kblocks->whichblock[y*cr+x]] &= mask;
+ if (usage->diag) {
+ if (ondiag0(y*cr+x))
+ usage->diag[0] &= mask;
+ if (ondiag1(y*cr+x))
+ usage->diag[1] &= mask;
+ }
+ usage->grid[y*cr+x] = 0;
+}
+
+#define N_SINGLE 32
+
+/*
+ * The real recursive step in the generating function.
+ *
+ * Return values: 1 means solution found, 0 means no solution
+ * found on this branch.
+ */
+static int gridgen_real(struct gridgen_usage *usage, digit *grid, int *steps)
+{
+ int cr = usage->cr;
+ int i, j, n, sx, sy, bestm, bestr, ret;
+ int *digits;
+ unsigned int used;
+
+ /*
+ * Firstly, check for completion! If there are no spaces left
+ * in the grid, we have a solution.
+ */
+ if (usage->nspaces == 0)
+ return TRUE;
+
+ /*
+ * Next, abandon generation if we went over our steps limit.
+ */
+ if (*steps <= 0)
+ return FALSE;
+ (*steps)--;
+
+ /*
+ * Otherwise, there must be at least one space. Find the most
+ * constrained space, using the `r' field as a tie-breaker.
+ */
+ bestm = cr+1; /* so that any space will beat it */
+ bestr = 0;
+ used = ~0;
+ i = sx = sy = -1;
+ for (j = 0; j < usage->nspaces; j++) {
+ int x = usage->spaces[j].x, y = usage->spaces[j].y;
+ unsigned int used_xy;
+ int m;
+
+ m = usage->blocks->whichblock[y*cr+x];
+ used_xy = usage->row[y] | usage->col[x] | usage->blk[m];
+ if (usage->cge != NULL)
+ used_xy |= usage->cge[usage->kblocks->whichblock[y*cr+x]];
+ if (usage->cge != NULL)
+ used_xy |= usage->cge[usage->kblocks->whichblock[y*cr+x]];
+ if (usage->diag != NULL) {
+ if (ondiag0(y*cr+x))
+ used_xy |= usage->diag[0];
+ if (ondiag1(y*cr+x))
+ used_xy |= usage->diag[1];
+ }
+
+ /*
+ * Find the number of digits that could go in this space.
+ */
+ m = 0;
+ for (n = 1; n <= cr; n++) {
+ unsigned int bit = 1 << n;
+ if ((used_xy & bit) == 0)
+ m++;
+ }
+ if (m < bestm || (m == bestm && usage->spaces[j].r < bestr)) {
+ bestm = m;
+ bestr = usage->spaces[j].r;
+ sx = x;
+ sy = y;
+ i = j;
+ used = used_xy;
+ }
+ }
+
+ /*
+ * Swap that square into the final place in the spaces array,
+ * so that decrementing nspaces will remove it from the list.
+ */
+ if (i != usage->nspaces-1) {
+ struct gridgen_coord t;
+ t = usage->spaces[usage->nspaces-1];
+ usage->spaces[usage->nspaces-1] = usage->spaces[i];
+ usage->spaces[i] = t;
+ }
+
+ /*
+ * Now we've decided which square to start our recursion at,
+ * simply go through all possible values, shuffling them
+ * randomly first if necessary.
+ */
+ digits = snewn(bestm, int);
+
+ j = 0;
+ for (n = 1; n <= cr; n++) {
+ unsigned int bit = 1 << n;
+
+ if ((used & bit) == 0)
+ digits[j++] = n;
+ }
+
+ if (usage->rs)
+ shuffle(digits, j, sizeof(*digits), usage->rs);