+ 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 */
+}
+
+/*
+ * 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 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 c = usage->c, r = usage->r, cr = c*r;
+ 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, xblk, 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;
+ xblk = xx - (xx % r);
+ for (yt = yy % r; yt < cr; yt += r)
+ for (xt = xblk; xt < xblk+r; xt++)
+ neighbours[nneighbours++] = yt*cr+xt;
+
+ /*
+ * 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 ||
+ (xt / r == x / r && yt % r == y % r))) {
+#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+YUNTRANS(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+YUNTRANS(yt));
+ }
+#endif
+ cube(xt, yt, orign) = FALSE;
+ return 1;
+ }
+ }
+ }
+ }
+ }
+
+ return 0;
+}
+
+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(3*cr, int);
+ scratch->bfsqueue = snewn(cr*cr, int);
+#ifdef STANDALONE_SOLVER
+ scratch->bfsprev = snewn(cr*cr, int);
+#endif
+ 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);
+}
+
+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);
+ memset(usage->col, FALSE, cr * cr);
+ memset(usage->blk, FALSE, cr * cr);
+
+ scratch = solver_new_scratch(usage);
+
+ /*
+ * Place all the clue numbers we are given.
+ */
+ for (x = 0; x < cr; x++)