Hats tiling: more uniform parent selection.

This tweak improves the uniformity of the generated patches of hat
tiling, by selecting from (the closest 32-bit approximation I can get
to) the limiting probability distribution of finite patches in the
whole plane.

This shouldn't invalidate any grid description that contains enough
coordinates to uniquely specify a piece of tiling - in particular, any
generated by the game itself. But if anyone's been brave enough to
hand-type a grid description in the last two days and left off some of
the coordinates, then those might be invalidated.

auxiliary/hatgen.c

@@ -1476,10 +1476,6 @@ int main(int argc, char **argv)
         printf(" };\n\n");
 
         {
-            struct Parent {
-                MetatileType t;
-                unsigned index;
-            } parents[4][4*MT_MAXEXPAND];
             size_t psizes[4] = {0, 0, 0, 0};
             size_t csizes[4] = {0, 0, 0, 0};
 
@@ -1492,8 +1488,6 @@ int main(int argc, char **argv)
                        "   ", HTPF[i]);
                 for (j = 0; j < nt; j++) {
                     MetatileType c = t[j].type;
-                    parents[c][psizes[c]].t = i;
-                    parents[c][psizes[c]].index = j;
                     psizes[c]++;
                     csizes[i]++;
                     printf(" TT_%c,", HTPF[c]);
@@ -1509,26 +1503,6 @@ int main(int argc, char **argv)
             for (i = 0; i < 4; i++)
                 printf("    %u,\n", (unsigned)csizes[i]);
             printf("};\n\n");
-
-            for (i = 0; i < 4; i++) {
-                printf("static const MetatilePossibleParent "
-                       "permitted_parents_%c[] = {\n", HTPF[i]);
-                for (j = 0; j < psizes[i]; j++)
-                    printf("    { TT_%c, %u },\n", HTPF[parents[i][j].t],
-                           parents[i][j].index);
-                printf("};\n");
-            }
-
-            printf("static const MetatilePossibleParent *const "
-                   "permitted_parents[] = {\n");
-            for (i = 0; i < 4; i++)
-                printf("    permitted_parents_%c,\n", HTPF[i]);
-            printf("};\n");
-
-            printf("static const size_t n_permitted_parents[] = {\n");
-            for (i = 0; i < 4; i++)
-                printf("    %u,\n", (unsigned)psizes[i]);
-            printf("};\n\n");
         }
 
         {

hat-tables.h

@@ -31,69 +31,6 @@ static const size_t nchildren[] = {
     11,
 };
 
-static const MetatilePossibleParent permitted_parents_H[] = {
-    { TT_H, 0 },
-    { TT_H, 1 },
-    { TT_H, 2 },
-    { TT_T, 0 },
-    { TT_P, 0 },
-    { TT_P, 1 },
-    { TT_F, 0 },
-    { TT_F, 1 },
-};
-static const MetatilePossibleParent permitted_parents_T[] = {
-    { TT_H, 3 },
-};
-static const MetatilePossibleParent permitted_parents_P[] = {
-    { TT_H, 4 },
-    { TT_H, 5 },
-    { TT_H, 6 },
-    { TT_T, 1 },
-    { TT_T, 2 },
-    { TT_T, 3 },
-    { TT_P, 2 },
-    { TT_P, 3 },
-    { TT_P, 4 },
-    { TT_F, 2 },
-    { TT_F, 3 },
-};
-static const MetatilePossibleParent permitted_parents_F[] = {
-    { TT_H, 7 },
-    { TT_H, 8 },
-    { TT_H, 9 },
-    { TT_H, 10 },
-    { TT_H, 11 },
-    { TT_H, 12 },
-    { TT_T, 4 },
-    { TT_T, 5 },
-    { TT_T, 6 },
-    { TT_P, 5 },
-    { TT_P, 6 },
-    { TT_P, 7 },
-    { TT_P, 8 },
-    { TT_P, 9 },
-    { TT_P, 10 },
-    { TT_F, 4 },
-    { TT_F, 5 },
-    { TT_F, 6 },
-    { TT_F, 7 },
-    { TT_F, 8 },
-    { TT_F, 9 },
-    { TT_F, 10 },
-};
-static const MetatilePossibleParent *const permitted_parents[] = {
-    permitted_parents_H,
-    permitted_parents_T,
-    permitted_parents_P,
-    permitted_parents_F,
-};
-static const size_t n_permitted_parents[] = {
-    8,
-    1,
-    11,
-    22,
-};
-
 static const KitemapEntry kitemap_H[] = {
     /* hat #0 in metatile #0 (type H) */
     {1,0,0}, {7,3,0}, {3,0,4}, {4,0,4},

hat.c

@@ -318,11 +318,6 @@ static const char tilechars[] = "HTPF";
  * Definitions for the autogenerated hat-tables.h header file that
  * defines all the lookup tables.
  */
-typedef struct MetatilePossibleParent {
-    TileType type;
-    unsigned index;
-} MetatilePossibleParent;
-
 typedef struct KitemapEntry {
     int kite, hat, meta;               /* all -1 if impossible */
 } KitemapEntry;
@@ -505,15 +500,189 @@ static HatCoords *initial_coords(HatCoordContext *ctx)
  */
 static void ensure_coords(HatCoordContext *ctx, HatCoords *hc, size_t n)
 {
+    /*
+     * One table that we write by hand: the permitted ways to extend
+     * the coordinate system outwards from a given metatile.
+     *
+     * One obvious approach would be to make a table of all the places
+     * each metatile can appear in the expansion of another (e.g. H
+     * can be subtile 0, 1 or 2 of another H, subtile 0 of a T, or 0
+     * or 1 of a P or an F), and when we need to decide what our
+     * current topmost tile turns out to be a subtile of, choose
+     * equiprobably at random from those options.
+     *
+     * That's what I did originally, but a better approach is to skew
+     * the probabilities. We'd like to generate our patch of actual
+     * tiling uniformly at random, in the sense that if you selected
+     * uniformly from a very large region of the plane, the
+     * distribution of possible finite patches of tiling would
+     * converge to some limit as that region tended to infinity, and
+     * we'd be picking from that limiting distribution on finite
+     * patches.
+     *
+     * For this we have to refer back to the original paper, which
+     * indicates the subset of each metatile's expansion that can be
+     * considered to 'belong' to that metatile, such that every
+     * subtile belongs to exactly one parent metatile, and the
+     * overlaps are eliminated. Reading out the diagrams from their
+     * Figure 2.8:
+     *
+     * - H: we discard three of the outer F subtiles, in the symmetric
+     *    positions index by our coordinates as 7, 10, 11. So we keep
+     *    the remaining subtiles {0,1,2,3,4,5,6,8,9,12}, which consist
+     *    of three H, one T, three P and three F.
+     *
+     * - T: only the central H expanded from a T is considered to
+     *   belong to it, so we just keep {0}, a single H.
+     *
+     * - P: we discard everything intersected by a long edge of the
+     *   parallelogram, leaving the central three tiles and the
+     *   endmost pair of F. That is, we keep {0,1,4,5,10}, consisting
+     *   of two H, one P and two F.
+     *
+     * - F: looks like P at one end, and we retain the corresponding
+     *   set of tiles there, but at the other end we keep the two F on
+     *   either side of the endmost one. So we keep {0,1,3,6,8,10},
+     *   consisting of two H, one P and _three_ F.
+     *
+     * Adding up the tile numbers gives us this matrix system:
+     *
+     * (H_1)   (3 1 2 2)(H_0)
+     * (T_1) = (1 0 0 0)(T_0)
+     * (P_1)   (3 0 1 1)(P_0)
+     * (F_1)   (3 0 2 3)(F_0)
+     *
+     * which says that if you have a patch of metatiling consisting of
+     * H_0 H tiles, T_0 T tiles etc, then this matrix shows the number
+     * H_1 of smaller H tiles, etc, expanded from it.
+     *
+     * If you expand _many_ times, that's equivalent to raising the
+     * matrix to a power:
+     *
+     *                  n
+     * (H_n)   (3 1 2 2) (H_0)
+     * (T_n) = (1 0 0 0) (T_0)
+     * (P_n)   (3 0 1 1) (P_0)
+     * (F_n)   (3 0 2 3) (F_0)
+     *
+     * The limiting distribution of metatiles is obtained by looking
+     * at the four-way ratio between H_n, T_n, P_n and F_n as n tends
+     * to infinity. To calculate this, we find the eigenvalues and
+     * eigenvectors of the matrix, and extract the eigenvector
+     * corresponding to the eigenvalue of largest magnitude. (Things
+     * get more complicated in cases where that's not unique, but
+     * here, it is.)
+     *
+     * That eigenvector is
+     *
+     *   [          1          ]      [ 1                      ]
+     *   [ (7 - 3 sqrt(5)) / 2 ]  ~=  [ 0.14589803375031545538 ]
+     *   [    3 sqrt(5) - 6    ]      [ 0.70820393249936908922 ]
+     *   [ (9 - 3 sqrt(5)) / 2 ]      [ 1.14589803375031545538 ]
+     *
+     * So those are the limiting relative proportions of metatiles.
+     *
+     * So if we have a particular metatile, how likely is it for its
+     * parent to be one of those? We have to adjust by the number of
+     * metatiles of each type that each tile has as its children. For
+     * example, the P and F tiles have one P child each, but the H has
+     * three P children. So if we have a P, the proportion of H in its
+     * potential ancestry is three times what's shown here. (And T
+     * can't occur at all as a parent.)
+     *
+     * In other words, we should choose _each coordinate_ with
+     * probability corresponding to one of those numbers (scaled down
+     * so they all sum to 1). Continuing to use P as an example, it
+     * will be:
+     *
+     *  - child 4 of H with relative probability 1
+     *  - child 5 of H with relative probability 1
+     *  - child 6 of H with relative probability 1
+     *  - child 4 of P with relative probability 0.70820393249936908922
+     *  - child 3 of F with relative probability 1.14589803375031545538
+     *
+     * and then we obtain the true probabilities by scaling those
+     * values down so that they sum to 1.
+     *
+     * The tables below give a reasonable approximation in 32-bit
+     * integers to these proportions.
+     */
+
+    typedef struct MetatilePossibleParent {
+        TileType type;
+        unsigned index;
+        unsigned long probability;
+    } MetatilePossibleParent;
+
+    /* The above probabilities scaled up by 10000000 */
+    #define PROB_H 10000000
+    #define PROB_T  1458980
+    #define PROB_P  7082039
+    #define PROB_F 11458980
+
+    static const MetatilePossibleParent parents_H[] = {
+        { TT_H,  0, PROB_H },
+        { TT_H,  1, PROB_H },
+        { TT_H,  2, PROB_H },
+        { TT_T,  0, PROB_T },
+        { TT_P,  0, PROB_P },
+        { TT_P,  1, PROB_P },
+        { TT_F,  0, PROB_F },
+        { TT_F,  1, PROB_F },
+    };
+    static const MetatilePossibleParent parents_T[] = {
+        { TT_H,  3, PROB_H },
+    };
+    static const MetatilePossibleParent parents_P[] = {
+        { TT_H,  4, PROB_H },
+        { TT_H,  5, PROB_H },
+        { TT_H,  6, PROB_H },
+        { TT_P,  4, PROB_P },
+        { TT_F,  3, PROB_F },
+    };
+    static const MetatilePossibleParent parents_F[] = {
+        { TT_H,  8, PROB_H },
+        { TT_H,  9, PROB_H },
+        { TT_H, 12, PROB_H },
+        { TT_P,  5, PROB_P },
+        { TT_P, 10, PROB_P },
+        { TT_F,  6, PROB_F },
+        { TT_F,  8, PROB_F },
+        { TT_F, 10, PROB_F },
+    };
+
+    #undef PROB_H
+    #undef PROB_T
+    #undef PROB_P
+    #undef PROB_F
+
+    static const MetatilePossibleParent *const possible_parents[] = {
+        parents_H, parents_T, parents_P, parents_F,
+    };
+    static const size_t n_possible_parents[] = {
+        lenof(parents_H), lenof(parents_T), lenof(parents_P), lenof(parents_F),
+    };
+
     if (ctx->prototype->nc < n) {
         hc_make_space(ctx->prototype, n);
         while (ctx->prototype->nc < n) {
             TileType type = ctx->prototype->c[ctx->prototype->nc - 1].type;
             assert(ctx->prototype->c[ctx->prototype->nc - 1].index == -1);
-            const MetatilePossibleParent *parents = permitted_parents[type];
+            const MetatilePossibleParent *parents = possible_parents[type];
             size_t parent_index;
             if (ctx->rs) {
-                parent_index = random_upto(ctx->rs, n_permitted_parents[type]);
+                unsigned long limit = 0, value;
+                size_t nparents = n_possible_parents[type], i;
+                for (i = 0; i < nparents; i++)
+                    limit += parents[i].probability;
+                value = random_upto(ctx->rs, limit);
+                for (i = 0; i < nparents; i++) {
+                    if (value < parents[i].probability)
+                        break;
+                    value -= parents[i].probability;
+                }
+                assert(i < nparents);
+                parent_index = i;
             } else {
                 parent_index = 0;
             }