11 files changed, 2525 insertions, 496 deletions

diff --git a/apps/plugins/puzzles/SOURCES b/apps/plugins/puzzles/SOURCES
index c9611548f9..5301ceef3d 100644
--- a/apps/plugins/puzzles/SOURCES
+++ b/apps/plugins/puzzles/SOURCES
@@ -19,6 +19,7 @@ src/misc.c
 src/penrose.c
 src/printing.c
 src/random.c
+src/sort.c
 src/tdq.c
 src/tree234.c
 src/version.c
diff --git a/apps/plugins/puzzles/src/benchmark.pl b/apps/plugins/puzzles/src/benchmark.pl
index 98763859e8..7ac48abc25 100755
--- a/apps/plugins/puzzles/src/benchmark.pl
+++ b/apps/plugins/puzzles/src/benchmark.pl
@@ -9,7 +9,7 @@ my @presets = ();
 my %presets = ();
 my $maxval = 0;
 
-while (<>) {
+while (<<>>) {
     chomp;
     if (/^(.*)(#.*): ([\d\.]+)$/) {
         push @presets, $1 unless defined $presets{$1};
diff --git a/apps/plugins/puzzles/src/dominosa.R b/apps/plugins/puzzles/src/dominosa.R
index 99218366e6..b85e7dc1e0 100644
--- a/apps/plugins/puzzles/src/dominosa.R
+++ b/apps/plugins/puzzles/src/dominosa.R
@@ -1,6 +1,6 @@
 # -*- makefile -*-
 
-DOMINOSA_EXTRA = laydomino
+DOMINOSA_EXTRA = laydomino dsf sort findloop
 
 dominosa : [X] GTK COMMON dominosa DOMINOSA_EXTRA dominosa-icon|no-icon
 
@@ -8,6 +8,9 @@ dominosa : [G] WINDOWS COMMON dominosa DOMINOSA_EXTRA dominosa.res|noicon.res
 
 ALL += dominosa[COMBINED] DOMINOSA_EXTRA
 
+dominosasolver :   [U] dominosa[STANDALONE_SOLVER] DOMINOSA_EXTRA STANDALONE
+dominosasolver :   [C] dominosa[STANDALONE_SOLVER] DOMINOSA_EXTRA STANDALONE
+
 !begin am gtk
 GAMES += dominosa
 !end
diff --git a/apps/plugins/puzzles/src/dominosa.c b/apps/plugins/puzzles/src/dominosa.c
index 5f035e9250..67a1d69c91 100644
--- a/apps/plugins/puzzles/src/dominosa.c
+++ b/apps/plugins/puzzles/src/dominosa.c
@@ -5,65 +5,41 @@
  */
 
 /*
- * TODO:
- * 
- *  - improve solver so as to use more interesting forms of
- *    deduction
+ * Further possible deduction types in the solver:
  *
- *     * rule out a domino placement if it would divide an unfilled
- *       region such that at least one resulting region had an odd
- *       area
- *        + Tarjan's bridge-finding algorithm would be a way to find
- *          domino placements that split a connected region in two:
- *          form the graph whose vertices are unpaired squares and
- *          whose edges are potential (not placed but also not ruled
- *          out) dominoes covering two of them, and any bridge in that
- *          graph is a candidate.
- *        + Then, finding any old spanning forest of the unfilled
- *          squares should be sufficient to determine the area parity
- *          of the region that any such placement would cut off.
+ *  * possibly an advanced form of deduce_parity via 2-connectedness.
+ *    We currently deal with areas of the graph with exactly one way
+ *    in and out; but if you have an area with exactly _two_ routes in
+ *    and out of it, then you can at least decide on the _relative_
+ *    parity of the two (either 'these two edges both bisect dominoes
+ *    or neither do', or 'exactly one of these edges bisects a
+ *    domino'). And occasionally that can be enough to let you rule
+ *    out one of the two remaining choices.
+ *     + For example, if both those edges bisect a domino, then those
+ *       two dominoes would also be both the same.
+ *     + Or perhaps between them they rule out all possibilities for
+ *       some other square.
+ *     + Or perhaps they themselves would be duplicates!
+ *     + Or perhaps, on purely geometric grounds, they would box in a
+ *       square to the point where it ended up having to be an
+ *       isolated singleton.
+ *     + The tricky part of this is how you do the graph theory.
+ *       Perhaps a modified form of Tarjan's bridge-finding algorithm
+ *       would work, but I haven't thought through the details.
  *
- *     * set analysis
- *        + look at all unclaimed squares containing a given number
- *        + for each one, find the set of possible numbers that it
- *          can connect to (i.e. each neighbouring tile such that
- *          the placement between it and that neighbour has not yet
- *          been ruled out)
- *        + now proceed similarly to Solo set analysis: try to find
- *          a subset of the squares such that the union of their
- *          possible numbers is the same size as the subset. If so,
- *          rule out those possible numbers for all other squares.
- *           * important wrinkle: the double dominoes complicate
- *             matters. Connecting a number to itself uses up _two_
- *             of the unclaimed squares containing a number. Thus,
- *             when finding the initial subset we must never
- *             include two adjacent squares; and also, when ruling
- *             things out after finding the subset, we must be
- *             careful that we don't rule out precisely the domino
- *             placement that was _included_ in our set!
- *
- *     * playing off the two ends of one potential domino, by
- *       considering the alternatives to that domino that each end
- *       might otherwise be part of.
- *        + if not playing this domino would require each end to be
- *          part of an identical domino, play it. (e.g. the middle of
- *          5-4-4-5)
- *        + if not playing this domino would guarantee that the two
- *          ends between them used up all of some other square's
- *          choices, play it. (e.g. the middle of 2-3-3-1 if another 3
- *          cell can only link to a 2 or a 1)
- *
- *     * identify 'forcing chains', in the sense of any path of cells
- *       each of which has only two possible dominoes to be part of,
- *       and each of those rules out one of the choices for the next
- *       cell. Such a chain has the property that either all the odd
- *       dominoes are placed, or all the even ones are placed; so if
- *       either set of those introduces a conflict (e.g. a dupe within
- *       the chain, or using up all of some other square's choices),
- *       then the whole set can be ruled out, and the other set played
- *       immediately.
- *        + this is of course a generalisation of the previous idea,
- *          which is simply a forcing chain of length 3.
+ *  * possibly an advanced version of set analysis which doesn't have
+ *    to start from squares all having the same number? For example,
+ *    if you have three mutually non-adjacent squares labelled 1,2,3
+ *    such that the numbers adjacent to each are precisely the other
+ *    two, then set analysis can work just fine in principle, and
+ *    tells you that those three squares must overlap the three
+ *    dominoes 1-2, 2-3 and 1-3 in some order, so you can rule out any
+ *    placements of those elsewhere.
+ *     + the difficulty with this is how you avoid it going painfully
+ *       exponential-time. You can't iterate over all the subsets, so
+ *       you'd need some kind of more sophisticated directed search.
+ *     + and the adjacency allowance has to be similarly accounted
+ *       for, which could get tricky to keep track of.
  */
 
 #include <stdio.h>
@@ -84,6 +60,26 @@
 
 #define FLASH_TIME 0.13F
 
+/*
+ * Difficulty levels. I do some macro ickery here to ensure that my
+ * enum and the various forms of my name list always match up.
+ */
+#define DIFFLIST(X)                             \
+    X(TRIVIAL,Trivial,t)                        \
+    X(BASIC,Basic,b)                            \
+    X(HARD,Hard,h)                              \
+    X(EXTREME,Extreme,e)                        \
+    X(AMBIGUOUS,Ambiguous,a)                    \
+    /* end of list */
+#define ENUM(upper,title,lower) DIFF_ ## upper,
+#define TITLE(upper,title,lower) #title,
+#define ENCODE(upper,title,lower) #lower
+#define CONFIG(upper,title,lower) ":" #title
+enum { DIFFLIST(ENUM) DIFFCOUNT };
+static char const *const dominosa_diffnames[] = { DIFFLIST(TITLE) };
+static char const dominosa_diffchars[] = DIFFLIST(ENCODE);
+#define DIFFCONFIG DIFFLIST(CONFIG)
+
 enum {
     COL_BACKGROUND,
     COL_TEXT,
@@ -98,7 +94,7 @@ enum {
 
 struct game_params {
     int n;
-    bool unique;
+    int diff;
 };
 
 struct game_numbers {
@@ -125,35 +121,41 @@ static game_params *default_params(void)
     game_params *ret = snew(game_params);
 
     ret->n = 6;
-    ret->unique = true;
+    ret->diff = DIFF_BASIC;
 
     return ret;
 }
 
-static bool game_fetch_preset(int i, char **name, game_params **params)
+static const struct game_params dominosa_presets[] = {
+    {  3, DIFF_TRIVIAL },
+    {  4, DIFF_TRIVIAL },
+    {  5, DIFF_TRIVIAL },
+    {  6, DIFF_TRIVIAL },
+    {  4, DIFF_BASIC   },
+    {  5, DIFF_BASIC   },
+    {  6, DIFF_BASIC   },
+    {  7, DIFF_BASIC   },
+    {  8, DIFF_BASIC   },
+    {  9, DIFF_BASIC   },
+    {  6, DIFF_HARD    },
+    {  6, DIFF_EXTREME },
+};
+
+static bool game_fetch_preset(int i, char **name, game_params **params_out)
 {
-    game_params *ret;
-    int n;
+    game_params *params;
     char buf[80];
 
-    switch (i) {
-      case 0: n = 3; break;
-      case 1: n = 4; break;
-      case 2: n = 5; break;
-      case 3: n = 6; break;
-      case 4: n = 7; break;
-      case 5: n = 8; break;
-      case 6: n = 9; break;
-      default: return false;
-    }
+    if (i < 0 || i >= lenof(dominosa_presets))
+        return false;
 
-    sprintf(buf, "Up to double-%d", n);
-    *name = dupstr(buf);
+    params = snew(game_params);
+    *params = dominosa_presets[i]; /* structure copy */
 
-    *params = ret = snew(game_params);
-    ret->n = n;
-    ret->unique = true;
+    sprintf(buf, "Order %d, %s", params->n, dominosa_diffnames[params->diff]);
 
+    *name = dupstr(buf);
+    *params_out = params;
     return true;
 }
 
@@ -171,18 +173,36 @@ static game_params *dup_params(const game_params *params)
 
 static void decode_params(game_params *params, char const *string)
 {
-    params->n = atoi(string);
-    while (*string && isdigit((unsigned char)*string)) string++;
-    if (*string == 'a')
-        params->unique = false;
+    const char *p = string;
+
+    params->n = atoi(p);
+    while (*p && isdigit((unsigned char)*p)) p++;
+
+    while (*p) {
+        char c = *p++;
+        if (c == 'a') {
+            /* Legacy encoding from before the difficulty system */
+            params->diff = DIFF_AMBIGUOUS;
+        } else if (c == 'd') {
+            int i;
+            params->diff = DIFFCOUNT+1; /* ...which is invalid */
+            if (*p) {
+                for (i = 0; i < DIFFCOUNT; i++) {
+                    if (*p == dominosa_diffchars[i])
+                        params->diff = i;
+                }
+                p++;
+            }
+        }
+    }
 }
 
 static char *encode_params(const game_params *params, bool full)
 {
     char buf[80];
-    sprintf(buf, "%d", params->n);
-    if (full && !params->unique)
-        strcat(buf, "a");
+    int len = sprintf(buf, "%d", params->n);
+    if (full)
+        len += sprintf(buf + len, "d%c", dominosa_diffchars[params->diff]);
     return dupstr(buf);
 }
 
@@ -198,9 +218,10 @@ static config_item *game_configure(const game_params *params)
     sprintf(buf, "%d", params->n);
     ret[0].u.string.sval = dupstr(buf);
 
-    ret[1].name = "Ensure unique solution";
-    ret[1].type = C_BOOLEAN;
-    ret[1].u.boolean.bval = params->unique;
+    ret[1].name = "Difficulty";
+    ret[1].type = C_CHOICES;
+    ret[1].u.choices.choicenames = DIFFCONFIG;
+    ret[1].u.choices.selected = params->diff;
 
     ret[2].name = NULL;
     ret[2].type = C_END;
@@ -213,7 +234,7 @@ static game_params *custom_params(const config_item *cfg)
     game_params *ret = snew(game_params);
 
     ret->n = atoi(cfg[0].u.string.sval);
-    ret->unique = cfg[1].u.boolean.bval;
+    ret->diff = cfg[1].u.choices.selected;
 
     return ret;
 }
@@ -222,6 +243,8 @@ static const char *validate_params(const game_params *params, bool full)
 {
     if (params->n < 1)
         return "Maximum face number must be at least one";
+    if (params->diff >= DIFFCOUNT)
+        return "Unknown difficulty rating";
     return NULL;
 }
 
@@ -229,361 +252,1993 @@ static const char *validate_params(const game_params *params, bool full)
  * Solver.
  */
 
-static int find_overlaps(int w, int h, int placement, int *set)
+#ifdef STANDALONE_SOLVER
+#define SOLVER_DIAGNOSTICS
+bool solver_diagnostics = false;
+#elif defined SOLVER_DIAGNOSTICS
+const bool solver_diagnostics = true;
+#endif
+
+struct solver_domino;
+struct solver_placement;
+
+/*
+ * Information about a particular domino.
+ */
+struct solver_domino {
+    /* The numbers on the domino, and its index in the dominoes array. */
+    int lo, hi, index;
+
+    /* List of placements not yet ruled out for this domino. */
+    int nplacements;
+    struct solver_placement **placements;
+
+#ifdef SOLVER_DIAGNOSTICS
+    /* A textual name we can easily reuse in solver diagnostics. */
+    char *name;
+#endif
+};
+
+/*
+ * Information about a particular 'placement' (i.e. specific location
+ * that a domino might go in).
+ */
+struct solver_placement {
+    /* The index of this placement in sc->placements. */
+    int index;
+
+    /* The two squares that make up this placement. */
+    struct solver_square *squares[2];
+
+    /* The domino that has to go in this position, if any. */
+    struct solver_domino *domino;
+
+    /* The index of this placement in each square's placements array,
+     * and in that of the domino. */
+    int spi[2], dpi;
+
+    /* Whether this is still considered a possible placement. */
+    bool active;
+
+    /* Other domino placements that overlap with this one. (Maximum 6:
+     * three overlapping each square of the placement.) */
+    int noverlaps;
+    struct solver_placement *overlaps[6];
+
+#ifdef SOLVER_DIAGNOSTICS
+    /* A textual name we can easily reuse in solver diagnostics. */
+    char *name;
+#endif
+};
+
+/*
+ * Information about a particular solver square.
+ */
+struct solver_square {
+    /* The coordinates of the square, and its index in a normal grid array. */
+    int x, y, index;
+
+    /* List of domino placements not yet ruled out for this square. */
+    int nplacements;
+    struct solver_placement *placements[4];
+
+    /* The number in the square. */
+    int number;
+
+#ifdef SOLVER_DIAGNOSTICS
+    /* A textual name we can easily reuse in solver diagnostics. */
+    char *name;
+#endif
+};
+
+struct solver_scratch {
+    int n, dc, pc, w, h, wh;
+    int max_diff_used;
+    struct solver_domino *dominoes;
+    struct solver_placement *placements;
+    struct solver_square *squares;
+    struct solver_placement **domino_placement_lists;
+    struct solver_square **squares_by_number;
+    struct findloopstate *fls;
+    bool squares_by_number_initialised;
+    int *wh_scratch, *pc_scratch, *pc_scratch2, *dc_scratch;
+};
+
+static struct solver_scratch *solver_make_scratch(int n)
 {
-    int x, y, n;
+    int dc = DCOUNT(n), w = n+2, h = n+1, wh = w*h;
+    int pc = (w-1)*h + w*(h-1);
+    struct solver_scratch *sc = snew(struct solver_scratch);
+    int hi, lo, di, x, y, pi, si;
+
+    sc->n = n;
+    sc->dc = dc;
+    sc->pc = pc;
+    sc->w = w;
+    sc->h = h;
+    sc->wh = wh;
+
+    sc->dominoes = snewn(dc, struct solver_domino);
+    sc->placements = snewn(pc, struct solver_placement);
+    sc->squares = snewn(wh, struct solver_square);
+    sc->domino_placement_lists = snewn(pc, struct solver_placement *);
+    sc->fls = findloop_new_state(wh);
+
+    for (di = hi = 0; hi <= n; hi++) {
+        for (lo = 0; lo <= hi; lo++) {
+            assert(di == DINDEX(hi, lo));
+            sc->dominoes[di].hi = hi;
+            sc->dominoes[di].lo = lo;
+            sc->dominoes[di].index = di;
 
-    n = 0;                             /* number of returned placements */
+#ifdef SOLVER_DIAGNOSTICS
+            {
+                char buf[128];
+                sprintf(buf, "%d-%d", hi, lo);
+                sc->dominoes[di].name = dupstr(buf);
+            }
+#endif
 
-    x = placement / 2;
-    y = x / w;
-    x %= w;
+            di++;
+        }
+    }
 
-    if (placement & 1) {
-        /*
-         * Horizontal domino, indexed by its left end.
-         */
-        if (x > 0)
-            set[n++] = placement-2;    /* horizontal domino to the left */
-        if (y > 0)
-            set[n++] = placement-2*w-1;/* vertical domino above left side */
-        if (y+1 < h)
-            set[n++] = placement-1;    /* vertical domino below left side */
-        if (x+2 < w)
-            set[n++] = placement+2;    /* horizontal domino to the right */
-        if (y > 0)
-            set[n++] = placement-2*w+2-1;/* vertical domino above right side */
-        if (y+1 < h)
-            set[n++] = placement+2-1;  /* vertical domino below right side */
-    } else {
-        /*
-         * Vertical domino, indexed by its top end.
-         */
-        if (y > 0)
-            set[n++] = placement-2*w;  /* vertical domino above */
-        if (x > 0)
-            set[n++] = placement-2+1;  /* horizontal domino left of top */
-        if (x+1 < w)
-            set[n++] = placement+1;    /* horizontal domino right of top */
-        if (y+2 < h)
-            set[n++] = placement+2*w;  /* vertical domino below */
-        if (x > 0)
-            set[n++] = placement-2+2*w+1;/* horizontal domino left of bottom */
-        if (x+1 < w)
-            set[n++] = placement+2*w+1;/* horizontal domino right of bottom */
+    for (y = 0; y < h; y++) {
+        for (x = 0; x < w; x++) {
+            struct solver_square *sq = &sc->squares[y*w+x];
+            sq->x = x;
+            sq->y = y;
+            sq->index = y * w + x;
+            sq->nplacements = 0;
+
+#ifdef SOLVER_DIAGNOSTICS
+            {
+                char buf[128];
+                sprintf(buf, "(%d,%d)", x, y);
+                sq->name = dupstr(buf);
+            }
+#endif
+        }
     }
 
-    return n;
+    pi = 0;
+    for (y = 0; y < h-1; y++) {
+        for (x = 0; x < w; x++) {
+            assert(pi < pc);
+            sc->placements[pi].squares[0] = &sc->squares[y*w+x];
+            sc->placements[pi].squares[1] = &sc->squares[(y+1)*w+x];
+#ifdef SOLVER_DIAGNOSTICS
+            {
+                char buf[128];
+                sprintf(buf, "(%d,%d-%d)", x, y, y+1);
+                sc->placements[pi].name = dupstr(buf);
+            }
+#endif
+            pi++;
+        }
+    }
+    for (y = 0; y < h; y++) {
+        for (x = 0; x < w-1; x++) {
+            assert(pi < pc);
+            sc->placements[pi].squares[0] = &sc->squares[y*w+x];
+            sc->placements[pi].squares[1] = &sc->squares[y*w+(x+1)];
+#ifdef SOLVER_DIAGNOSTICS
+            {
+                char buf[128];
+                sprintf(buf, "(%d-%d,%d)", x, x+1, y);
+                sc->placements[pi].name = dupstr(buf);
+            }
+#endif
+            pi++;
+        }
+    }
+    assert(pi == pc);
+
+    /* Set up the full placement lists for all squares, temporarily,
+     * so as to use them to calculate the overlap lists */
+    for (si = 0; si < wh; si++)
+        sc->squares[si].nplacements = 0;
+    for (pi = 0; pi < pc; pi++) {
+        struct solver_placement *p = &sc->placements[pi];
+        for (si = 0; si < 2; si++) {
+            struct solver_square *sq = p->squares[si];
+            p->spi[si] = sq->nplacements;
+            sq->placements[sq->nplacements++] = p;
+        }
+    }
+
+    /* Actually calculate the overlap lists */
+    for (pi = 0; pi < pc; pi++) {
+        struct solver_placement *p = &sc->placements[pi];
+        p->noverlaps = 0;
+        for (si = 0; si < 2; si++) {
+            struct solver_square *sq = p->squares[si];
+            int j;
+            for (j = 0; j < sq->nplacements; j++) {
+                struct solver_placement *q = sq->placements[j];
+                if (q != p)
+                    p->overlaps[p->noverlaps++] = q;
+            }
+        }
+    }
+
+    /* Fill in the index field of the placements */
+    for (pi = 0; pi < pc; pi++)
+        sc->placements[pi].index = pi;
+
+    /* Lazily initialised by particular solver techniques that might
+     * never be needed */
+    sc->squares_by_number = NULL;
+    sc->squares_by_number_initialised = false;
+    sc->wh_scratch = NULL;
+    sc->pc_scratch = sc->pc_scratch2 = NULL;
+    sc->dc_scratch = NULL;
+
+    return sc;
+}
+
+static void solver_free_scratch(struct solver_scratch *sc)
+{
+#ifdef SOLVER_DIAGNOSTICS
+    {
+        int i;
+        for (i = 0; i < sc->dc; i++)
+            sfree(sc->dominoes[i].name);
+        for (i = 0; i < sc->pc; i++)
+            sfree(sc->placements[i].name);
+        for (i = 0; i < sc->wh; i++)
+            sfree(sc->squares[i].name);
+    }
+#endif
+    sfree(sc->dominoes);
+    sfree(sc->placements);
+    sfree(sc->squares);
+    sfree(sc->domino_placement_lists);
+    sfree(sc->squares_by_number);
+    findloop_free_state(sc->fls);
+    sfree(sc->wh_scratch);
+    sfree(sc->pc_scratch);
+    sfree(sc->pc_scratch2);
+    sfree(sc->dc_scratch);
+    sfree(sc);
+}
+
+static void solver_setup_grid(struct solver_scratch *sc, const int *numbers)
+{
+    int i, j;
+
+    for (i = 0; i < sc->wh; i++) {
+        sc->squares[i].nplacements = 0;
+        sc->squares[i].number = numbers[sc->squares[i].index];
+    }
+
+    for (i = 0; i < sc->pc; i++) {
+        struct solver_placement *p = &sc->placements[i];
+        int di = DINDEX(p->squares[0]->number, p->squares[1]->number);
+        p->domino = &sc->dominoes[di];
+    }
+
+    for (i = 0; i < sc->dc; i++)
+        sc->dominoes[i].nplacements = 0;
+    for (i = 0; i < sc->pc; i++)
+        sc->placements[i].domino->nplacements++;
+    for (i = j = 0; i < sc->dc; i++) {
+        sc->dominoes[i].placements = sc->domino_placement_lists + j;
+        j += sc->dominoes[i].nplacements;
+        sc->dominoes[i].nplacements = 0;
+    }
+    for (i = 0; i < sc->pc; i++) {
+        struct solver_placement *p = &sc->placements[i];
+        p->dpi = p->domino->nplacements;
+        p->domino->placements[p->domino->nplacements++] = p;
+        p->active = true;
+    }
+
+    for (i = 0; i < sc->wh; i++)
+        sc->squares[i].nplacements = 0;
+    for (i = 0; i < sc->pc; i++) {
+        struct solver_placement *p = &sc->placements[i];
+        for (j = 0; j < 2; j++) {
+            struct solver_square *sq = p->squares[j];
+            p->spi[j] = sq->nplacements;
+            sq->placements[sq->nplacements++] = p;
+        }
+    }
+
+    sc->max_diff_used = DIFF_TRIVIAL;
+    sc->squares_by_number_initialised = false;
+}
+
+/* Given two placements p,q that overlap, returns si such that
+ * p->squares[si] is the square also in q */
+static int common_square_index(struct solver_placement *p,
+                               struct solver_placement *q)
+{
+    return (p->squares[0] == q->squares[0] ||
+            p->squares[0] == q->squares[1]) ? 0 : 1;
+}
+
+/* Sort function used to set up squares_by_number */
+static int squares_by_number_cmpfn(const void *av, const void *bv)
+{
+    struct solver_square *a = *(struct solver_square *const *)av;
+    struct solver_square *b = *(struct solver_square *const *)bv;
+    return (a->number < b->number ? -1 : a->number > b->number ? +1 :
+            a->index  < b->index  ? -1 : a->index  > b->index  ? +1 : 0);
+}
+
+static void rule_out_placement(
+    struct solver_scratch *sc, struct solver_placement *p)
+{
+    struct solver_domino *d = p->domino;
+    int i, j, si;
+
+#ifdef SOLVER_DIAGNOSTICS
+    if (solver_diagnostics)
+        printf("  ruling out placement %s for domino %s\n", p->name, d->name);
+#endif
+
+    p->active = false;
+
+    i = p->dpi;
+    assert(d->placements[i] == p);
+    if (--d->nplacements != i) {
+        d->placements[i] = d->placements[d->nplacements];
+        d->placements[i]->dpi = i;
+    }
+
+    for (si = 0; si < 2; si++) {
+        struct solver_square *sq = p->squares[si];
+        i = p->spi[si];
+        assert(sq->placements[i] == p);
+        if (--sq->nplacements != i) {
+            sq->placements[i] = sq->placements[sq->nplacements];
+            j = (sq->placements[i]->squares[0] == sq ? 0 : 1);
+            sq->placements[i]->spi[j] = i;
+        }
+    }
 }
 
 /*
- * Returns 0, 1 or 2 for number of solutions. 2 means `any number
- * more than one', or more accurately `we were unable to prove
- * there was only one'.
- * 
- * Outputs in a `placements' array, indexed the same way as the one
- * within this function (see below); entries in there are <0 for a
- * placement ruled out, 0 for an uncertain placement, and 1 for a
- * definite one.
+ * If a domino has only one placement remaining, rule out all other
+ * placements that overlap it.
  */
-static int solver(int w, int h, int n, int *grid, int *output)
+static bool deduce_domino_single_placement(struct solver_scratch *sc, int di)
 {
-    int wh = w*h, dc = DCOUNT(n);
-    int *placements, *heads;
-    int i, j, x, y, ret;
+    struct solver_domino *d = &sc->dominoes[di];
+    struct solver_placement *p, *q;
+    int oi;
+    bool done_something = false;
+
+    if (d->nplacements != 1)
+        return false;
+    p = d->placements[0];
+
+    for (oi = 0; oi < p->noverlaps; oi++) {
+        q = p->overlaps[oi];
+        if (q->active) {
+            if (!done_something) {
+                done_something = true;
+#ifdef SOLVER_DIAGNOSTICS
+                if (solver_diagnostics)
+                    printf("domino %s has unique placement %s\n",
+                           d->name, p->name);
+#endif
+            }
+            rule_out_placement(sc, q);
+        }
+    }
+
+    return done_something;
+}
+
+/*
+ * If a square has only one placement remaining, rule out all other
+ * placements of its domino.
+ */
+static bool deduce_square_single_placement(struct solver_scratch *sc, int si)
+{
+    struct solver_square *sq = &sc->squares[si];
+    struct solver_placement *p;
+    struct solver_domino *d;
+
+    if (sq->nplacements != 1)
+        return false;
+    p = sq->placements[0];
+    d = p->domino;
+
+    if (d->nplacements <= 1)
+        return false;   /* we already knew everything this would tell us */
+
+#ifdef SOLVER_DIAGNOSTICS
+    if (solver_diagnostics)
+        printf("square %s has unique placement %s (domino %s)\n",
+               sq->name, p->name, p->domino->name);
+#endif
+
+    while (d->nplacements > 1)
+        rule_out_placement(
+            sc, d->placements[0] == p ? d->placements[1] : d->placements[0]);
+
+    return true;
+}
+
+/*
+ * If all placements for a square involve the same domino, rule out
+ * all other placements of that domino.
+ */
+static bool deduce_square_single_domino(struct solver_scratch *sc, int si)
+{
+    struct solver_square *sq = &sc->squares[si];
+    struct solver_domino *d;
+    int i;
 
     /*
-     * This array has one entry for every possible domino
-     * placement. Vertical placements are indexed by their top
-     * half, at (y*w+x)*2; horizontal placements are indexed by
-     * their left half at (y*w+x)*2+1.
-     * 
-     * This array is used to link domino placements together into
-     * linked lists, so that we can track all the possible
-     * placements of each different domino. It's also used as a
-     * quick means of looking up an individual placement to see
-     * whether we still think it's possible. Actual values stored
-     * in this array are -2 (placement not possible at all), -1
-     * (end of list), or the array index of the next item.
-     * 
-     * Oh, and -3 for `not even valid', used for array indices
-     * which don't even represent a plausible placement.
+     * We only bother with this if the square has at least _two_
+     * placements. If it only has one, then a simpler deduction will
+     * have handled it already, or will do so the next time round the
+     * main solver loop - and we should let the simpler deduction do
+     * it, because that will give a less overblown diagnostic.
      */
-    placements = snewn(2*wh, int);
-    for (i = 0; i < 2*wh; i++)
-        placements[i] = -3;            /* not even valid */
+    if (sq->nplacements < 2)
+        return false;
+
+    d = sq->placements[0]->domino;
+    for (i = 1; i < sq->nplacements; i++)
+        if (sq->placements[i]->domino != d)
+            return false;              /* not all the same domino */
+
+    if (d->nplacements <= sq->nplacements)
+        return false;       /* no other placements of d to rule out */
+
+#ifdef SOLVER_DIAGNOSTICS
+    if (solver_diagnostics)
+        printf("square %s can only contain domino %s\n", sq->name, d->name);
+#endif
+
+    for (i = d->nplacements; i-- > 0 ;) {
+        struct solver_placement *p = d->placements[i];
+        if (p->squares[0] != sq && p->squares[1] != sq)
+            rule_out_placement(sc, p);
+    }
+
+    return true;
+}
+
+/*
+ * If any placement is overlapped by _all_ possible placements of a
+ * given domino, rule that placement out.
+ */
+static bool deduce_domino_must_overlap(struct solver_scratch *sc, int di)
+{
+    struct solver_domino *d = &sc->dominoes[di];
+    struct solver_placement *intersection[6], *p;
+    int nintersection = 0;
+    int i, j, k;
 
     /*
-     * This array has one entry for every domino, and it is an
-     * index into `placements' denoting the head of the placement
-     * list for that domino.
+     * As in deduce_square_single_domino, we only bother with this
+     * deduction if the domino has at least two placements.
      */
-    heads = snewn(dc, int);
-    for (i = 0; i < dc; i++)
-        heads[i] = -1;
+    if (d->nplacements < 2)
+        return false;
+
+    /* Initialise our set of overlapped placements with all the active
+     * ones overlapped by placements[0]. */
+    p = d->placements[0];
+    for (i = 0; i < p->noverlaps; i++)
+        if (p->overlaps[i]->active)
+            intersection[nintersection++] = p->overlaps[i];
+
+    /* Now loop over the other placements of d, winnowing that set. */
+    for (j = 1; j < d->nplacements; j++) {
+        int old_n;
+
+        p = d->placements[j];
+
+        old_n = nintersection;
+        nintersection = 0;
+
+        for (k = 0; k < old_n; k++) {
+            for (i = 0; i < p->noverlaps; i++)
+                if (p->overlaps[i] == intersection[k])
+                    goto found;
+            /* If intersection[k] isn't in p->overlaps, exclude it
+             * from our set of placements overlapped by everything */
+            continue;
+          found:
+            intersection[nintersection++] = intersection[k];
+        }
+    }
+
+    if (nintersection == 0)
+        return false;                  /* no new exclusions */
+
+    for (i = 0; i < nintersection; i++) {
+        p = intersection[i];
+
+#ifdef SOLVER_DIAGNOSTICS
+        if (solver_diagnostics) {
+            printf("placement %s of domino %s overlaps all placements "
+                   "of domino %s:", p->name, p->domino->name, d->name);
+            for (j = 0; j < d->nplacements; j++)
+                printf(" %s", d->placements[j]->name);
+            printf("\n");
+        }
+#endif
+        rule_out_placement(sc, p);
+    }
+
+    return true;
+}
+
+/*
+ * If a placement of domino D overlaps the only remaining placement
+ * for some square S which is not also for domino D, then placing D
+ * here would require another copy of it in S, so we can rule it out.
+ */
+static bool deduce_local_duplicate(struct solver_scratch *sc, int pi)
+{
+    struct solver_placement *p = &sc->placements[pi];
+    struct solver_domino *d = p->domino;
+    int i, j;
+
+    if (!p->active)
+        return false;
+
+    for (i = 0; i < p->noverlaps; i++) {
+        struct solver_placement *q = p->overlaps[i];
+        struct solver_square *sq;
+
+        if (!q->active)
+            continue;
+
+        /* Find the square of q that _isn't_ part of p */
+        sq = q->squares[1 - common_square_index(q, p)];
+
+        for (j = 0; j < sq->nplacements; j++)
+            if (sq->placements[j] != q && sq->placements[j]->domino != d)
+                goto no;
+
+        /* If we get here, every possible placement for sq is either q
+         * itself, or another copy of d. Success! We can rule out p. */
+#ifdef SOLVER_DIAGNOSTICS
+        if (solver_diagnostics) {
+            printf("placement %s of domino %s would force another copy of %s "
+                   "in square %s\n", p->name, d->name, d->name, sq->name);
+        }
+#endif
+
+        rule_out_placement(sc, p);
+        return true;
+
+      no:;
+    }
+
+    return false;
+}
+
+/*
+ * If placement P overlaps one placement for each of two squares S,T
+ * such that all the remaining placements for both S and T are the
+ * same domino D (and none of those placements joins S and T to each
+ * other), then P can't be placed, because it would leave S,T each
+ * having to be a copy of D, i.e. duplicates.
+ */
+static bool deduce_local_duplicate_2(struct solver_scratch *sc, int pi)
+{
+    struct solver_placement *p = &sc->placements[pi];
+    int i, j, k;
+
+    if (!p->active)
+        return false;
 
     /*
-     * Set up the initial possibility lists by scanning the grid.
+     * Iterate over pairs of placements qi,qj overlapping p.
      */
-    for (y = 0; y < h-1; y++)
-        for (x = 0; x < w; x++) {
-            int di = DINDEX(grid[y*w+x], grid[(y+1)*w+x]);
-            placements[(y*w+x)*2] = heads[di];
-            heads[di] = (y*w+x)*2;
+    for (i = 0; i < p->noverlaps; i++) {
+        struct solver_placement *qi = p->overlaps[i];
+        struct solver_square *sqi;
+        struct solver_domino *di = NULL;
+
+        if (!qi->active)
+            continue;
+
+        /* Find the square of qi that _isn't_ part of p */
+        sqi = qi->squares[1 - common_square_index(qi, p)];
+
+        /*
+         * Identify the unique domino involved in all possible
+         * placements of sqi other than qi. If there isn't a unique
+         * one (either too many or too few), move on and try the next
+         * qi.
+         */
+        for (k = 0; k < sqi->nplacements; k++) {
+            struct solver_placement *pk = sqi->placements[k];
+            if (sqi->placements[k] == qi)
+                continue;              /* not counting qi itself */
+            if (!di)
+                di = pk->domino;
+            else if (di != pk->domino)
+                goto done_qi;
         }
-    for (y = 0; y < h; y++)
-        for (x = 0; x < w-1; x++) {
-            int di = DINDEX(grid[y*w+x], grid[y*w+(x+1)]);
-            placements[(y*w+x)*2+1] = heads[di];
-            heads[di] = (y*w+x)*2+1;
+        if (!di)
+            goto done_qi;
+
+        /*
+         * Now find an appropriate qj != qi.
+         */
+        for (j = 0; j < p->noverlaps; j++) {
+            struct solver_placement *qj = p->overlaps[j];
+            struct solver_square *sqj;
+            bool found_di = false;
+
+            if (j == i || !qj->active)
+                continue;
+
+            sqj = qj->squares[1 - common_square_index(qj, p)];
+
+            /*
+             * As above, we want the same domino di to be the only one
+             * sqj can be if placement qj is ruled out. But also we
+             * need no placement of sqj to overlap sqi.
+             */
+            for (k = 0; k < sqj->nplacements; k++) {
+                struct solver_placement *pk = sqj->placements[k];
+                if (pk == qj)
+                    continue;          /* not counting qj itself */
+                if (pk->domino != di)
+                    goto done_qj;      /* found a different domino */
+                if (pk->squares[0] == sqi || pk->squares[1] == sqi)
+                    goto done_qj; /* sqi,sqj can be joined to each other */
+                found_di = true;
+            }
+            if (!found_di)
+                goto done_qj;
+
+            /* If we get here, then every placement for either of sqi
+             * and sqj is a copy of di, except for the ones that
+             * overlap p. Success! We can rule out p. */
+#ifdef SOLVER_DIAGNOSTICS
+            if (solver_diagnostics) {
+                printf("placement %s of domino %s would force squares "
+                       "%s and %s to both be domino %s\n",
+                       p->name, p->domino->name,
+                       sqi->name, sqj->name, di->name);
+            }
+#endif
+            rule_out_placement(sc, p);
+            return true;
+
+          done_qj:;
         }
 
+      done_qi:;
+    }
+
+    return false;
+}
+
+struct parity_findloop_ctx {
+    struct solver_scratch *sc;
+    struct solver_square *sq;
+    int i;
+};
+
+int parity_neighbour(int vertex, void *vctx)
+{
+    struct parity_findloop_ctx *ctx = (struct parity_findloop_ctx *)vctx;
+    struct solver_placement *p;
+
+    if (vertex >= 0) {
+        ctx->sq = &ctx->sc->squares[vertex];
+        ctx->i = 0;
+    } else {
+        assert(ctx->sq);
+    }
+
+    if (ctx->i >= ctx->sq->nplacements) {
+        ctx->sq = NULL;
+        return -1;
+    }
+
+    p = ctx->sq->placements[ctx->i++];
+    return p->squares[0]->index + p->squares[1]->index - ctx->sq->index;
+}
+
+/*
+ * Look for dominoes whose placement would disconnect the unfilled
+ * area of the grid into pieces with odd area. Such a domino can't be
+ * placed, because then the area on each side of it would be
+ * untileable.
+ */
+static bool deduce_parity(struct solver_scratch *sc)
+{
+    struct parity_findloop_ctx pflctx;
+    bool done_something = false;
+    int pi;
+
+    /*
+     * Run findloop, aka Tarjan's bridge-finding algorithm, on the
+     * graph whose vertices are squares, with two vertices separated
+     * by an edge iff some not-yet-ruled-out domino placement covers
+     * them both. (So each edge itself corresponds to a domino
+     * placement.)
+     *
+     * The effect is that any bridge in this graph is a domino whose
+     * placement would separate two previously connected areas of the
+     * unfilled squares of the grid.
+     *
+     * Placing that domino would not just disconnect those areas from
+     * each other, but also use up one square of each. So if we want
+     * to avoid leaving two odd areas after placing the domino, it
+     * follows that we want to avoid the bridge having an _even_
+     * number of vertices on each side.
+     */
+    pflctx.sc = sc;
+    findloop_run(sc->fls, sc->wh, parity_neighbour, &pflctx);
+
+    for (pi = 0; pi < sc->pc; pi++) {
+        struct solver_placement *p = &sc->placements[pi];
+        int size0, size1;
+
+        if (!p->active)
+            continue;
+        if (!findloop_is_bridge(
+                sc->fls, p->squares[0]->index, p->squares[1]->index,
+                &size0, &size1))
+            continue;
+        /* To make a deduction, size0 and size1 must both be even,
+         * i.e. after placing this domino decrements each by 1 they
+         * would both become odd and untileable areas. */
+        if ((size0 | size1) & 1)
+            continue;
+
 #ifdef SOLVER_DIAGNOSTICS
-    printf("before solver:\n");
-    for (i = 0; i <= n; i++)
-        for (j = 0; j <= i; j++) {
-            int k, m;
-            m = 0;
-            printf("%2d [%d %d]:", DINDEX(i, j), i, j);
-            for (k = heads[DINDEX(i,j)]; k >= 0; k = placements[k])
-                printf(" %3d [%d,%d,%c]", k, k/2%w, k/2/w, k%2?'h':'v');
-            printf("\n");
+        if (solver_diagnostics) {
+            printf("placement %s of domino %s would create two odd-sized "
+                   "areas\n", p->name, p->domino->name);
         }
 #endif
+        rule_out_placement(sc, p);
+        done_something = true;
+    }
 
-    while (1) {
-        bool done_something = false;
+    return done_something;
+}
 
+/*
+ * Try to find a set of squares all containing the same number, such
+ * that the set of possible dominoes for all the squares in that set
+ * is small enough to let us rule out placements of those dominoes
+ * elsewhere.
+ */
+static bool deduce_set(struct solver_scratch *sc, bool doubles)
+{
+    struct solver_square **sqs, **sqp, **sqe;
+    int num, nsq, i, j;
+    unsigned long domino_sets[16], adjacent[16];
+    struct solver_domino *ds[16];
+    bool done_something = false;
+
+    if (!sc->squares_by_number)
+        sc->squares_by_number = snewn(sc->wh, struct solver_square *);
+    if (!sc->wh_scratch)
+        sc->wh_scratch = snewn(sc->wh, int);
+
+    if (!sc->squares_by_number_initialised) {
         /*
-         * For each domino, look at its possible placements, and
-         * for each placement consider the placements (of any
-         * domino) it overlaps. Any placement overlapped by all
-         * placements of this domino can be ruled out.
-         * 
-         * Each domino placement overlaps only six others, so we
-         * need not do serious set theory to work this out.
+         * If this is the first call to this function for a given
+         * grid, start by sorting the squares by their containing
+         * number.
          */
-        for (i = 0; i < dc; i++) {
-            int permset[6], permlen = 0, p;
-            
+        for (i = 0; i < sc->wh; i++)
+            sc->squares_by_number[i] = &sc->squares[i];
+        qsort(sc->squares_by_number, sc->wh, sizeof(*sc->squares_by_number),
+              squares_by_number_cmpfn);
+    }
 
-            if (heads[i] == -1) {      /* no placement for this domino */
-                ret = 0;               /* therefore puzzle is impossible */
-                goto done;
-            }
-            for (j = heads[i]; j >= 0; j = placements[j]) {
-                assert(placements[j] != -2);
+    sqp = sc->squares_by_number;
+    sqe = sc->squares_by_number + sc->wh;
+    for (num = 0; num <= sc->n; num++) {
+        unsigned long squares;
+        unsigned long squares_done;
 
-                if (j == heads[i]) {
-                    permlen = find_overlaps(w, h, j, permset);
-                } else {
-                    int tempset[6], templen, m, n, k;
+        /* Find the bounds of the subinterval of squares_by_number
+         * containing squares with this particular number. */
+        sqs = sqp;
+        while (sqp < sqe && (*sqp)->number == num)
+            sqp++;
+        nsq = sqp - sqs;
+
+        /*
+         * Now sqs[0], ..., sqs[nsq-1] are the squares containing 'num'.
+         */
+
+        if (nsq > lenof(domino_sets)) {
+            /*
+             * Abort this analysis if we're trying to enumerate all
+             * the subsets of a too-large base set.
+             *
+             * This _shouldn't_ happen, at the time of writing this
+             * code, because the largest puzzle we support is only
+             * supposed to have 10 instances of each number, and part
+             * of our input grid validation checks that each number
+             * does appear the right number of times. But just in case
+             * weird test input makes its way to this function, or the
+             * puzzle sizes are expanded later, it's easy enough to
+             * just rule out doing this analysis for overlarge sets of
+             * numbers.
+             */
+            continue;
+        }
+
+        /*
+         * Index the squares in wh_scratch, which we're using as a
+         * lookup table to map the official index of a square back to
+         * its value in our local indexing scheme.
+         */
+        for (i = 0; i < nsq; i++)
+            sc->wh_scratch[sqs[i]->index] = i;
+
+        /*
+         * For each square, make a bit mask of the dominoes that can
+         * overlap it, by finding the number at the other end of each
+         * one.
+         *
+         * Also, for each square, make a bit mask of other squares in
+         * the current list that might occupy the _same_ domino
+         * (because a possible placement of a double overlaps both).
+         * We'll need that for evaluating whether sets are properly
+         * exhaustive.
+         */
+        for (i = 0; i < nsq; i++)
+            adjacent[i] = 0;
 
-                    templen = find_overlaps(w, h, j, tempset);
+        for (i = 0; i < nsq; i++) {
+            struct solver_square *sq = sqs[i];
+            unsigned long mask = 0;
 
+            for (j = 0; j < sq->nplacements; j++) {
+                struct solver_placement *p = sq->placements[j];
+                int othernum = p->domino->lo + p->domino->hi - num;
+                mask |= 1UL << othernum;
+                ds[othernum] = p->domino; /* so we can find them later */
+
+                if (othernum == num) {
                     /*
-                     * Pathetically primitive set intersection
-                     * algorithm, which I'm only getting away with
-                     * because I know my sets are bounded by a very
-                     * small size.
+                     * Special case: this is a double, so it gives
+                     * rise to entries in adjacent[].
                      */
-                    for (m = n = 0; m < permlen; m++) {
-                        for (k = 0; k < templen; k++)
-                            if (tempset[k] == permset[m])
-                                break;
-                        if (k < templen)
-                            permset[n++] = permset[m];
-                    }
-                    permlen = n;
+                    int i2 = sc->wh_scratch[p->squares[0]->index +
+                                            p->squares[1]->index - sq->index];
+                    adjacent[i] |= 1UL << i2;
+                    adjacent[i2] |= 1UL << i;
                 }
             }
-            for (p = 0; p < permlen; p++) {
-                j = permset[p];
-                if (placements[j] != -2) {
-                    int p1, p2, di;
 
-                    done_something = true;
+            domino_sets[i] = mask;
+
+        }
+
+        squares_done = 0;
+
+        for (squares = 0; squares < (1UL << nsq); squares++) {
+            unsigned long dominoes = 0;
+            int bitpos, nsquares, ndominoes;
+            bool got_adj_squares = false;
+            bool reported = false;
+            bool rule_out_nondoubles;
+            int min_nused_for_double;
+#ifdef SOLVER_DIAGNOSTICS
+            const char *rule_out_text;
+#endif
+
+            /*
+             * We don't do set analysis on the same square of the grid
+             * more than once in this loop. Otherwise you generate
+             * pointlessly overcomplicated diagnostics for simpler
+             * follow-up deductions. For example, suppose squares
+             * {A,B} must go with dominoes {X,Y}. So you rule out X,Y
+             * elsewhere, and then it turns out square C (from which
+             * one of those was eliminated) has only one remaining
+             * possibility Z. What you _don't_ want to do is
+             * triumphantly report a second case of set elimination
+             * where you say 'And also, squares {A,B,C} have to be
+             * {X,Y,Z}!' You'd prefer to give 'now C has to be Z' as a
+             * separate deduction later, more simply phrased.
+             */
+            if (squares & squares_done)
+                continue;
+
+            nsquares = 0;
+
+            /* Make the set of dominoes that these squares can inhabit. */
+            for (bitpos = 0; bitpos < nsq; bitpos++) {
+                if (!(1 & (squares >> bitpos)))
+                    continue;          /* this bit isn't set in the mask */
+
+                if (adjacent[bitpos] & squares)
+                    got_adj_squares = true;
 
+                dominoes |= domino_sets[bitpos];
+                nsquares++;
+            }
+
+            /* Count them. */
+            ndominoes = 0;
+            for (bitpos = 0; bitpos < nsq; bitpos++)
+                ndominoes += 1 & (dominoes >> bitpos);
+
+            /*
+             * Do the two sets have the right relative size?
+             */
+            if (!got_adj_squares) {
+                /*
+                 * The normal case, in which every possible domino
+                 * placement involves at most _one_ of these squares.
+                 *
+                 * This is exactly analogous to the set analysis
+                 * deductions in many other puzzles: if our N squares
+                 * between them have to account for N distinct
+                 * dominoes, with exactly one of those dominoes to
+                 * each square, then all those dominoes correspond to
+                 * all those squares and we can rule out any
+                 * placements of the same dominoes appearing
+                 * elsewhere.
+                 */
+                if (ndominoes != nsquares)
+                    continue;
+                rule_out_nondoubles = true;
+                min_nused_for_double = 1;
+#ifdef SOLVER_DIAGNOSTICS
+                rule_out_text = "all of them elsewhere";
+#endif
+            } else {
+                if (!doubles)
+                    continue;          /* not at this difficulty level */
+
+                /*
+                 * But in Dominosa, there's a special case if _two_
+                 * squares in this set can possibly both be covered by
+                 * the same double domino. (I.e. if they are adjacent,
+                 * and moreover, the double-domino placement
+                 * containing both is not yet ruled out.)
+                 *
+                 * In that situation, the simple argument doesn't hold
+                 * up, because the N squares might be covered by N-1
+                 * dominoes - or, put another way, if you list the
+                 * containing domino for each of the squares, they
+                 * might not be all distinct.
+                 *
+                 * In that situation, we can still do something, but
+                 * the details vary, and there are two further cases.
+                 */
+                if (ndominoes == nsquares-1) {
                     /*
-                     * Rule out this placement. First find what
-                     * domino it is...
+                     * Suppose there is one _more_ square in our set
+                     * than there are dominoes it can involve. For
+                     * example, suppose we had four '0' squares which
+                     * between them could contain only the 0-0, 0-1
+                     * and 0-2 dominoes.
+                     *
+                     * Then that can only work at all if the 0-0
+                     * covers two of those squares - and in that
+                     * situation that _must_ be what's happened.
+                     *
+                     * So we can rule out the 0-1 and 0-2 dominoes (in
+                     * this example) in any placement that doesn't use
+                     * one of the squares in this set. And we can rule
+                     * out a placement of the 0-0 even if it uses
+                     * _one_ square from this set: in this situation,
+                     * we have to insist on it using _two_.
                      */
-                    p1 = j / 2;
-                    p2 = (j & 1) ? p1 + 1 : p1 + w;
-                    di = DINDEX(grid[p1], grid[p2]);
+                    rule_out_nondoubles = true;
+                    min_nused_for_double = 2;
 #ifdef SOLVER_DIAGNOSTICS
-                    printf("considering domino %d: ruling out placement %d"
-                           " for %d\n", i, j, di);
+                    rule_out_text = "all of them elsewhere "
+                        "(including the double if it fails to use both)";
 #endif
-
+                } else if (ndominoes == nsquares) {
+                    /*
+                     * A restricted form of the deduction is still
+                     * possible if we have the same number of dominoes
+                     * as squares.
+                     *
+                     * If we have _three_ '0' squares none of which
+                     * can be any domino other than 0-0, 0-1 and 0-2,
+                     * and there's still a possibility of an 0-0
+                     * domino using up two of them, then we can't rule
+                     * out 0-1 or 0-2 anywhere else, because it's
+                     * possible that these three squares only use two
+                     * of the dominoes between them.
+                     *
+                     * But we _can_ rule out the double 0-0, in any
+                     * placement that uses _none_ of our three
+                     * squares. Because we do know that _at least one_
+                     * of our squares must be involved in the 0-0, or
+                     * else the three of them would only have the
+                     * other two dominoes left.
+                     */
+                    rule_out_nondoubles = false;
+                    min_nused_for_double = 1;
+#ifdef SOLVER_DIAGNOSTICS
+                    rule_out_text = "the double elsewhere";
+#endif
+                } else {
                     /*
-                     * ... then walk that domino's placement list,
-                     * removing this placement when we find it.
+                     * If none of those cases has happened, then our
+                     * set admits no deductions at all.
                      */
-                    if (heads[di] == j)
-                        heads[di] = placements[j];
-                    else {
-                        int k = heads[di];
-                        while (placements[k] != -1 && placements[k] != j)
-                            k = placements[k];
-                        assert(placements[k] == j);
-                        placements[k] = placements[j];
+                    continue;
+                }
+            }
+
+            /* Skip sets of size 1, or whose complement has size 1.
+             * Those can be handled by a simpler analysis, and should
+             * be, for more sensible solver diagnostics. */
+            if (ndominoes <= 1 || ndominoes >= nsq-1)
+                continue;
+
+            /*
+             * We've found a set! That means we can rule out any
+             * placement of any domino in that set which would leave
+             * the squares in the set with too few dominoes between
+             * them.
+             *
+             * We may or may not actually end up ruling anything out
+             * here. But even if we don't, we should record that these
+             * squares form a self-contained set, so that we don't
+             * pointlessly report a superset of them later which could
+             * instead be reported as just the other ones.
+             *
+             * Or rather, we do that for the main cases that let us
+             * rule out lots of dominoes. We only do this with the
+             * borderline case where we can only rule out a double if
+             * we _actually_ rule something out. Otherwise we'll never
+             * even _find_ a larger set with the same number of
+             * dominoes!
+             */
+            if (rule_out_nondoubles)
+                squares_done |= squares;
+
+            for (bitpos = 0; bitpos < nsq; bitpos++) {
+                struct solver_domino *d;
+
+                if (!(1 & (dominoes >> bitpos)))
+                    continue;
+                d = ds[bitpos];
+
+                for (i = d->nplacements; i-- > 0 ;) {
+                    struct solver_placement *p = d->placements[i];
+                    int si, nused;
+
+                    /* Count how many of our squares this placement uses. */
+                    for (si = nused = 0; si < 2; si++) {
+                        struct solver_square *sq2 = p->squares[si];
+                        if (sq2->number == num &&
+                            (1 & (squares >> sc->wh_scratch[sq2->index])))
+                            nused++;
+                    }
+
+                    /* See if that's too many to rule it out. */
+                    if (d->lo == d->hi) {
+                        if (nused >= min_nused_for_double)
+                            continue;
+                    } else {
+                        if (nused > 0 || !rule_out_nondoubles)
+                            continue;
+                    }
+
+                    if (!reported) {
+                        reported = true;
+                        done_something = true;
+
+                        /* In case we didn't do this above */
+                        squares_done |= squares;
+
+#ifdef SOLVER_DIAGNOSTICS
+                        if (solver_diagnostics) {
+                            int b;
+                            const char *sep;
+                            printf("squares {");
+                            for (sep = "", b = 0; b < nsq; b++)
+                                if (1 & (squares >> b)) {
+                                    printf("%s%s", sep, sqs[b]->name);
+                                    sep = ",";
+                                }
+                            printf("} can contain only dominoes {");
+                            for (sep = "", b = 0; b < nsq; b++)
+                                if (1 & (dominoes >> b)) {
+                                    printf("%s%s", sep, ds[b]->name);
+                                    sep = ",";
+                                }
+                            printf("}, so rule out %s", rule_out_text);
+                            printf("\n");
+                        }
+#endif
                     }
-                    placements[j] = -2;
+
+                    rule_out_placement(sc, p);
                 }
             }
         }
 
+    }
+
+    return done_something;
+}
+
+static int forcing_chain_dup_cmp(const void *av, const void *bv, void *ctx)
+{
+    struct solver_scratch *sc = (struct solver_scratch *)ctx;
+    int a = *(const int *)av, b = *(const int *)bv;
+    int ac, bc;
+
+    ac = sc->pc_scratch[a];
+    bc = sc->pc_scratch[b];
+    if (ac != bc) return ac > bc ? +1 : -1;
+
+    ac = sc->placements[a].domino->index;
+    bc = sc->placements[b].domino->index;
+    if (ac != bc) return ac > bc ? +1 : -1;
+
+    return 0;
+}
+
+static int forcing_chain_sq_cmp(const void *av, const void *bv, void *ctx)
+{
+    struct solver_scratch *sc = (struct solver_scratch *)ctx;
+    int a = *(const int *)av, b = *(const int *)bv;
+    int ac, bc;
+
+    ac = sc->placements[a].domino->index;
+    bc = sc->placements[b].domino->index;
+    if (ac != bc) return ac > bc ? +1 : -1;
+
+    ac = sc->pc_scratch[a];
+    bc = sc->pc_scratch[b];
+    if (ac != bc) return ac > bc ? +1 : -1;
+
+    return 0;
+}
+
+static bool deduce_forcing_chain(struct solver_scratch *sc)
+{
+    int si, pi, di, j, k, m;
+    bool done_something = false;
+
+    if (!sc->wh_scratch)
+        sc->wh_scratch = snewn(sc->wh, int);
+    if (!sc->pc_scratch)
+        sc->pc_scratch = snewn(sc->pc, int);
+    if (!sc->pc_scratch2)
+        sc->pc_scratch2 = snewn(sc->pc, int);
+    if (!sc->dc_scratch)
+        sc->dc_scratch = snewn(sc->dc, int);
+
+    /*
+     * Start by identifying chains of placements which must all occur
+     * together if any of them occurs. We do this by making
+     * pc_scratch2 an edsf binding the placements into an equivalence
+     * class for each entire forcing chain, with the two possible sets
+     * of dominoes for the chain listed as inverses.
+     */
+    dsf_init(sc->pc_scratch2, sc->pc);
+    for (si = 0; si < sc->wh; si++) {
+        struct solver_square *sq = &sc->squares[si];
+        if (sq->nplacements == 2)
+            edsf_merge(sc->pc_scratch2,
+                       sq->placements[0]->index,
+                       sq->placements[1]->index, true);
+    }
+    /*
+     * Now read out the whole dsf into pc_scratch, flattening its
+     * structured data into a simple integer id per chain of dominoes
+     * that must occur together.
+     *
+     * The integer ids have the property that any two that differ only
+     * in the lowest bit (i.e. of the form {2n,2n+1}) represent
+     * complementary chains, each of which rules out the other.
+     */
+    for (pi = 0; pi < sc->pc; pi++) {
+        bool inv;
+        int c = edsf_canonify(sc->pc_scratch2, pi, &inv);
+        sc->pc_scratch[pi] = c * 2 + (inv ? 1 : 0);
+    }
+
+    /*
+     * Identify chains that contain a duplicate domino, and rule them
+     * out. We do this by making a list of the placement indices in
+     * pc_scratch2, sorted by (chain id, domino id), so that dupes
+     * become adjacent.
+     */
+    for (pi = 0; pi < sc->pc; pi++)
+        sc->pc_scratch2[pi] = pi;
+    arraysort(sc->pc_scratch2, sc->pc, forcing_chain_dup_cmp, sc);
+
+    for (j = 0; j < sc->pc ;) {
+        struct solver_domino *duplicated_domino = NULL;
+
         /*
-         * For each square, look at the available placements
-         * involving that square. If all of them are for the same
-         * domino, then rule out any placements for that domino
-         * _not_ involving this square.
+         * This loop iterates once per contiguous segment of the same
+         * value in pc_scratch2, i.e. once per chain.
          */
-        for (i = 0; i < wh; i++) {
-            int list[4], k, n, adi;
-
-            x = i % w;
-            y = i / w;
-
-            j = 0;
-            if (x > 0)
-                list[j++] = 2*(i-1)+1;
-            if (x+1 < w)
-                list[j++] = 2*i+1;
-            if (y > 0)
-                list[j++] = 2*(i-w);
-            if (y+1 < h)
-                list[j++] = 2*i;
-
-            for (n = k = 0; k < j; k++)
-                if (placements[list[k]] >= -1)
-                    list[n++] = list[k];
-
-            adi = -1;
-
-            for (j = 0; j < n; j++) {
-                int p1, p2, di;
-                k = list[j];
-
-                p1 = k / 2;
-                p2 = (k & 1) ? p1 + 1 : p1 + w;
-                di = DINDEX(grid[p1], grid[p2]);
-
-                if (adi == -1)
-                    adi = di;
-                if (adi != di)
-                    break;
+        int ci = sc->pc_scratch[sc->pc_scratch2[j]];
+        int climit, cstart = j;
+        while (j < sc->pc && sc->pc_scratch[sc->pc_scratch2[j]] == ci)
+            j++;
+        climit = j;
+
+        /*
+         * Now look for a duplicate domino within that chain.
+         */
+        for (k = cstart; k + 1 < climit; k++) {
+            struct solver_placement *p = &sc->placements[sc->pc_scratch2[k]];
+            struct solver_placement *q = &sc->placements[sc->pc_scratch2[k+1]];
+            if (p->domino == q->domino) {
+                duplicated_domino = p->domino;
+                break;
             }
+        }
 
-            if (j == n) {
-                int nn;
+        if (!duplicated_domino)
+            continue;
 
-                assert(adi >= 0);
-                /*
-                 * We've found something. All viable placements
-                 * involving this square are for domino `adi'. If
-                 * the current placement list for that domino is
-                 * longer than n, reduce it to precisely this
-                 * placement list and we've done something.
-                 */
-                nn = 0;
-                for (k = heads[adi]; k >= 0; k = placements[k])
-                    nn++;
-                if (nn > n) {
-                    done_something = true;
 #ifdef SOLVER_DIAGNOSTICS
-                    printf("considering square %d,%d: reducing placements "
-                           "of domino %d\n", x, y, adi);
+        if (solver_diagnostics) {
+            printf("domino %s occurs more than once in forced chain:",
+                   duplicated_domino->name);
+            for (k = cstart; k < climit; k++)
+                printf(" %s", sc->placements[sc->pc_scratch2[k]].name);
+            printf("\n");
+        }
 #endif
+
+        for (k = cstart; k < climit; k++)
+            rule_out_placement(sc, &sc->placements[sc->pc_scratch2[k]]);
+
+        done_something = true;
+    }
+
+    if (done_something)
+        return true;
+
+    /*
+     * A second way in which a whole forcing chain can be ruled out is
+     * if it contains all the dominoes that can occupy some other
+     * square, so that if the domnioes in the chain were all laid, the
+     * other square would be left without any choices.
+     *
+     * To detect this, we sort the placements again, this time by
+     * (domino index, chain index), so that we can easily find a
+     * sorted list of chains per domino. That allows us to iterate
+     * over the squares and check for a chain id common to all the
+     * placements of that square.
+     */
+    for (pi = 0; pi < sc->pc; pi++)
+        sc->pc_scratch2[pi] = pi;
+    arraysort(sc->pc_scratch2, sc->pc, forcing_chain_sq_cmp, sc);
+
+    /* Store a lookup table of the first entry in pc_scratch2
+     * corresponding to each domino. */
+    for (di = j = 0; j < sc->pc; j++) {
+        while (di <= sc->placements[sc->pc_scratch2[j]].domino->index) {
+            assert(di < sc->dc);
+            sc->dc_scratch[di++] = j;
+        }
+    }
+    assert(di == sc->dc);
+
+    for (si = 0; si < sc->wh; si++) {
+        struct solver_square *sq = &sc->squares[si];
+        int listpos = 0, listsize = 0, listout = 0;
+        int exclude[4];
+        int n_exclude;
+
+        if (sq->nplacements < 2)
+            continue;              /* too simple to be worth trying */
+
+        /*
+         * Start by checking for chains this square can actually form
+         * part of. We won't consider those. (The aim is to find a
+         * completely _different_ square whose placements are all
+         * ruled out by a chain.)
+         */
+        assert(sq->nplacements <= lenof(exclude));
+        for (j = n_exclude = 0; j < sq->nplacements; j++)
+            exclude[n_exclude++] = sc->pc_scratch[sq->placements[j]->index];
+
+        for (j = 0; j < sq->nplacements; j++) {
+            struct solver_domino *d = sq->placements[j]->domino;
+
+            listout = listpos = 0;
+
+            for (k = sc->dc_scratch[d->index];
+                 k < sc->pc && sc->placements[sc->pc_scratch2[k]].domino == d;
+                 k++) {
+                int chain = sc->pc_scratch[sc->pc_scratch2[k]];
+                bool keep;
+
+                if (!sc->placements[sc->pc_scratch2[k]].active)
+                    continue;
+
+                if (j == 0) {
+                    keep = true;
+                } else {
+                    while (listpos < listsize &&
+                           sc->wh_scratch[listpos] < chain)
+                        listpos++;
+                    keep = (listpos < listsize &&
+                            sc->wh_scratch[listpos] == chain);
+                }
+
+                for (m = 0; m < n_exclude; m++)
+                    if (chain == exclude[m])
+                        keep = false;
+
+                if (keep)
+                    sc->wh_scratch[listout++] = chain;
+            }
+
+            listsize = listout;
+            if (listsize == 0)
+                break; /* ruled out all chains; terminate loop early */
+        }
+
+        for (listpos = 0; listpos < listsize; listpos++) {
+            int chain = sc->wh_scratch[listpos];
+
+            /*
+             * We've found a chain we can rule out.
+             */
+#ifdef SOLVER_DIAGNOSTICS
+            if (solver_diagnostics) {
+                printf("all choices for square %s would be ruled out "
+                       "by forced chain:", sq->name);
+                for (pi = 0; pi < sc->pc; pi++)
+                    if (sc->pc_scratch[pi] == chain)
+                        printf(" %s", sc->placements[pi].name);
+                printf("\n");
+            }
+#endif
+
+            for (pi = 0; pi < sc->pc; pi++)
+                if (sc->pc_scratch[pi] == chain)
+                    rule_out_placement(sc, &sc->placements[pi]);
+
+            done_something = true;
+        }
+    }
+
+    /*
+     * Another thing you can do with forcing chains, besides ruling
+     * out a whole one at a time, is to look at each pair of chains
+     * that overlap each other. Each such pair gives you two sets of
+     * domino placements, such that if either set is not placed, then
+     * the other one must be.
+     *
+     * This means that any domino which has a placement in _both_
+     * chains of a pair must occupy one of those two placements, i.e.
+     * we can rule that domino out anywhere else it might appear.
+     */
+    for (di = 0; di < sc->dc; di++) {
+        struct solver_domino *d = &sc->dominoes[di];
+
+        if (d->nplacements <= 2)
+            continue;      /* not enough placements to rule one out */
+
+        for (j = 0; j+1 < d->nplacements; j++) {
+            int ij = d->placements[j]->index;
+            int cj = sc->pc_scratch[ij];
+            for (k = j+1; k < d->nplacements; k++) {
+                int ik = d->placements[k]->index;
+                int ck = sc->pc_scratch[ik];
+                if ((cj ^ ck) == 1) {
                     /*
-                     * Set all other placements on the list to
-                     * impossible.
+                     * Placements j,k of domino d are in complementary
+                     * chains, so we can rule out all the others.
                      */
-                    k = heads[adi];
-                    while (k >= 0) {
-                        int tmp = placements[k];
-                        placements[k] = -2;
-                        k = tmp;
+                    int i;
+
+#ifdef SOLVER_DIAGNOSTICS
+                    if (solver_diagnostics) {
+                        printf("domino %s occurs in both complementary "
+                               "forced chains:", d->name);
+                        for (i = 0; i < sc->pc; i++)
+                            if (sc->pc_scratch[i] == cj)
+                                printf(" %s", sc->placements[i].name);
+                        printf(" and");
+                        for (i = 0; i < sc->pc; i++)
+                            if (sc->pc_scratch[i] == ck)
+                                printf(" %s", sc->placements[i].name);
+                        printf("\n");
                     }
-                    /*
-                     * Set up the new list.
-                     */
-                    heads[adi] = list[0];
-                    for (k = 0; k < n; k++)
-                        placements[list[k]] = (k+1 == n ? -1 : list[k+1]);
+#endif
+
+                    for (i = d->nplacements; i-- > 0 ;)
+                        if (i != j && i != k)
+                            rule_out_placement(sc, d->placements[i]);
+
+                    done_something = true;
+                    goto done_this_domino;
                 }
             }
         }
 
-        if (!done_something)
-            break;
+      done_this_domino:;
     }
 
+    return done_something;
+}
+
+/*
+ * Run the solver until it can't make any more progress.
+ *
+ * Return value is:
+ *   0 = no solution exists (puzzle clues are unsatisfiable)
+ *   1 = unique solution found (success!)
+ *   2 = multiple possibilities remain (puzzle is ambiguous or solver is not
+ *                                      smart enough)
+ */
+static int run_solver(struct solver_scratch *sc, int max_diff_allowed)
+{
+    int di, si, pi;
+    bool done_something;
+
 #ifdef SOLVER_DIAGNOSTICS
-    printf("after solver:\n");
-    for (i = 0; i <= n; i++)
-        for (j = 0; j <= i; j++) {
-            int k, m;
-            m = 0;
-            printf("%2d [%d %d]:", DINDEX(i, j), i, j);
-            for (k = heads[DINDEX(i,j)]; k >= 0; k = placements[k])
-                printf(" %3d [%d,%d,%c]", k, k/2%w, k/2/w, k%2?'h':'v');
+    if (solver_diagnostics) {
+        int di, j;
+        printf("Initial possible placements:\n");
+        for (di = 0; di < sc->dc; di++) {
+            struct solver_domino *d = &sc->dominoes[di];
+            printf("  %s:", d->name);
+            for (j = 0; j < d->nplacements; j++)
+                printf(" %s", d->placements[j]->name);
             printf("\n");
         }
+    }
 #endif
 
-    ret = 1;
-    for (i = 0; i < wh*2; i++) {
-        if (placements[i] == -2) {
-            if (output)
-                output[i] = -1;        /* ruled out */
-        } else if (placements[i] != -3) {
-            int p1, p2, di;
-
-            p1 = i / 2;
-            p2 = (i & 1) ? p1 + 1 : p1 + w;
-            di = DINDEX(grid[p1], grid[p2]);
-
-            if (i == heads[di] && placements[i] == -1) {
-                if (output)
-                    output[i] = 1;     /* certain */
-            } else {
-                if (output)
-                    output[i] = 0;     /* uncertain */
-                ret = 2;
+    do {
+        done_something = false;
+
+        for (di = 0; di < sc->dc; di++)
+            if (deduce_domino_single_placement(sc, di))
+                done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_TRIVIAL);
+            continue;
+        }
+
+        for (si = 0; si < sc->wh; si++)
+            if (deduce_square_single_placement(sc, si))
+                done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_TRIVIAL);
+            continue;
+        }
+
+        if (max_diff_allowed <= DIFF_TRIVIAL)
+            continue;
+
+        for (si = 0; si < sc->wh; si++)
+            if (deduce_square_single_domino(sc, si))
+                done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_BASIC);
+            continue;
+        }
+
+        for (di = 0; di < sc->dc; di++)
+            if (deduce_domino_must_overlap(sc, di))
+                done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_BASIC);
+            continue;
+        }
+
+        for (pi = 0; pi < sc->pc; pi++)
+            if (deduce_local_duplicate(sc, pi))
+                done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_BASIC);
+            continue;
+        }
+
+        for (pi = 0; pi < sc->pc; pi++)
+            if (deduce_local_duplicate_2(sc, pi))
+                done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_BASIC);
+            continue;
+        }
+
+        if (deduce_parity(sc))
+            done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_BASIC);
+            continue;
+        }
+
+        if (max_diff_allowed <= DIFF_BASIC)
+            continue;
+
+        if (deduce_set(sc, false))
+            done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_HARD);
+            continue;
+        }
+
+        if (max_diff_allowed <= DIFF_HARD)
+            continue;
+
+        if (deduce_set(sc, true))
+            done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_EXTREME);
+            continue;
+        }
+
+        if (deduce_forcing_chain(sc))
+            done_something = true;
+        if (done_something) {
+            sc->max_diff_used = max(sc->max_diff_used, DIFF_EXTREME);
+            continue;
+        }
+
+    } while (done_something);
+
+#ifdef SOLVER_DIAGNOSTICS
+    if (solver_diagnostics) {
+        int di, j;
+        printf("Final possible placements:\n");
+        for (di = 0; di < sc->dc; di++) {
+            struct solver_domino *d = &sc->dominoes[di];
+            printf("  %s:", d->name);
+            for (j = 0; j < d->nplacements; j++)
+                printf(" %s", d->placements[j]->name);
+            printf("\n");
+        }
+    }
+#endif
+
+    for (di = 0; di < sc->dc; di++)
+        if (sc->dominoes[di].nplacements == 0)
+            return 0;
+    for (di = 0; di < sc->dc; di++)
+        if (sc->dominoes[di].nplacements > 1)
+            return 2;
+    return 1;
+}
+
+/* ----------------------------------------------------------------------
+ * Functions for generating a candidate puzzle (before we run the
+ * solver to check it's soluble at the right difficulty level).
+ */
+
+struct alloc_val;
+struct alloc_loc;
+
+struct alloc_scratch {
+    /* Game parameters. */
+    int n, w, h, wh, dc;
+
+    /* The domino layout. Indexed by squares in the usual y*w+x raster
+     * order: layout[i] gives the index of the other square in the
+     * same domino as square i. */
+    int *layout;
+
+    /* The output array, containing a number in every square. */
+    int *numbers;
+
+    /* List of domino values (i.e. number pairs), indexed by DINDEX. */
+    struct alloc_val *vals;
+
+    /* List of domino locations, indexed arbitrarily. */
+    struct alloc_loc *locs;
+
+    /* Preallocated scratch spaces. */
+    int *wh_scratch;                   /* size wh */
+    int *wh2_scratch;                  /* size 2*wh */
+};
+
+struct alloc_val {
+    int lo, hi;
+    bool confounder;
+};
+
+struct alloc_loc {
+    int sq[2];
+};
+
+static struct alloc_scratch *alloc_make_scratch(int n)
+{
+    struct alloc_scratch *as = snew(struct alloc_scratch);
+    int lo, hi;
+
+    as->n = n;
+    as->w = n+2;
+    as->h = n+1;
+    as->wh = as->w * as->h;
+    as->dc = DCOUNT(n);
+
+    as->layout = snewn(as->wh, int);
+    as->numbers = snewn(as->wh, int);
+    as->vals = snewn(as->dc, struct alloc_val);
+    as->locs = snewn(as->dc, struct alloc_loc);
+    as->wh_scratch = snewn(as->wh, int);
+    as->wh2_scratch = snewn(as->wh * 2, int);
+
+    for (hi = 0; hi <= n; hi++)
+        for (lo = 0; lo <= hi; lo++) {
+            struct alloc_val *v = &as->vals[DINDEX(hi, lo)];
+            v->lo = lo;
+            v->hi = hi;
+        }
+
+    return as;
+}
+
+static void alloc_free_scratch(struct alloc_scratch *as)
+{
+    sfree(as->layout);
+    sfree(as->numbers);
+    sfree(as->vals);
+    sfree(as->locs);
+    sfree(as->wh_scratch);
+    sfree(as->wh2_scratch);
+    sfree(as);
+}
+
+static void alloc_make_layout(struct alloc_scratch *as, random_state *rs)
+{
+    int i, pos;
+
+    domino_layout_prealloc(as->w, as->h, rs,
+                           as->layout, as->wh_scratch, as->wh2_scratch);
+
+    for (i = pos = 0; i < as->wh; i++) {
+        if (as->layout[i] > i) {
+            struct alloc_loc *loc;
+            assert(pos < as->dc);
+
+            loc = &as->locs[pos++];
+            loc->sq[0] = i;
+            loc->sq[1] = as->layout[i];
+        }
+    }
+    assert(pos == as->dc);
+}
+
+static void alloc_trivial(struct alloc_scratch *as, random_state *rs)
+{
+    int i;
+    for (i = 0; i < as->dc; i++)
+        as->wh_scratch[i] = i;
+    shuffle(as->wh_scratch, as->dc, sizeof(*as->wh_scratch), rs);
+
+    for (i = 0; i < as->dc; i++) {
+        struct alloc_val *val = &as->vals[as->wh_scratch[i]];
+        struct alloc_loc *loc = &as->locs[i];
+        int which_lo = random_upto(rs, 2), which_hi = 1 - which_lo;
+        as->numbers[loc->sq[which_lo]] = val->lo;
+        as->numbers[loc->sq[which_hi]] = val->hi;
+    }
+}
+
+/*
+ * Given a domino location in the form of two square indices, compute
+ * the square indices of the domino location that would lie on one
+ * side of it. Returns false if the location would be outside the
+ * grid, or if it isn't actually a domino in the layout.
+ */
+static bool alloc_find_neighbour(
+    struct alloc_scratch *as, int p0, int p1, int *n0, int *n1)
+{
+    int x0 = p0 % as->w, y0 = p0 / as->w, x1 = p1 % as->w, y1 = p1 / as->w;
+    int dy = y1-y0, dx = x1-x0;
+    int nx0 = x0 + dy, ny0 = y0 - dx, nx1 = x1 + dy, ny1 = y1 - dx;
+    int np0, np1;
+
+    if (!(nx0 >= 0 && nx0 < as->w && ny0 >= 0 && ny0 < as->h &&
+          nx1 >= 1 && nx1 < as->w && ny1 >= 1 && ny1 < as->h))
+        return false;                  /* out of bounds */
+
+    np0 = ny0 * as->w + nx0;
+    np1 = ny1 * as->w + nx1;
+    if (as->layout[np0] != np1)
+        return false;                  /* not a domino */
+
+    *n0 = np0;
+    *n1 = np1;
+    return true;
+}
+
+static bool alloc_try_unique(struct alloc_scratch *as, random_state *rs)
+{
+    int i;
+    for (i = 0; i < as->dc; i++)
+        as->wh_scratch[i] = i;
+    shuffle(as->wh_scratch, as->dc, sizeof(*as->wh_scratch), rs);
+    for (i = 0; i < as->dc; i++)
+        as->wh2_scratch[i] = i;
+    shuffle(as->wh2_scratch, as->dc, sizeof(*as->wh2_scratch), rs);
+
+    for (i = 0; i < as->wh; i++)
+        as->numbers[i] = -1;
+
+    for (i = 0; i < as->dc; i++) {
+        struct alloc_val *val = &as->vals[as->wh_scratch[i]];
+        struct alloc_loc *loc = &as->locs[as->wh2_scratch[i]];
+        int which_lo, which_hi;
+        bool can_lo_0 = true, can_lo_1 = true;
+        int n0, n1;
+
+        /*
+         * This is basically the same strategy as alloc_trivial:
+         * simply iterate through the locations and values in random
+         * relative order and pair them up. But we make sure to avoid
+         * the most common, and also simplest, cause of a non-unique
+         * solution:two dominoes side by side, sharing a number at
+         * opposite ends. Any section of that form automatically leads
+         * to an alternative solution:
+         *
+         *  +-------+         +---+---+
+         *  | 1   2 |         | 1 | 2 |
+         *  +-------+   <->   |   |   |
+         *  | 2   3 |         | 2 | 3 |
+         *  +-------+         +---+---+
+         *
+         * So as we place each domino, we check for a neighbouring
+         * domino on each side, and if there is one, rule out any
+         * placement of _this_ domino that places a number diagonally
+         * opposite the same number in the neighbour.
+         *
+         * Sometimes this can fail completely, if a domino on each
+         * side is already placed and between them they rule out all
+         * placements of this one. But it happens rarely enough that
+         * it's fine to just abort and try the layout again.
+         */
+
+        if (alloc_find_neighbour(as, loc->sq[0], loc->sq[1], &n0, &n1) &&
+            (as->numbers[n0] == val->hi || as->numbers[n1] == val->lo))
+            can_lo_0 = false;
+        if (alloc_find_neighbour(as, loc->sq[1], loc->sq[0], &n0, &n1) &&
+            (as->numbers[n0] == val->hi || as->numbers[n1] == val->lo))
+            can_lo_1 = false;
+
+        if (!can_lo_0 && !can_lo_1)
+            return false;              /* layout failed */
+        else if (can_lo_0 && can_lo_1)
+            which_lo = random_upto(rs, 2);
+        else
+            which_lo = can_lo_0 ? 0 : 1;
+
+        which_hi = 1 - which_lo;
+        as->numbers[loc->sq[which_lo]] = val->lo;
+        as->numbers[loc->sq[which_hi]] = val->hi;
+    }
+
+    return true;
+}
+
+static bool alloc_try_hard(struct alloc_scratch *as, random_state *rs)
+{
+    int i, x, y, hi, lo, vals, locs, confounders_needed;
+    bool ok;
+
+    for (i = 0; i < as->wh; i++)
+        as->numbers[i] = -1;
+
+    /*
+     * Shuffle the location indices.
+     */
+    for (i = 0; i < as->dc; i++)
+        as->wh2_scratch[i] = i;
+    shuffle(as->wh2_scratch, as->dc, sizeof(*as->wh2_scratch), rs);
+
+    /*
+     * Start by randomly placing the double dominoes, to give a
+     * starting instance of every number to try to put other things
+     * next to.
+     */
+    for (i = 0; i <= as->n; i++)
+        as->wh_scratch[i] = DINDEX(i, i);
+    shuffle(as->wh_scratch, i, sizeof(*as->wh_scratch), rs);
+    for (i = 0; i <= as->n; i++) {
+        struct alloc_loc *loc = &as->locs[as->wh2_scratch[i]];
+        as->numbers[loc->sq[0]] = as->numbers[loc->sq[1]] = i;
+    }
+
+    /*
+     * Find all the dominoes that don't yet have a _wrong_ placement
+     * somewhere in the grid.
+     */
+    for (i = 0; i < as->dc; i++)
+        as->vals[i].confounder = false;
+    for (y = 0; y < as->h; y++) {
+        for (x = 0; x < as->w; x++) {
+            int p = y * as->w + x;
+            if (as->numbers[p] == -1)
+                continue;
+
+            if (x+1 < as->w) {
+                int p1 = y * as->w + (x+1);
+                if (as->layout[p] != p1 && as->numbers[p1] != -1)
+                    as->vals[DINDEX(as->numbers[p], as->numbers[p1])]
+                        .confounder = true;
+            }
+            if (y+1 < as->h) {
+                int p1 = (y+1) * as->w + x;
+                if (as->layout[p] != p1 && as->numbers[p1] != -1)
+                    as->vals[DINDEX(as->numbers[p], as->numbers[p1])]
+                        .confounder = true;
             }
         }
     }
 
-    done:
+    for (i = confounders_needed = 0; i < as->dc; i++)
+        if (!as->vals[i].confounder)
+            confounders_needed++;
+
     /*
-     * Free working data.
+     * Make a shuffled list of all the unplaced dominoes, and go
+     * through it trying to find a placement for each one that also
+     * fills in at least one of the needed confounders.
      */
-    sfree(placements);
-    sfree(heads);
+    vals = 0;
+    for (hi = 0; hi <= as->n; hi++)
+        for (lo = 0; lo < hi; lo++)
+            as->wh_scratch[vals++] = DINDEX(hi, lo);
+    shuffle(as->wh_scratch, vals, sizeof(*as->wh_scratch), rs);
+
+    locs = as->dc;
+
+    while (vals > 0) {
+        int valpos, valout, oldvals = vals;
+
+        for (valpos = valout = 0; valpos < vals; valpos++) {
+            int validx = as->wh_scratch[valpos];
+            struct alloc_val *val = &as->vals[validx];
+            struct alloc_loc *loc;
+            int locpos, si, which_lo;
+
+            for (locpos = 0; locpos < locs; locpos++) {
+                int locidx = as->wh2_scratch[locpos];
+                int wi, flip;
+
+                loc = &as->locs[locidx];
+                if (as->numbers[loc->sq[0]] != -1)
+                    continue;              /* this location is already filled */
+
+                flip = random_upto(rs, 2);
+
+                /* Try this location both ways round. */
+                for (wi = 0; wi < 2; wi++) {
+                    int n0, n1;
+
+                    which_lo = wi ^ flip;
+
+                    /* First, do the same check as in alloc_try_unique, to
+                     * avoid making an obviously insoluble puzzle. */
+                    if (alloc_find_neighbour(as, loc->sq[which_lo],
+                                             loc->sq[1-which_lo], &n0, &n1) &&
+                        (as->numbers[n0] == val->hi ||
+                         as->numbers[n1] == val->lo))
+                        break;             /* can't place it this way round */
+
+                    if (confounders_needed == 0)
+                        goto place_ok;
+
+                    /* Look to see if we're adding at least one
+                     * previously absent confounder. */
+                    for (si = 0; si < 2; si++) {
+                        int x = loc->sq[si] % as->w, y = loc->sq[si] / as->w;
+                        int n = (si == which_lo ? val->lo : val->hi);
+                        int d;
+                        for (d = 0; d < 4; d++) {
+                            int dx = d==0 ? +1 : d==2 ? -1 : 0;
+                            int dy = d==1 ? +1 : d==3 ? -1 : 0;
+                            int x1 = x+dx, y1 = y+dy, p1 = y1 * as->w + x1;
+                            if (x1 >= 0 && x1 < as->w &&
+                                y1 >= 0 && y1 < as->h &&
+                                as->numbers[p1] != -1 &&
+                                !(as->vals[DINDEX(n, as->numbers[p1])]
+                                  .confounder)) {
+                                /*
+                                 * Place this domino.
+                                 */
+                                goto place_ok;
+                            }
+                        }
+                    }
+                }
+            }
 
-    return ret;
-}
+            /* If we get here without executing 'goto place_ok', we
+             * didn't find anywhere useful to put this domino. Put it
+             * back on the list for the next pass. */
+            as->wh_scratch[valout++] = validx;
+            continue;
+
+          place_ok:;
+
+            /* We've found a domino to place. Place it, and fill in
+             * all the confounders it adds. */
+            as->numbers[loc->sq[which_lo]] = val->lo;
+            as->numbers[loc->sq[1 - which_lo]] = val->hi;
+
+            for (si = 0; si < 2; si++) {
+                int p = loc->sq[si];
+                int n = as->numbers[p];
+                int x = p % as->w, y = p / as->w;
+                int d;
+                for (d = 0; d < 4; d++) {
+                    int dx = d==0 ? +1 : d==2 ? -1 : 0;
+                    int dy = d==1 ? +1 : d==3 ? -1 : 0;
+                    int x1 = x+dx, y1 = y+dy, p1 = y1 * as->w + x1;
+
+                    if (x1 >= 0 && x1 < as->w && y1 >= 0 && y1 < as->h &&
+                        p1 != loc->sq[1-si] && as->numbers[p1] != -1) {
+                        int di = DINDEX(n, as->numbers[p1]);
+                        if (!as->vals[di].confounder)
+                            confounders_needed--;
+                        as->vals[di].confounder = true;
+                    }
+                }
+            }
+        }
 
-/* ----------------------------------------------------------------------
- * End of solver code.
- */
+        vals = valout;
+
+        if (oldvals == vals)
+            break;
+    }
+
+    ok = true;
+
+    for (i = 0; i < as->dc; i++)
+        if (!as->vals[i].confounder)
+            ok = false;
+    for (i = 0; i < as->wh; i++)
+        if (as->numbers[i] == -1)
+            ok = false;
+
+    return ok;
+}
 
 static char *new_game_desc(const game_params *params, random_state *rs,
                            char **aux, bool interactive)
 {
-    int n = params->n, w = n+2, h = n+1, wh = w*h;
-    int *grid, *grid2, *list;
+    int n = params->n, w = n+2, h = n+1, wh = w*h, diff = params->diff;
+    struct solver_scratch *sc;
+    struct alloc_scratch *as;
     int i, j, k, len;
     char *ret;
 
+#ifndef OMIT_DIFFICULTY_CAP
+    /*
+     * Cap the difficulty level for small puzzles which would
+     * otherwise become impossible to generate.
+     *
+     * Under an #ifndef, to make it easy to remove this cap for the
+     * purpose of re-testing what it ought to be.
+     */
+    if (diff != DIFF_AMBIGUOUS) {
+        if (n == 1 && diff > DIFF_TRIVIAL)
+            diff = DIFF_TRIVIAL;
+        if (n == 2 && diff > DIFF_BASIC)
+            diff = DIFF_BASIC;
+    }
+#endif /* OMIT_DIFFICULTY_CAP */
+
     /*
      * Allocate space in which to lay the grid out.
      */
-    grid = snewn(wh, int);
-    grid2 = snewn(wh, int);
-    list = snewn(2*wh, int);
+    sc = solver_make_scratch(n);
+    as = alloc_make_scratch(n);
 
     /*
      * I haven't been able to think of any particularly clever
@@ -614,91 +2269,75 @@ static char *new_game_desc(const game_params *params, random_state *rs,
      * and 26 respectively, which is a lot more sensible.
      */
 
-    do {
-        domino_layout_prealloc(w, h, rs, grid, grid2, list);
+    while (1) {
+        alloc_make_layout(as, rs);
+
+        if (diff == DIFF_AMBIGUOUS) {
+            /* Just assign numbers to each domino completely at random. */
+            alloc_trivial(as, rs);
+        } else if (diff < DIFF_HARD) {
+            /* Try to rule out the most common case of a non-unique solution */
+            if (!alloc_try_unique(as, rs))
+                continue;
+        } else {
+            /*
+             * For Hard puzzles and above, we'd like there not to be
+             * any easy toehold to start with.
+             *
+             * Mostly, that's arranged by alloc_try_hard, which will
+             * ensure that no domino starts off with only one
+             * potential placement. But a few other deductions
+             * possible at Basic level can still sneak through the
+             * cracks - for example, if the only two placements of one
+             * domino overlap in a square, and you therefore rule out
+             * some other domino that can use that square, you might
+             * then find that _that_ domino now has only one
+             * placement, and you've made a start.
+             *
+             * Of course, the main difficulty-level check will still
+             * guarantee that you have to do a harder deduction
+             * _somewhere_ in the grid. But it's more elegant if
+             * there's nowhere obvious to get started at all.
+             */
+            int di;
+            bool ok;
 
-        /*
-         * Now we have a complete layout covering the whole
-         * rectangle with dominoes. So shuffle the actual domino
-         * values and fill the rectangle with numbers.
-         */
-        k = 0;
-        for (i = 0; i <= params->n; i++)
-            for (j = 0; j <= i; j++) {
-                list[k++] = i;
-                list[k++] = j;
-            }
-        shuffle(list, k/2, 2*sizeof(*list), rs);
-        j = 0;
-        for (i = 0; i < wh; i++)
-            if (grid[i] > i) {
-                /* Optionally flip the domino round. */
-                int flip = -1;
+            if (!alloc_try_hard(as, rs))
+                continue;
 
-                if (params->unique) {
-                    int t1, t2;
-                    /*
-                     * If we're after a unique solution, we can do
-                     * something here to improve the chances. If
-                     * we're placing a domino so that it forms a
-                     * 2x2 rectangle with one we've already placed,
-                     * and if that domino and this one share a
-                     * number, we can try not to put them so that
-                     * the identical numbers are diagonally
-                     * separated, because that automatically causes
-                     * non-uniqueness:
-                     * 
-                     * +---+      +-+-+
-                     * |2 3|      |2|3|
-                     * +---+  ->  | | |
-                     * |4 2|      |4|2|
-                     * +---+      +-+-+
-                     */
-                    t1 = i;
-                    t2 = grid[i];
-                    if (t2 == t1 + w) {  /* this domino is vertical */
-                        if (t1 % w > 0 &&/* and not on the left hand edge */
-                            grid[t1-1] == t2-1 &&/* alongside one to left */
-                            (grid2[t1-1] == list[j] ||   /* and has a number */
-                             grid2[t1-1] == list[j+1] ||   /* in common */
-                             grid2[t2-1] == list[j] ||
-                             grid2[t2-1] == list[j+1])) {
-                            if (grid2[t1-1] == list[j] ||
-                                grid2[t2-1] == list[j+1])
-                                flip = 0;
-                            else
-                                flip = 1;
-                        }
-                    } else {           /* this domino is horizontal */
-                        if (t1 / w > 0 &&/* and not on the top edge */
-                            grid[t1-w] == t2-w &&/* alongside one above */
-                            (grid2[t1-w] == list[j] ||   /* and has a number */
-                             grid2[t1-w] == list[j+1] ||   /* in common */
-                             grid2[t2-w] == list[j] ||
-                             grid2[t2-w] == list[j+1])) {
-                            if (grid2[t1-w] == list[j] ||
-                                grid2[t2-w] == list[j+1])
-                                flip = 0;
-                            else
-                                flip = 1;
-                        }
-                    }
-                }
+            solver_setup_grid(sc, as->numbers);
+            if (run_solver(sc, DIFF_BASIC) < 2)
+                continue;
 
-                if (flip < 0)
-                    flip = random_upto(rs, 2);
+            ok = true;
+            for (di = 0; di < sc->dc; di++)
+                if (sc->dominoes[di].nplacements <= 1) {
+                    ok = false;
+                    break;
+                }
 
-                grid2[i] = list[j + flip];
-                grid2[grid[i]] = list[j + 1 - flip];
-                j += 2;
+            if (!ok) {
+                continue;
             }
-        assert(j == k);
-    } while (params->unique && solver(w, h, n, grid2, NULL) > 1);
+        }
+
+        if (diff != DIFF_AMBIGUOUS) {
+            int solver_result;
+            solver_setup_grid(sc, as->numbers);
+            solver_result = run_solver(sc, diff);
+            if (solver_result > 1)
+                continue; /* puzzle couldn't be solved at this difficulty */
+            if (sc->max_diff_used < diff)
+                continue; /* puzzle _could_ be solved at easier difficulty */
+        }
+
+        break;
+    }
 
 #ifdef GENERATION_DIAGNOSTICS
     for (j = 0; j < h; j++) {
         for (i = 0; i < w; i++) {
-            putchar('0' + grid2[j*w+i]);
+            putchar('0' + as->numbers[j*w+i]);
         }
         putchar('\n');
     }
@@ -738,7 +2377,7 @@ static char *new_game_desc(const game_params *params, random_state *rs,
     ret = snewn(len+1, char);
     j = 0;
     for (i = 0; i < wh; i++) {
-        k = grid2[i];
+        k = as->numbers[i];
         if (k < 10)
             ret[j++] = '0' + k;
         else
@@ -755,7 +2394,7 @@ static char *new_game_desc(const game_params *params, random_state *rs,
         char *auxinfo = snewn(wh+1, char);
 
         for (i = 0; i < wh; i++) {
-            int v = grid[i];
+            int v = as->layout[i];
             auxinfo[i] = (v == i+1 ? 'L' : v == i-1 ? 'R' :
                           v == i+w ? 'T' : v == i-w ? 'B' : '.');
         }
@@ -764,9 +2403,8 @@ static char *new_game_desc(const game_params *params, random_state *rs,
         *aux = auxinfo;
     }
 
-    sfree(list);
-    sfree(grid2);
-    sfree(grid);
+    solver_free_scratch(sc);
+    alloc_free_scratch(as);
 
     return ret;
 }
@@ -898,14 +2536,62 @@ static void free_game(game_state *state)
     sfree(state);
 }
 
+static char *solution_move_string(struct solver_scratch *sc)
+{
+    char *ret;
+    int retlen, retsize;
+    int i, pass;
+
+    /*
+     * First make a pass putting in edges for -1, then make a pass
+     * putting in dominoes for +1.
+     */
+    retsize = 256;
+    ret = snewn(retsize, char);
+    retlen = sprintf(ret, "S");
+
+    for (pass = 0; pass < 2; pass++) {
+        char type = "ED"[pass];
+
+        for (i = 0; i < sc->pc; i++) {
+            struct solver_placement *p = &sc->placements[i];
+            char buf[80];
+            int extra;
+
+            if (pass == 0) {
+                /* Emit a barrier if this placement is ruled out for
+                 * the domino. */
+                if (p->active)
+                    continue;
+            } else {
+                /* Emit a domino if this placement is the only one not
+                 * ruled out. */
+                if (!p->active || p->domino->nplacements > 1)
+                    continue;
+            }
+
+            extra = sprintf(buf, ";%c%d,%d", type,
+                            p->squares[0]->index, p->squares[1]->index);
+
+            if (retlen + extra + 1 >= retsize) {
+                retsize = retlen + extra + 256;
+                ret = sresize(ret, retsize, char);
+            }
+            strcpy(ret + retlen, buf);
+            retlen += extra;
+        }
+    }
+
+    return ret;
+}
+
 static char *solve_game(const game_state *state, const game_state *currstate,
                         const char *aux, const char **error)
 {
     int n = state->params.n, w = n+2, h = n+1, wh = w*h;
-    int *placements;
     char *ret;
     int retlen, retsize;
-    int i, v;
+    int i;
     char buf[80];
     int extra;
 
@@ -931,38 +2617,11 @@ static char *solve_game(const game_state *state, const game_state *currstate,
         }
 
     } else {
-
-        placements = snewn(wh*2, int);
-        for (i = 0; i < wh*2; i++)
-            placements[i] = -3;
-        solver(w, h, n, state->numbers->numbers, placements);
-
-        /*
-         * First make a pass putting in edges for -1, then make a pass
-         * putting in dominoes for +1.
-         */
-        retsize = 256;
-        ret = snewn(retsize, char);
-        retlen = sprintf(ret, "S");
-
-        for (v = -1; v <= +1; v += 2)
-            for (i = 0; i < wh*2; i++)
-                if (placements[i] == v) {
-                    int p1 = i / 2;
-                    int p2 = (i & 1) ? p1+1 : p1+w;
-
-                    extra = sprintf(buf, ";%c%d,%d",
-                                    (int)(v==-1 ? 'E' : 'D'), p1, p2);
-
-                    if (retlen + extra + 1 >= retsize) {
-                        retsize = retlen + extra + 256;
-                        ret = sresize(ret, retsize, char);
-                    }
-                    strcpy(ret + retlen, buf);
-                    retlen += extra;
-                }
-
-        sfree(placements);
+        struct solver_scratch *sc = solver_make_scratch(n);
+        solver_setup_grid(sc, state->numbers->numbers);
+        run_solver(sc, DIFFCOUNT);
+        ret = solution_move_string(sc);
+        solver_free_scratch(sc);
     }
 
     return ret;
@@ -1770,5 +3429,98 @@ const struct game thegame = {
     0,                                 /* flags */
 };
 
+#ifdef STANDALONE_SOLVER
+
+int main(int argc, char **argv)
+{
+    game_params *p;
+    game_state *s, *s2;
+    char *id = NULL, *desc;
+    int maxdiff = DIFFCOUNT;
+    const char *err;
+    bool grade = false, diagnostics = false;
+    struct solver_scratch *sc;
+    int retd;
+
+    while (--argc > 0) {
+        char *p = *++argv;
+        if (!strcmp(p, "-v")) {
+            diagnostics = true;
+        } else if (!strcmp(p, "-g")) {
+            grade = true;
+        } else if (!strncmp(p, "-d", 2) && p[2] && !p[3]) {
+            int i;
+            bool bad = true;
+            for (i = 0; i < lenof(dominosa_diffchars); i++)
+                if (dominosa_diffchars[i] != DIFF_AMBIGUOUS &&
+                    dominosa_diffchars[i] == p[2]) {
+                    bad = false;
+                    maxdiff = i;
+                    break;
+                }
+            if (bad) {
+                fprintf(stderr, "%s: unrecognised difficulty `%c'\n",
+                        argv[0], p[2]);
+                return 1;
+            }
+        } else if (*p == '-') {
+            fprintf(stderr, "%s: unrecognised option `%s'\n", argv[0], p);
+            return 1;
+        } else {
+            id = p;
+        }
+    }
+
+    if (!id) {
+        fprintf(stderr, "usage: %s [-v | -g] <game_id>\n", argv[0]);
+        return 1;
+    }
+
+    desc = strchr(id, ':');
+    if (!desc) {
+        fprintf(stderr, "%s: game id expects a colon in it\n", argv[0]);
+        return 1;
+    }
+    *desc++ = '\0';
+
+    p = default_params();
+    decode_params(p, id);
+    err = validate_desc(p, desc);
+    if (err) {
+        fprintf(stderr, "%s: %s\n", argv[0], err);
+        return 1;
+    }
+    s = new_game(NULL, p, desc);
+
+    solver_diagnostics = diagnostics;
+    sc = solver_make_scratch(p->n);
+    solver_setup_grid(sc, s->numbers->numbers);
+    retd = run_solver(sc, maxdiff);
+    if (retd == 0) {
+        printf("Puzzle is inconsistent\n");
+    } else if (grade) {
+        printf("Difficulty rating: %s\n",
+               dominosa_diffnames[sc->max_diff_used]);
+    } else {
+        char *move, *text;
+        move = solution_move_string(sc);
+        s2 = execute_move(s, move);
+        text = game_text_format(s2);
+        sfree(move);
+        fputs(text, stdout);
+        sfree(text);
+        free_game(s2);
+        if (retd > 1)
+            printf("Could not deduce a unique solution\n");
+    }
+    solver_free_scratch(sc);
+    free_game(s);
+    free_params(p);
+
+    return 0;
+}
+
+#endif
+
 /* vim: set shiftwidth=4 :set textwidth=80: */
 
diff --git a/apps/plugins/puzzles/src/emccpre.js b/apps/plugins/puzzles/src/emccpre.js
index 5082555617..56f69721f7 100644
--- a/apps/plugins/puzzles/src/emccpre.js
+++ b/apps/plugins/puzzles/src/emccpre.js
@@ -212,28 +212,51 @@ function initPuzzle() {
     // button down (our puzzles don't want those events).
     mousedown = Module.cwrap('mousedown', 'void',
                              ['number', 'number', 'number']);
-    buttons_down = 0;
+
+    button_phys2log = [null, null, null];
+    buttons_down = function() {
+        var i, toret = 0;
+        for (i = 0; i < 3; i++)
+            if (button_phys2log[i] !== null)
+                toret |= 1 << button_phys2log[i];
+        return toret;
+    };
+
     onscreen_canvas.onmousedown = function(event) {
+        if (event.button >= 3)
+            return;
+
         var xy = relative_mouse_coords(event, onscreen_canvas);
-        mousedown(xy.x, xy.y, event.button);
-        buttons_down |= 1 << event.button;
+        var logbutton = event.button;
+        if (event.shiftKey)
+            logbutton = 1;   // Shift-click overrides to middle button
+        else if (event.ctrlKey)
+            logbutton = 2;   // Ctrl-click overrides to right button
+
+        mousedown(xy.x, xy.y, logbutton);
+        button_phys2log[event.button] = logbutton;
+
         onscreen_canvas.setCapture(true);
     };
     mousemove = Module.cwrap('mousemove', 'void',
                              ['number', 'number', 'number']);
     onscreen_canvas.onmousemove = function(event) {
-        if (buttons_down) {
+        var down = buttons_down();
+        if (down) {
             var xy = relative_mouse_coords(event, onscreen_canvas);
-            mousemove(xy.x, xy.y, buttons_down);
+            mousemove(xy.x, xy.y, down);
         }
     };
     mouseup = Module.cwrap('mouseup', 'void',
                            ['number', 'number', 'number']);
     onscreen_canvas.onmouseup = function(event) {
-        if (buttons_down & (1 << event.button)) {
-            buttons_down ^= 1 << event.button;
+        if (event.button >= 3)
+            return;
+
+        if (button_phys2log[event.button] !== null) {
             var xy = relative_mouse_coords(event, onscreen_canvas);
-            mouseup(xy.x, xy.y, event.button);
+            mouseup(xy.x, xy.y, button_phys2log[event.button]);
+            button_phys2log[event.button] = null;
         }
     };
 
diff --git a/apps/plugins/puzzles/src/findloop.c b/apps/plugins/puzzles/src/findloop.c
index ffda12d716..4ebdea1f85 100644
--- a/apps/plugins/puzzles/src/findloop.c
+++ b/apps/plugins/puzzles/src/findloop.c
@@ -14,7 +14,7 @@
 #include "puzzles.h"
 
 struct findloopstate {
-    int parent, child, sibling;
+    int parent, child, sibling, component_root;
     bool visited;
     int index, minindex, maxindex;
     int minreachable, maxreachable;
@@ -57,6 +57,33 @@ bool findloop_is_loop_edge(struct findloopstate *pv, int u, int v)
     return !(pv[u].bridge == v || pv[v].bridge == u);
 }
 
+static bool findloop_is_bridge_oneway(
+    struct findloopstate *pv, int u, int v, int *u_vertices, int *v_vertices)
+{
+    int r, total, below;
+
+    if (pv[u].bridge != v)
+        return false;
+
+    r = pv[u].component_root;
+    total = pv[r].maxindex - pv[r].minindex + 1;
+    below = pv[u].maxindex - pv[u].minindex + 1;
+
+    if (u_vertices)
+        *u_vertices = below;
+    if (v_vertices)
+        *v_vertices = total - below;
+
+    return true;
+}
+
+bool findloop_is_bridge(
+    struct findloopstate *pv, int u, int v, int *u_vertices, int *v_vertices)
+{
+    return (findloop_is_bridge_oneway(pv, u, v, u_vertices, v_vertices) ||
+            findloop_is_bridge_oneway(pv, v, u, v_vertices, u_vertices));
+}
+
 bool findloop_run(struct findloopstate *pv, int nvertices,
                   neighbour_fn_t neighbour, void *ctx)
 {
@@ -94,6 +121,7 @@ bool findloop_run(struct findloopstate *pv, int nvertices,
              */
             pv[v].sibling = pv[root].child;
             pv[root].child = v;
+            pv[v].component_root = v;
             debug(("%d is new child of root\n", v));
 
             u = v;
@@ -116,6 +144,7 @@ bool findloop_run(struct findloopstate *pv, int nvertices,
                             pv[w].child = -1;
                             pv[w].sibling = pv[u].child;
                             pv[w].parent = u;
+                            pv[w].component_root = pv[u].component_root;
                             pv[u].child = w;
                         }
 
diff --git a/apps/plugins/puzzles/src/galaxies.c b/apps/plugins/puzzles/src/galaxies.c
index 0cc3198ae0..fe7cd24ecf 100644
--- a/apps/plugins/puzzles/src/galaxies.c
+++ b/apps/plugins/puzzles/src/galaxies.c
@@ -346,29 +346,44 @@ static void add_assoc(const game_state *state, space *tile, space *dot) {
            tile->x, tile->y, dot->x, dot->y, dot->nassoc));*/
 }
 
-static void add_assoc_with_opposite(game_state *state, space *tile, space *dot) {
+static bool ok_to_add_assoc_with_opposite_internal(
+    const game_state *state, space *tile, space *opposite)
+{
     int *colors;
-    space *opposite = space_opposite_dot(state, tile, dot);
+    bool toret;
 
-    if (opposite == NULL) {
-        return;
-    }
-    if (opposite->flags & F_DOT) {
-        return;
-    }
+    if (tile->flags & F_DOT)
+        return false;
+    if (opposite == NULL)
+        return false;
+    if (opposite->flags & F_DOT)
+        return false;
 
+    toret = true;
     colors = snewn(state->w * state->h, int);
     check_complete(state, NULL, colors);
-    if (colors[(tile->y - 1)/2 * state->w + (tile->x - 1)/2]) {
-        sfree(colors);
-        return;
-    }
-    if (colors[(opposite->y - 1)/2 * state->w + (opposite->x - 1)/2]) {
-        sfree(colors);
-        return;
-    }
+
+    if (colors[(tile->y - 1)/2 * state->w + (tile->x - 1)/2])
+        toret = false;
+    if (colors[(opposite->y - 1)/2 * state->w + (opposite->x - 1)/2])
+        toret = false;
 
     sfree(colors);
+    return toret;
+}
+
+static bool ok_to_add_assoc_with_opposite(
+    const game_state *state, space *tile, space *dot)
+{
+    space *opposite = space_opposite_dot(state, tile, dot);
+    return ok_to_add_assoc_with_opposite_internal(state, tile, opposite);
+}
+
+static void add_assoc_with_opposite(game_state *state, space *tile, space *dot) {
+    space *opposite = space_opposite_dot(state, tile, dot);
+
+    assert(ok_to_add_assoc_with_opposite_internal(state, tile, opposite));
+
     remove_assoc_with_opposite(state, tile);
     add_assoc(state, tile, dot);
     remove_assoc_with_opposite(state, opposite);
@@ -2596,8 +2611,15 @@ static char *interpret_move(const game_state *state, game_ui *ui,
          */
         if (INUI(state, px, py)) {
             sp = &SPACE(state, px, py);
+            dot = &SPACE(state, ui->dotx, ui->doty);
 
-            if (!(sp->flags & F_DOT))
+            /*
+             * Exception: if it's not actually legal to add an arrow
+             * and its opposite at this position, we don't try,
+             * because otherwise we'd append an empty entry to the
+             * undo chain.
+             */
+            if (ok_to_add_assoc_with_opposite(state, sp, dot))
                 sprintf(buf + strlen(buf), "%sA%d,%d,%d,%d",
                         sep, px, py, ui->dotx, ui->doty);
         }
diff --git a/apps/plugins/puzzles/src/midend.c b/apps/plugins/puzzles/src/midend.c
index 036c8569c7..a8dd179690 100644
--- a/apps/plugins/puzzles/src/midend.c
+++ b/apps/plugins/puzzles/src/midend.c
@@ -171,6 +171,8 @@ midend *midend_new(frontend *fe, const game *ourgame,
     me->params = ourgame->default_params();
     me->game_id_change_notify_function = NULL;
     me->game_id_change_notify_ctx = NULL;
+    me->encoded_presets = NULL;
+    me->n_encoded_presets = 0;
 
     /*
      * Allow environment-based changing of the default settings by
@@ -261,8 +263,13 @@ static void midend_free_preset_menu(midend *me, struct preset_menu *menu)
 
 void midend_free(midend *me)
 {
+    int i;
+
     midend_free_game(me);
 
+    for (i = 0; i < me->n_encoded_presets; i++)
+        sfree(me->encoded_presets[i]);
+    sfree(me->encoded_presets);
     if (me->drawing)
         drawing_free(me->drawing);
     random_free(me->random);
diff --git a/apps/plugins/puzzles/src/pegs.c b/apps/plugins/puzzles/src/pegs.c
index 32673d56e7..db9caf298f 100644
--- a/apps/plugins/puzzles/src/pegs.c
+++ b/apps/plugins/puzzles/src/pegs.c
@@ -792,6 +792,12 @@ static void game_changed_state(game_ui *ui, const game_state *oldstate,
      * unoccupied.
      */
     ui->dragging = false;
+
+    /*
+     * Also, cancel a keyboard-driven jump if one is half way to being
+     * input.
+     */
+    ui->cur_jumping = false;
 }
 
 #define PREFERRED_TILE_SIZE 33
diff --git a/apps/plugins/puzzles/src/puzzles.h b/apps/plugins/puzzles/src/puzzles.h
index 48d3d83b6e..1732abe3e9 100644
--- a/apps/plugins/puzzles/src/puzzles.h
+++ b/apps/plugins/puzzles/src/puzzles.h
@@ -595,6 +595,32 @@ bool findloop_run(struct findloopstate *state, int nvertices,
 bool findloop_is_loop_edge(struct findloopstate *state, int u, int v);
 
 /*
+ * Alternative query function, which returns true if the u-v edge is a
+ * _bridge_, i.e. a non-loop edge, i.e. an edge whose removal would
+ * disconnect a currently connected component of the graph.
+ *
+ * If the return value is true, then the numbers of vertices that
+ * would be in the new components containing u and v are written into
+ * u_vertices and v_vertices respectively.
+ */
+bool findloop_is_bridge(
+    struct findloopstate *pv, int u, int v, int *u_vertices, int *v_vertices);
+
+/*
+ * Helper function to sort an array. Differs from standard qsort in
+ * that it takes a context parameter that is passed to the compare
+ * function.
+ *
+ * I wrap it in a macro so that you only need to give the element
+ * count of the array. The element size is determined by sizeof.
+ */
+typedef int (*arraysort_cmpfn_t)(const void *av, const void *bv, void *ctx);
+void arraysort_fn(void *array, size_t nmemb, size_t size,
+                  arraysort_cmpfn_t cmp, void *ctx);
+#define arraysort(array, nmemb, cmp, ctx) \
+    arraysort_fn(array, nmemb, sizeof(*(array)), cmp, ctx)
+
+/*
  * Data structure containing the function calls and data specific
  * to a particular game. This is enclosed in a data structure so
  * that a particular platform can choose, if it wishes, to compile
diff --git a/apps/plugins/puzzles/src/sort.c b/apps/plugins/puzzles/src/sort.c
new file mode 100644
index 0000000000..d1897b6fdf
--- /dev/null
+++ b/apps/plugins/puzzles/src/sort.c
@@ -0,0 +1,160 @@
+/*
+ * Implement arraysort() defined in puzzles.h.
+ *
+ * Strategy: heapsort.
+ */
+
+#include <stddef.h>
+#include <string.h>
+
+#include "puzzles.h"
+
+static void memswap(void *av, void *bv, size_t size)
+{
+    char t[4096];
+    char *a = (char *)av, *b = (char *)bv;
+
+    while (size > 0) {
+        size_t thissize = size < sizeof(t) ? size : sizeof(t);
+
+        memcpy(t, a, thissize);
+        memcpy(a, b, thissize);
+        memcpy(b, t, thissize);
+
+        size -= thissize;
+        a += thissize;
+        b += thissize;
+    }
+}
+
+#define PTR(i) ((char *)array + size * (i))
+#define SWAP(i,j) memswap(PTR(i), PTR(j), size)
+#define CMP(i,j) cmp(PTR(i), PTR(j), ctx)
+
+#define LCHILD(i) (2*(i)+1)
+#define RCHILD(i) (2*(i)+2)
+#define PARENT(i) (((i)-1)/2)
+
+static void downheap(void *array, size_t nmemb, size_t size,
+                     arraysort_cmpfn_t cmp, void *ctx, size_t i)
+{
+    while (LCHILD(i) < nmemb) {
+        /* Identify the smallest element out of i and its children. */
+        size_t j = i;
+        if (CMP(j, LCHILD(i)) < 0)
+            j = LCHILD(i);
+        if (RCHILD(i) < nmemb &&
+            CMP(j, RCHILD(i)) < 0)
+            j = RCHILD(i);
+
+        if (j == i)
+            return; /* smallest element is already where it should be */
+
+        SWAP(j, i);
+        i = j;
+    }
+}
+
+void arraysort_fn(void *array, size_t nmemb, size_t size,
+                  arraysort_cmpfn_t cmp, void *ctx)
+{
+    size_t i;
+
+    if (nmemb < 2)
+        return;                        /* trivial */
+
+    /*
+     * Stage 1: build the heap.
+     *
+     * Linear-time if we do it by downheaping the elements in
+     * decreasing order of index, instead of the more obvious approach
+     * of upheaping in increasing order. (Also, it means we don't need
+     * the upheap function at all.)
+     *
+     * We don't need to downheap anything in the second half of the
+     * array, because it can't have any children to swap with anyway.
+     */
+    for (i = PARENT(nmemb-1) + 1; i-- > 0 ;)
+        downheap(array, nmemb, size, cmp, ctx, i);
+
+    /*
+     * Stage 2: dismantle the heap by repeatedly swapping the root
+     * element (at index 0) into the last position and then
+     * downheaping the new root.
+     */
+    for (i = nmemb-1; i > 0; i--) {
+        SWAP(0, i);
+        downheap(array, i, size, cmp, ctx, 0);
+    }
+}
+
+#ifdef SORT_TEST
+
+#include <stdlib.h>
+#include <time.h>
+
+int testcmp(const void *av, const void *bv, void *ctx)
+{
+    int a = *(const int *)av, b = *(const int *)bv;
+    const int *keys = (const int *)ctx;
+    return keys[a] < keys[b] ? -1 : keys[a] > keys[b] ? +1 : 0;
+}
+
+int resetcmp(const void *av, const void *bv)
+{
+    int a = *(const int *)av, b = *(const int *)bv;
+    return a < b ? -1 : a > b ? +1 : 0;
+}
+
+int main(int argc, char **argv)
+{
+    typedef int Array[3723];
+    Array data, keys;
+    int iteration;
+    unsigned seed;
+
+    seed = (argc > 1 ? strtoul(argv[1], NULL, 0) : time(NULL));
+    printf("Random seed = %u\n", seed);
+    srand(seed);
+
+    for (iteration = 0; iteration < 10000; iteration++) {
+        int j;
+        const char *fail = NULL;
+
+        for (j = 0; j < lenof(data); j++) {
+            data[j] = j;
+            keys[j] = rand();
+        }
+
+        arraysort(data, lenof(data), testcmp, keys);
+
+        for (j = 1; j < lenof(data); j++) {
+            if (keys[data[j]] < keys[data[j-1]])
+                fail = "output misordered";
+        }
+        if (!fail) {
+            Array reset;
+            memcpy(reset, data, sizeof(data));
+            qsort(reset, lenof(reset), sizeof(*reset), resetcmp);
+            for (j = 0; j < lenof(reset); j++)
+                if (reset[j] != j)
+                    fail = "output not permuted";
+        }
+
+        if (fail) {
+            printf("Failed at iteration %d: %s\n", iteration, fail);
+            printf("Key values:\n");
+            for (j = 0; j < lenof(keys); j++)
+                printf("  [%2d] %10d\n", j, keys[j]);
+            printf("Output sorted order:\n");
+            for (j = 0; j < lenof(data); j++)
+                printf("  [%2d] %10d\n", data[j], keys[data[j]]);
+            return 1;
+        }
+    }
+
+    printf("OK\n");
+    return 0;
+}
+
+#endif /* SORT_TEST */