+ /*
+ * 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 = 0; i < cr; i++) {
+ int p = start2+i*step2;
+ if (usage->cube[p] &&
+ !(p >= start1 && p < start1+cr*step1 && (p - start1) % step1 == 0))
+ {
+#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;
+ py = p / cr;
+ px = py / cr;
+ py %= cr;
+
+ printf("%*s ruling out %d at (%d,%d)\n",
+ solver_recurse_depth*4, "", pn, 1+px, 1+YUNTRANS(py));
+ }
+#endif
+ ret = +1; /* we did something */
+ usage->cube[p] = 0;
+ }
+ }
+
+ return ret;
+}
+
+struct solver_scratch {
+ unsigned char *grid, *rowidx, *colidx, *set;
+ int *mne;
+};
+
+static int solver_set(struct solver_usage *usage,
+ struct solver_scratch *scratch,
+ int start, int step1, int step2
+#ifdef STANDALONE_SOLVER
+ , char *fmt, ...
+#endif
+ )
+{
+ int c = usage->c, r = usage->r, cr = c*r;
+ 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[start+i*step1+j*step2])
+ 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[start+rowidx[i]*step1+colidx[j]*step2];
+
+ /*
+ * 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 = (start+rowidx[i]*step1+
+ colidx[j]*step2);
+#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;
+ py = fpos / cr;
+ px = py / cr;
+ py %= cr;
+
+ printf("%*s ruling out %d at (%d,%d)\n",
+ solver_recurse_depth*4, "",
+ pn, 1+px, 1+YUNTRANS(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;
+}
+
+/*
+ * Try to find a number in the possible set of (x1,y1) which can be
+ * ruled out because it would leave no possibilities for (x2,y2).
+ */
+static int solver_mne(struct solver_usage *usage,
+ struct solver_scratch *scratch,
+ int x1, int y1, int x2, int y2)
+{
+ int c = usage->c, r = usage->r, cr = c*r;
+ int *nb[2];
+ unsigned char *set = scratch->set;
+ unsigned char *numbers = scratch->rowidx;
+ unsigned char *numbersleft = scratch->colidx;
+ int nnb, count;
+ int i, j, n, nbi;
+
+ nb[0] = scratch->mne;
+ nb[1] = scratch->mne + cr;
+
+ /*
+ * First, work out the mutual neighbour squares of the two. We
+ * can assert that they're not actually in the same block,
+ * which leaves two possibilities: they're in different block
+ * rows _and_ different block columns (thus their mutual
+ * neighbours are precisely the other two corners of the
+ * rectangle), or they're in the same row (WLOG) and different
+ * columns, in which case their mutual neighbours are the
+ * column of each block aligned with the other square.
+ *
+ * We divide the mutual neighbours into two separate subsets
+ * nb[0] and nb[1]; squares in the same subset are not only
+ * adjacent to both our key squares, but are also always
+ * adjacent to one another.
+ */
+ if (x1 / r != x2 / r && y1 % r != y2 % r) {
+ /* Corners of the rectangle. */
+ nnb = 1;
+ nb[0][0] = cubepos(x2, y1, 1);
+ nb[1][0] = cubepos(x1, y2, 1);
+ } else if (x1 / r != x2 / r) {
+ /* Same row of blocks; different blocks within that row. */
+ int x1b = x1 - (x1 % r);
+ int x2b = x2 - (x2 % r);
+
+ nnb = r;
+ for (i = 0; i < r; i++) {
+ nb[0][i] = cubepos(x2b+i, y1, 1);
+ nb[1][i] = cubepos(x1b+i, y2, 1);
+ }
+ } else {
+ /* Same column of blocks; different blocks within that column. */
+ int y1b = y1 % r;
+ int y2b = y2 % r;
+
+ assert(y1 % r != y2 % r);
+
+ nnb = c;
+ for (i = 0; i < c; i++) {
+ nb[0][i] = cubepos(x2, y1b+i*r, 1);
+ nb[1][i] = cubepos(x1, y2b+i*r, 1);
+ }
+ }
+
+ /*
+ * Right. Now loop over each possible number.
+ */
+ for (n = 1; n <= cr; n++) {
+ if (!cube(x1, y1, n))
+ continue;
+ for (j = 0; j < cr; j++)
+ numbersleft[j] = cube(x2, y2, j+1);
+
+ /*
+ * Go over every possible subset of each neighbour list,
+ * and see if its union of possible numbers minus n has the
+ * same size as the subset. If so, add the numbers in that
+ * subset to the set of things which would be ruled out
+ * from (x2,y2) if n were placed at (x1,y1).
+ */
+ memset(set, 0, nnb);
+ count = 0;
+ while (1) {
+ /*
+ * Binary increment: change the rightmost 0 to a 1, and
+ * change all 1s to the right of it to 0s.
+ */
+ i = nnb;
+ while (i > 0 && set[i-1])
+ set[--i] = 0, count--;
+ if (i > 0)
+ set[--i] = 1, count++;
+ else
+ break; /* done */
+
+ /*
+ * Examine this subset of each neighbour set.
+ */
+ for (nbi = 0; nbi < 2; nbi++) {
+ int *nbs = nb[nbi];
+
+ memset(numbers, 0, cr);
+
+ for (i = 0; i < nnb; i++)
+ if (set[i])
+ for (j = 0; j < cr; j++)
+ if (j != n-1 && usage->cube[nbs[i] + j])
+ numbers[j] = 1;
+
+ for (i = j = 0; j < cr; j++)
+ i += numbers[j];
+
+ if (i == count) {
+ /*
+ * Got one. This subset of nbs, in the absence
+ * of n, would definitely contain all the
+ * numbers listed in `numbers'. Rule them out
+ * of `numbersleft'.
+ */
+ for (j = 0; j < cr; j++)
+ if (numbers[j])
+ numbersleft[j] = 0;
+ }
+ }
+ }
+
+ /*
+ * If we've got nothing left in `numbersleft', we have a
+ * successful mutual neighbour elimination.
+ */
+ for (j = 0; j < cr; j++)
+ if (numbersleft[j])
+ break;
+
+ if (j == cr) {
+#ifdef STANDALONE_SOLVER
+ if (solver_show_working) {
+ printf("%*smutual neighbour elimination, (%d,%d) vs (%d,%d):\n",
+ solver_recurse_depth*4, "",
+ 1+x1, 1+YUNTRANS(y1), 1+x2, 1+YUNTRANS(y2));
+ printf("%*s ruling out %d at (%d,%d)\n",
+ solver_recurse_depth*4, "",
+ n, 1+x1, 1+YUNTRANS(y1));
+ }
+#endif
+ cube(x1, y1, n) = FALSE;
+ return +1;
+ }
+ }
+
+ return 0; /* nothing found */
+}
+
+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->mne = snewn(2*cr, int);
+ return scratch;
+}
+
+static void solver_free_scratch(struct solver_scratch *scratch)
+{
+ sfree(scratch->mne);
+ sfree(scratch->set);
+ sfree(scratch->colidx);
+ sfree(scratch->rowidx);
+ sfree(scratch->grid);
+ sfree(scratch);
+}
+
+static int solver(int c, int r, digit *grid, int maxdiff)
+{
+ struct solver_usage *usage;
+ struct solver_scratch *scratch;
+ int cr = c*r;
+ int x, y, x2, y2, n, ret;
+ int diff = DIFF_BLOCK;
+
+ /*
+ * Set up a usage structure as a clean slate (everything
+ * possible).
+ */
+ usage = snew(struct solver_usage);
+ usage->c = c;
+ usage->r = r;
+ usage->cr = cr;
+ usage->cube = snewn(cr*cr*cr, unsigned char);
+ usage->grid = grid; /* write straight back to the input */
+ 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);