shithub: puzzles

ref: 1af1204b9c33c9c03b2e0fe66c2f07d9729cbc72
dir: /hat.c/

View raw version
/*
 * Code to generate patches of the aperiodic 'hat' tiling discovered
 * in 2023.
 *
 * auxiliary/doc/hats.html contains an explanation of the basic ideas
 * of this algorithm, which can't really be put in a source file
 * because it just has too many complicated diagrams. So read that
 * first, because the comments in here will refer to it.
 *
 * Discoverers' website: https://cs.uwaterloo.ca/~csk/hat/
 * Preprint of paper:    https://arxiv.org/abs/2303.10798
 */

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "puzzles.h"
#include "hat.h"

/*
 * Coordinate system:
 *
 * The output of this code lives on the tiling known to grid.c as
 * 'Kites', which can be viewed as a tiling of hexagons each of which
 * is subdivided into six kites sharing their pointy vertex, or
 * (equivalently) a tiling of equilateral triangles each subdivided
 * into three kits sharing their blunt vertex.
 *
 * We express coordinates in this system relative to the basis (1, r)
 * where r = (1 + sqrt(3)i) / 2 is a primitive 6th root of unity. This
 * gives us a system in which two integer coordinates can address any
 * grid point, provided we scale up so that the side length of the
 * equilateral triangles in the tiling is 6.
 */

typedef struct Point {
    int x, y;                          /* represents x + yr */
} Point;

static inline Point pointscale(int scale, Point a)
{
    Point r = { scale * a.x, scale * a.y };
    return r;
}

static inline Point pointadd(Point a, Point b)
{
    Point r = { a.x + b.x, a.y + b.y };
    return r;
}

/*
 * We identify a single kite by the coordinates of its four vertices.
 * This allows us to construct the coordinates of an adjacent kite by
 * taking affine transformations of the original kite's vertices.
 *
 * This is a useful way to do it because it means that if you reflect
 * the kite (by swapping its left and right vertices) then these
 * transformations also perform in a reflected way. This will be
 * useful in the code below that outputs the coordinates of each hat,
 * because this way it can work by walking around its 8 kites using a
 * fixed set of steps, and if the hat is reflected, then we just
 * reflect the starting kite before doing that, and everything still
 * works.
 */

typedef struct Kite {
    Point centre, left, right, outer;
} Kite;

static inline Kite kite_left(Kite k)
{
    Kite r;
    r.centre = k.centre;
    r.right = k.left;
    r.outer = pointadd(pointscale(2, k.left), pointscale(-1, k.outer));
    r.left = pointadd(pointadd(k.centre, k.left), pointscale(-1, k.right));
    return r;
}

static inline Kite kite_right(Kite k)
{
    Kite r;
    r.centre = k.centre;
    r.left = k.right;
    r.outer = pointadd(pointscale(2, k.right), pointscale(-1, k.outer));
    r.right = pointadd(pointadd(k.centre, k.right), pointscale(-1, k.left));
    return r;
}

static inline Kite kite_forward_left(Kite k)
{
    Kite r;
    r.outer = k.outer;
    r.right = k.left;
    r.centre = pointadd(pointscale(2, k.left), pointscale(-1, k.centre));
    r.left = pointadd(pointadd(k.right, k.left), pointscale(-1, k.centre));
    return r;
}

static inline Kite kite_forward_right(Kite k)
{
    Kite r;
    r.outer = k.outer;
    r.left = k.right;
    r.centre = pointadd(pointscale(2, k.right), pointscale(-1, k.centre));
    r.right = pointadd(pointadd(k.left, k.right), pointscale(-1, k.centre));
    return r;
}

typedef enum KiteStep { KS_LEFT, KS_RIGHT, KS_F_LEFT, KS_F_RIGHT } KiteStep;

static inline Kite kite_step(Kite k, KiteStep step)
{
    switch (step) {
      case KS_LEFT: return kite_left(k);
      case KS_RIGHT: return kite_right(k);
      case KS_F_LEFT: return kite_forward_left(k);
      default /* case KS_F_RIGHT */: return kite_forward_right(k);
    }
}

/*
 * Function to enumerate the kites in a rectangular region, in a
 * serpentine-raster fashion so that every kite delivered shares an
 * edge with a recent previous one.
 */
#define KE_NKEEP 3
typedef struct KiteEnum {
    /* Fields private to the enumerator */
    int state;
    int x, y, w, h;
    unsigned curr_index;

    /* Fields the client can legitimately read out */
    Kite *curr;
    Kite recent[KE_NKEEP];
    unsigned last_index;
    KiteStep last_step; /* step that got curr from recent[last_index] */
} KiteEnum;

static void first_kite(KiteEnum *s, int w, int h)
{
    Kite start = { {0,0}, {0, 3}, {3, 0}, {2, 2} };
    size_t i;

    for (i = 0; i < KE_NKEEP; i++)
        s->recent[i] = start;          /* initialise to *something* */
    s->curr_index = 0;
    s->curr = &s->recent[s->curr_index];
    s->state = 1;
    s->w = w;
    s->h = h;
    s->x = 0;
    s->y = 0;
}
static bool next_kite(KiteEnum *s)
{
    unsigned lastbut1 = s->last_index;
    s->last_index = s->curr_index;
    s->curr_index = (s->curr_index + 1) % KE_NKEEP;
    s->curr = &s->recent[s->curr_index];

    switch (s->state) {
        /* States 1,2,3 walk rightwards along the upper side of a
         * horizontal grid line with a pointy kite end at the start
         * point */
      case 1:
        s->last_step = KS_F_RIGHT;
        s->state = 2;
        break;

      case 2:
        if (s->x+1 >= s->w) {
            s->last_step = KS_F_RIGHT;
            s->state = 4;
            break;
        }
        s->last_step = KS_RIGHT;
        s->state = 3;
        s->x++;
        break;

      case 3:
        s->last_step = KS_RIGHT;
        s->state = 1;
        break;

        /* State 4 is special: we've just moved up into a row below a
         * grid line, but we can't produce the rightmost tile of that
         * row because it's not adjacent any tile so far emitted. So
         * instead, emit the second-rightmost tile, and next time,
         * we'll emit the rightmost. */
      case 4:
        s->last_step = KS_LEFT;
        s->state = 5;
        break;

        /* And now we have to emit the third-rightmost tile relative
         * to the last but one tile we emitted (the one from state 2,
         * not state 4). */
      case 5:
        s->last_step = KS_RIGHT;
        s->last_index = lastbut1;
        s->state = 6;
        break;

        /* Now states 6-8 handle the general case of walking leftwards
         * along the lower side of a line, starting from a
         * right-angled kite end. */
      case 6:
        if (s->x <= 0) {
            if (s->y+1 >= s->h) {
                s->state = 0;
                return false;
            }
            s->last_step = KS_RIGHT;
            s->state = 9;
            s->y++;
            break;
        }
        s->last_step = KS_F_RIGHT;
        s->state = 7;
        s->x--;
        break;

      case 7:
        s->last_step = KS_RIGHT;
        s->state = 8;
        break;

      case 8:
        s->last_step = KS_RIGHT;
        s->state = 6;
        break;

        /* States 9,10,11 walk rightwards along the upper side of a
         * horizontal grid line with a right-angled kite end at the
         * start point. This time there's no awkward transition from
         * the previous row. */
      case 9:
        s->last_step = KS_RIGHT;
        s->state = 10;
        break;

      case 10:
        s->last_step = KS_RIGHT;
        s->state = 11;
        break;

      case 11:
        if (s->x+1 >= s->w) {
            /* Another awkward transition to the next row, where we
             * have to generate it based on the previous state-9 tile.
             * But this time at least we generate the rightmost tile
             * of the new row, so the next states will be simple. */
            s->last_step = KS_F_RIGHT;
            s->last_index = lastbut1;
            s->state = 12;
            break;
        }
        s->last_step = KS_F_RIGHT;
        s->state = 9;
        s->x++;
        break;

        /* States 12,13,14 walk leftwards along the upper edge of a
         * horizontal grid line with a pointy kite end at the start
         * point */
      case 12:
        s->last_step = KS_F_RIGHT;
        s->state = 13;
        break;

      case 13:
        if (s->x <= 0) {
            if (s->y+1 >= s->h) {
                s->state = 0;
                return false;
            }
            s->last_step = KS_LEFT;
            s->state = 1;
            s->y++;
            break;
        }
        s->last_step = KS_RIGHT;
        s->state = 14;
        s->x--;
        break;

      case 14:
        s->last_step = KS_RIGHT;
        s->state = 12;
        break;

      default:
        return false;
    }

    *s->curr = kite_step(s->recent[s->last_index], s->last_step);
    return true;
}

/*
 * Assorted useful definitions.
 */
typedef enum TileType { TT_H, TT_T, TT_P, TT_F, TT_KITE, TT_HAT } TileType;
static const char tilechars[] = "HTPF";

#define HAT_KITES 8     /* number of kites in a hat */
#define MT_MAXEXPAND 13 /* largest number of metatiles in any expansion */

/*
 * Definitions for the autogenerated hat-tables.h header file that
 * defines all the lookup tables.
 */
typedef struct KitemapEntry {
    int kite, hat, meta;               /* all -1 if impossible */
} KitemapEntry;

typedef struct MetamapEntry {
    int meta, meta2;
} MetamapEntry;

static inline size_t kitemap_index(KiteStep step, unsigned kite,
                                   unsigned hat, unsigned meta)
{
    return step + 4 * (kite + 8 * (hat + 4 * meta));
}

static inline size_t metamap_index(unsigned meta, unsigned meta2)
{
    return meta2 * MT_MAXEXPAND + meta;
}

/*
 * The actual tables.
 */
#include "hat-tables.h"

/*
 * One set of tables 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 there isn't a _unique_ eigenvalue
 * of largest magnitude, but here, there 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 },
};

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),
};

/*
 * Similarly, we also want to choose our absolute starting hat with
 * close to uniform probability, which again we do by looking at the
 * limiting ratio of the metatile types, and this time, scaling by the
 * number of hats in each metatile.
 *
 * We cheatingly use the same MetatilePossibleParent struct, because
 * it's got all the right fields, even if it has an inappropriate
 * name.
 */
static const MetatilePossibleParent starting_hats[] = {
    { TT_H,  0, PROB_H },
    { TT_H,  1, PROB_H },
    { TT_H,  2, PROB_H },
    { TT_H,  3, PROB_H },
    { TT_T,  0, PROB_P },
    { TT_P,  0, PROB_P },
    { TT_P,  1, PROB_P },
    { TT_F,  0, PROB_F },
    { TT_F,  1, PROB_F },
};

#undef PROB_H
#undef PROB_T
#undef PROB_P
#undef PROB_F

/*
 * Coordinate system for tracking kites within a randomly selected
 * part of the recursively expanded hat tiling.
 *
 * HatCoords will store an array of HatCoord, in little-endian
 * arrangement. So hc->c[0] will always have type TT_KITE and index a
 * single kite within a hat; hc->c[1] will have type TT_HAT and index
 * a hat within a first-order metatile; hc->c[2] will be the smallest
 * metatile containing this hat, and hc->c[3, 4, 5, ...] will be
 * higher-order metatiles as needed.
 *
 * The last coordinate stored, hc->c[hc->nc-1], will have a tile type
 * but no index (represented by index==-1). This means "we haven't
 * decided yet what this level of metatile needs to be". If we need to
 * refer to this level during the step_coords algorithm, we make it up
 * at random, based on a table of what metatiles each type can
 * possibly be part of, at what index.
 */
typedef struct HatCoord {
    int index; /* index within that tile, or -1 if not yet known */
    TileType type;  /* type of this tile */
} HatCoord;

typedef struct HatCoords {
    HatCoord *c;
    size_t nc, csize;
} HatCoords;

static HatCoords *hc_new(void)
{
    HatCoords *hc = snew(HatCoords);
    hc->nc = hc->csize = 0;
    hc->c = NULL;
    return hc;
}

static void hc_free(HatCoords *hc)
{
    if (hc) {
        sfree(hc->c);
        sfree(hc);
    }
}

static void hc_make_space(HatCoords *hc, size_t size)
{
    if (hc->csize < size) {
        hc->csize = hc->csize * 5 / 4 + 16;
        if (hc->csize < size)
            hc->csize = size;
        hc->c = sresize(hc->c, hc->csize, HatCoord);
    }
}

static HatCoords *hc_copy(HatCoords *hc_in)
{
    HatCoords *hc_out = hc_new();
    hc_make_space(hc_out, hc_in->nc);
    memcpy(hc_out->c, hc_in->c, hc_in->nc * sizeof(*hc_out->c));
    hc_out->nc = hc_in->nc;
    return hc_out;
}

static const MetatilePossibleParent *choose_mpp(
    random_state *rs, const MetatilePossibleParent *parents, size_t nparents)
{
    /*
     * If we needed to do this _efficiently_, we'd rewrite all those
     * tables above as cumulative frequency tables and use binary
     * search. But this happens about log n times in a grid of area n,
     * so it hardly matters, and it's easier to keep the tables
     * legible.
     */
    unsigned long limit = 0, value;
    size_t i;

    for (i = 0; i < nparents; i++)
        limit += parents[i].probability;

    value = random_upto(rs, limit);

    for (i = 0; i+1 < nparents; i++) {
        if (value < parents[i].probability)
            return &parents[i];
        value -= parents[i].probability;
    }

    assert(i == nparents - 1);
    assert(value < parents[i].probability);
    return &parents[i];
}

/*
 * HatCoordContext is the shared context of a whole run of the
 * algorithm. Its 'prototype' HatCoords object represents the
 * coordinates of the starting kite, and is extended as necessary; any
 * other HatCoord that needs extending will copy the higher-order
 * values from ctx->prototype as needed, so that once each choice has
 * been made, it remains consistent.
 *
 * When we're inventing a random piece of tiling in the first place,
 * we append to ctx->prototype by choosing a random (but legal)
 * higher-level metatile for the current topmost one to turn out to be
 * part of. When we're replaying a generation whose parameters are
 * already stored, we don't have a random_state, and we make fixed
 * decisions if not enough coordinates were provided.
 *
 * (Of course another approach would be to reject grid descriptions
 * that didn't define enough coordinates! But that would involve a
 * whole extra iteration over the whole grid region just for
 * validation, and that seems like more timewasting than really
 * needed. So we tolerate short descriptions, and do something
 * deterministic with them.)
 */

typedef struct HatCoordContext {
    random_state *rs;
    HatCoords *prototype;
} HatCoordContext;

static void init_coords_random(HatCoordContext *ctx, random_state *rs)
{
    const MetatilePossibleParent *starting_hat = choose_mpp(
        rs, starting_hats, lenof(starting_hats));

    ctx->rs = rs;
    ctx->prototype = hc_new();
    hc_make_space(ctx->prototype, 3);
    ctx->prototype->c[2].type = starting_hat->type;
    ctx->prototype->c[2].index = -1;
    ctx->prototype->c[1].type = TT_HAT;
    ctx->prototype->c[1].index = starting_hat->index;
    ctx->prototype->c[0].type = TT_KITE;
    ctx->prototype->c[0].index = random_upto(rs, HAT_KITES);
    ctx->prototype->nc = 3;
}

static inline int metatile_char_to_enum(char metatile)
{
    return (metatile == 'H' ? TT_H :
            metatile == 'T' ? TT_T :
            metatile == 'P' ? TT_P :
            metatile == 'F' ? TT_F : -1);
}

static void init_coords_params(HatCoordContext *ctx,
                               const struct HatPatchParams *hp)
{
    size_t i;

    ctx->rs = NULL;
    ctx->prototype = hc_new();

    assert(hp->ncoords >= 3);

    hc_make_space(ctx->prototype, hp->ncoords + 1);
    ctx->prototype->nc = hp->ncoords + 1;

    for (i = 0; i < hp->ncoords; i++)
        ctx->prototype->c[i].index = hp->coords[i];

    ctx->prototype->c[hp->ncoords].type =
        metatile_char_to_enum(hp->final_metatile);
    ctx->prototype->c[hp->ncoords].index = -1;

    ctx->prototype->c[0].type = TT_KITE;
    ctx->prototype->c[1].type = TT_HAT;

    for (i = hp->ncoords - 1; i > 1; i--) {
        TileType metatile = ctx->prototype->c[i+1].type;
        assert(hp->coords[i] < nchildren[metatile]);
        ctx->prototype->c[i].type = children[metatile][hp->coords[i]];
    }

    assert(hp->coords[0] < 8);
}

static HatCoords *initial_coords(HatCoordContext *ctx)
{
    return hc_copy(ctx->prototype);
}

/*
 * Extend hc until it has at least n coordinates in, by copying from
 * ctx->prototype if needed, and extending ctx->prototype if needed in
 * order to do that.
 */
static void ensure_coords(HatCoordContext *ctx, HatCoords *hc, size_t n)
{
    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 *parent;

            if (ctx->rs)
                parent = choose_mpp(ctx->rs, possible_parents[type],
                                    n_possible_parents[type]);
            else
                parent = possible_parents[type];

            ctx->prototype->c[ctx->prototype->nc - 1].index = parent->index;
            ctx->prototype->c[ctx->prototype->nc].index = -1;
            ctx->prototype->c[ctx->prototype->nc].type = parent->type;
            ctx->prototype->nc++;
        }
    }

    hc_make_space(hc, n);
    while (hc->nc < n) {
        assert(hc->c[hc->nc - 1].index == -1);
        assert(hc->c[hc->nc - 1].type == ctx->prototype->c[hc->nc - 1].type);
        hc->c[hc->nc - 1].index = ctx->prototype->c[hc->nc - 1].index;
        hc->c[hc->nc].index = -1;
        hc->c[hc->nc].type = ctx->prototype->c[hc->nc].type;
        hc->nc++;
    }
}

static void cleanup_coords(HatCoordContext *ctx)
{
    hc_free(ctx->prototype);
}

#ifdef DEBUG_COORDS
static inline void debug_coords(const char *prefix, HatCoords *hc,
                                const char *suffix)
{
    const char *sep = "";
    static const char *const types[] = {"H","T","P","F","kite","hat"};

    fputs(prefix, stderr);
    for (size_t i = 0; i < hc->nc; i++) {
        fprintf(stderr, "%s %s ", sep, types[hc->c[i].type]);
        sep = " .";
        if (hc->c[i].index == -1)
            fputs("?", stderr);
        else
            fprintf(stderr, "%d", hc->c[i].index);
    }
    fputs(suffix, stderr);
}
#else
#define debug_coords(p,c,s) ((void)0)
#endif

/*
 * The actual system for finding the coordinates of an adjacent kite.
 */

/*
 * Kitemap step: ensure we have enough coordinates to know two levels
 * of meta-tiling, and use the kite map for the outer layer to move
 * around the individual kites. If this fails, return NULL.
 */
static HatCoords *try_step_coords_kitemap(
    HatCoordContext *ctx, HatCoords *hc_in, KiteStep step)
{
    ensure_coords(ctx, hc_in, 4);
    debug_coords("  try kitemap  ", hc_in, "\n");
    unsigned kite = hc_in->c[0].index;
    unsigned hat = hc_in->c[1].index;
    unsigned meta = hc_in->c[2].index;
    TileType meta2type = hc_in->c[3].type;
    const KitemapEntry *ke = &kitemap[meta2type][
        kitemap_index(step, kite, hat, meta)];
    if (ke->kite >= 0) {
        /*
         * Success! We've got coordinates for the next kite in this
         * direction.
         */
        HatCoords *hc_out = hc_copy(hc_in);

        hc_out->c[2].index = ke->meta;
        hc_out->c[2].type = children[meta2type][ke->meta];
        hc_out->c[1].index = ke->hat;
        hc_out->c[1].type = TT_HAT;
        hc_out->c[0].index = ke->kite;
        hc_out->c[0].type = TT_KITE;

        debug_coords("  success!     ", hc_out, "\n");
        return hc_out;
    }

    return NULL;
}

/*
 * Recursive metamap step. Try using the metamap to rewrite the
 * coordinates at hc->c[depth] and hc->c[depth+1] (using the metamap
 * for the tile type described in hc->c[depth+2]). If successful,
 * recurse back down to see if this led to a successful step via the
 * kitemap. If even that fails (so that we need to try a higher-order
 * metamap rewrite), return NULL.
 */
static HatCoords *try_step_coords_metamap(
    HatCoordContext *ctx, HatCoords *hc_in, KiteStep step, size_t depth)
{
    HatCoords *hc_tmp = NULL, *hc_out;

    ensure_coords(ctx, hc_in, depth+3);
#ifdef DEBUG_COORDS
    fprintf(stderr, "  try meta %-4d", (int)depth);
    debug_coords("", hc_in, "\n");
#endif
    unsigned meta_orig = hc_in->c[depth].index;
    unsigned meta2_orig = hc_in->c[depth+1].index;
    TileType meta3type = hc_in->c[depth+2].type;

    unsigned meta = meta_orig, meta2 = meta2_orig;

    while (true) {
        const MetamapEntry *me;
        HatCoords *hc_curr = hc_tmp ? hc_tmp : hc_in;

        if (depth > 2)
            hc_out = try_step_coords_metamap(ctx, hc_curr, step, depth - 1);
        else
            hc_out = try_step_coords_kitemap(ctx, hc_curr, step);
        if (hc_out) {
            hc_free(hc_tmp);
            return hc_out;
        }

        me = &metamap[meta3type][metamap_index(meta, meta2)];
        assert(me->meta != -1);
        if (me->meta == meta_orig && me->meta2 == meta2_orig) {
            hc_free(hc_tmp);
            return NULL;
        }

        meta = me->meta;
        meta2 = me->meta2;

        /*
         * We must do the rewrite in a copy of hc_in. It's not
         * _necessarily_ obvious that that's the case (any successful
         * rewrite leaves the coordinates still valid and still
         * referring to the same kite, right?). But the problem is
         * that we might do a rewrite at this level more than once,
         * and in between, a metamap rewrite at the next level down
         * might have modified _one_ of the two coordinates we're
         * messing about with. So it's easiest to let the recursion
         * just use a separate copy.
         */
        if (!hc_tmp)
            hc_tmp = hc_copy(hc_in);

        hc_tmp->c[depth+1].index = meta2;
        hc_tmp->c[depth+1].type = children[meta3type][meta2];
        hc_tmp->c[depth].index = meta;
        hc_tmp->c[depth].type = children[hc_tmp->c[depth+1].type][meta];

        debug_coords("  rewritten -> ", hc_tmp, "\n");
    }
}

/*
 * The top-level algorithm for finding the next tile.
 */
static HatCoords *step_coords(HatCoordContext *ctx, HatCoords *hc_in,
                              KiteStep step)
{
    HatCoords *hc_out;
    size_t depth;

#ifdef DEBUG_COORDS
    static const char *const directions[] = {
        " left\n", " right\n", " forward left\n", " forward right\n" };
    debug_coords("step start     ", hc_in, directions[step]);
#endif

    /*
     * First, just try a kitemap step immediately. If that succeeds,
     * we're done.
     */
    if ((hc_out = try_step_coords_kitemap(ctx, hc_in, step)) != NULL)
        return hc_out;

    /*
     * Otherwise, try metamap rewrites at successively higher layers
     * until one works. Each one will recurse back down to the
     * kitemap, as described above.
     */
    for (depth = 2;; depth++) {
        if ((hc_out = try_step_coords_metamap(
                 ctx, hc_in, step, depth)) != NULL)
            return hc_out;
    }
}

/*
 * Generate a random set of parameters for a tiling of a given size.
 * To do this, we iterate over the whole tiling via first_kite and
 * next_kite, and for each kite, calculate its coordinates. But then
 * we throw the coordinates away and don't do anything with them!
 *
 * But the side effect of _calculating_ all those coordinates is that
 * we found out how far ctx->prototype needed to be extended, and did
 * so, pulling random choices out of our random_state. So after this
 * iteration, ctx->prototype contains everything we need to replicate
 * the same piece of tiling next time.
 */
void hat_tiling_randomise(struct HatPatchParams *hp, int w, int h,
                          random_state *rs)
{
    HatCoordContext ctx[1];
    HatCoords *coords[KE_NKEEP];
    KiteEnum s[1];
    size_t i;

    init_coords_random(ctx, rs);
    for (i = 0; i < lenof(coords); i++)
        coords[i] = NULL;

    first_kite(s, w, h);
    coords[s->curr_index] = initial_coords(ctx);

    while (next_kite(s)) {
        hc_free(coords[s->curr_index]);
        coords[s->curr_index] = step_coords(
            ctx, coords[s->last_index], s->last_step);
    }

    hp->ncoords = ctx->prototype->nc - 1;
    hp->coords = snewn(hp->ncoords, unsigned char);
    for (i = 0; i < hp->ncoords; i++)
        hp->coords[i] = ctx->prototype->c[i].index;
    hp->final_metatile = tilechars[ctx->prototype->c[hp->ncoords].type];

    cleanup_coords(ctx);
    for (i = 0; i < lenof(coords); i++)
        hc_free(coords[i]);
}

const char *hat_tiling_params_invalid(const struct HatPatchParams *hp)
{
    TileType metatile;
    size_t i;

    if (hp->ncoords < 3)
        return "Grid parameters require at least three coordinates";
    if (metatile_char_to_enum(hp->final_metatile) < 0)
        return "Grid parameters contain an invalid final metatile";
    if (hp->coords[0] >= 8)
        return "Grid parameters contain an invalid kite index";

    metatile = metatile_char_to_enum(hp->final_metatile);
    for (i = hp->ncoords - 1; i > 1; i--) {
        if (hp->coords[i] >= nchildren[metatile])
            return "Grid parameters contain an invalid metatile index";
        metatile = children[metatile][hp->coords[i]];
    }

    if (hp->coords[1] >= hats_in_metatile[metatile])
        return "Grid parameters contain an invalid hat index";

    return NULL;
}

/*
 * For each kite generated by hat_tiling_generate, potentially
 * generate an output hat and give it to our caller.
 *
 * We do this by starting from kite #0 of each hat, and tracing round
 * the boundary. If the whole boundary is within the caller's bounding
 * region, we return it; if it goes off the edge, we don't.
 *
 * (Of course, every hat we _do_ want to return will have all its
 * kites inside the rectangle, so its kite #0 will certainly be caught
 * by this iteration.)
 */

typedef void (*internal_hat_callback_fn)(void *ctx, Kite kite0, HatCoords *hc,
                                         int *coords);

static void maybe_report_hat(int w, int h, Kite kite, HatCoords *hc,
                             internal_hat_callback_fn cb, void *cbctx)
{
    Kite kite0;
    Point vertices[14];
    size_t i, j;
    bool reversed = false;
    int coords[28];

    /* Only iterate from kite #0 of a hat */
    if (hc->c[0].index != 0)
        return;
    kite0 = kite;

    /*
     * Identify reflected hats: they are always hat #3 of an H
     * metatile. If we find one, reflect the starting kite so that the
     * kite_step operations below will go in the other direction.
     */
    if (hc->c[2].type == TT_H && hc->c[1].index == 3) {
        reversed = true;
        Point tmp = kite.left;
        kite.left = kite.right;
        kite.right = tmp;
    }

    vertices[0] = kite.centre;
    vertices[1] = kite.right;
    vertices[2] = kite.outer;
    vertices[3] = kite.left;
    kite = kite_left(kite);            /* now on kite #1 */
    kite = kite_forward_right(kite);   /* now on kite #2 */
    vertices[4] = kite.centre;
    kite = kite_right(kite);           /* now on kite #3 */
    vertices[5] = kite.right;
    vertices[6] = kite.outer;
    kite = kite_forward_left(kite);    /* now on kite #4 */
    vertices[7] = kite.left;
    vertices[8] = kite.centre;
    kite = kite_right(kite);           /* now on kite #5 */
    kite = kite_right(kite);           /* now on kite #6 */
    kite = kite_right(kite);           /* now on kite #7 */
    vertices[9] = kite.right;
    vertices[10] = kite.outer;
    vertices[11] = kite.left;
    kite = kite_left(kite);            /* now on kite #6 again */
    vertices[12] = kite.outer;
    vertices[13] = kite.left;

    if (reversed) {
        /* For a reversed kite, also reverse the vertex order, so that
         * we report every polygon in a consistent orientation */
        for (i = 0, j = 13; i < j; i++, j--) {
            Point tmp = vertices[i];
            vertices[i] = vertices[j];
            vertices[j] = tmp;
        }
    }

    /*
     * Convert from our internal coordinate system into the orthogonal
     * one used in this module's external API. In the same loop, we
     * might as well do the bounds check.
     */
    for (i = 0; i < 14; i++) {
        Point v = vertices[i];
        int x = (v.x * 2 + v.y) / 3, y = v.y;

        if (x < 0 || x > 4*w || y < 0 || y > 6*h)
            return;       /* a vertex of this kite is out of bounds */

        coords[2*i] = x;
        coords[2*i+1] = y;
    }

    cb(cbctx, kite0, hc, coords);
}

struct internal_ctx {
    hat_tile_callback_fn external_cb;
    void *external_cbctx;
};
static void report_hat(void *vctx, Kite kite0, HatCoords *hc, int *coords)
{
    struct internal_ctx *ctx = (struct internal_ctx *)vctx;
    ctx->external_cb(ctx->external_cbctx, 14, coords);
}

/*
 * Generate a hat tiling from a previously generated set of parameters.
 */
void hat_tiling_generate(const struct HatPatchParams *hp, int w, int h,
                         hat_tile_callback_fn cb, void *cbctx)
{
    HatCoordContext ctx[1];
    HatCoords *coords[KE_NKEEP];
    KiteEnum s[1];
    size_t i;
    struct internal_ctx report_hat_ctx[1];

    report_hat_ctx->external_cb = cb;
    report_hat_ctx->external_cbctx = cbctx;

    init_coords_params(ctx, hp);
    for (i = 0; i < lenof(coords); i++)
        coords[i] = NULL;

    first_kite(s, w, h);
    coords[s->curr_index] = initial_coords(ctx);
    maybe_report_hat(w, h, *s->curr, coords[s->curr_index],
                     report_hat, report_hat_ctx);

    while (next_kite(s)) {
        hc_free(coords[s->curr_index]);
        coords[s->curr_index] = step_coords(
            ctx, coords[s->last_index], s->last_step);
        maybe_report_hat(w, h, *s->curr, coords[s->curr_index],
                         report_hat, report_hat_ctx);
    }

    cleanup_coords(ctx);
    for (i = 0; i < lenof(coords); i++)
        hc_free(coords[i]);
}

#ifdef TEST_HAT

#include <stdarg.h>

static HatCoords *hc_construct_v(TileType type, va_list ap)
{
    HatCoords *hc = hc_new();
    while (true) {
        int index = va_arg(ap, int);

        hc_make_space(hc, hc->nc + 1);
        hc->c[hc->nc].type = type;
        hc->c[hc->nc].index = index;
        hc->nc++;

        if (index < 0)
            return hc;

        type = va_arg(ap, TileType);
    }
}

static HatCoords *hc_construct(TileType type, ...)
{
    HatCoords *hc;
    va_list ap;

    va_start(ap, type);
    hc = hc_construct_v(type, ap);
    va_end(ap);

    return hc;
}

static bool hc_equal(HatCoords *hc1, HatCoords *hc2)
{
    size_t i;

    if (hc1->nc != hc2->nc)
        return false;

    for (i = 0; i < hc1->nc; i++) {
        if (hc1->c[i].type != hc2->c[i].type ||
            hc1->c[i].index != hc2->c[i].index)
            return false;
    }

    return true;
}

static bool hc_expect(const char *file, int line, HatCoords *hc,
                      TileType type, ...)
{
    bool equal;
    va_list ap;
    HatCoords *hce;

    va_start(ap, type);
    hce = hc_construct_v(type, ap);
    va_end(ap);

    equal = hc_equal(hc, hce);

    if (!equal) {
        fprintf(stderr, "%s:%d: coordinate mismatch\n", file, line);
        debug_coords("  expected: ", hce, "\n");
        debug_coords("  actual:   ", hc, "\n");
    }

    hc_free(hce);
    return equal;
}

#define EXPECT(hc, ...) do {                                    \
        if (!hc_expect(__FILE__, __LINE__, hc, __VA_ARGS__))    \
            fails++;                                            \
    } while (0)

/*
 * For four-colouring the tiling: these tables give a colouring of
 * each kitemap, with colour 3 assigned to the reflected tiles in the
 * middle of the H, and 0,1,2 chosen arbitrarily.
 */

static const int fourcolours_H[] = {
    /*  0 */ 0,  2,  1,  3,
    /*  1 */ 1,  0,  2,  3,
    /*  2 */ 0,  2,  1,  3,
    /*  3 */ 1, -1, -1, -1,
    /*  4 */ 1,  2, -1, -1,
    /*  5 */ 1,  2, -1, -1,
    /*  6 */ 2,  1, -1, -1,
    /*  7 */ 0,  1, -1, -1,
    /*  8 */ 2,  0, -1, -1,
    /*  9 */ 2,  0, -1, -1,
    /* 10 */ 0,  1, -1, -1,
    /* 11 */ 0,  1, -1, -1,
    /* 12 */ 2,  0, -1, -1,
};
static const int fourcolours_T[] = {
    /*  0 */ 1,  2,  0,  3,
    /*  1 */ 2,  1, -1, -1,
    /*  2 */ 0,  1, -1, -1,
    /*  3 */ 0,  2, -1, -1,
    /*  4 */ 2,  0, -1, -1,
    /*  5 */ 0,  1, -1, -1,
    /*  6 */ 1,  2, -1, -1,
};
static const int fourcolours_P[] = {
    /*  0 */ 2,  1,  0,  3,
    /*  1 */ 1,  2,  0,  3,
    /*  2 */ 2,  1, -1, -1,
    /*  3 */ 0,  2, -1, -1,
    /*  4 */ 0,  1, -1, -1,
    /*  5 */ 1,  2, -1, -1,
    /*  6 */ 2,  0, -1, -1,
    /*  7 */ 0,  1, -1, -1,
    /*  8 */ 1,  0, -1, -1,
    /*  9 */ 2,  1, -1, -1,
    /* 10 */ 0,  2, -1, -1,
};
static const int fourcolours_F[] = {
    /*  0 */ 2,  0,  1,  3,
    /*  1 */ 0,  2,  1,  3,
    /*  2 */ 1,  2, -1, -1,
    /*  3 */ 1,  0, -1, -1,
    /*  4 */ 0,  2, -1, -1,
    /*  5 */ 2,  1, -1, -1,
    /*  6 */ 2,  0, -1, -1,
    /*  7 */ 0,  1, -1, -1,
    /*  8 */ 0,  1, -1, -1,
    /*  9 */ 2,  0, -1, -1,
    /* 10 */ 1,  2, -1, -1,
};
static const int *const fourcolours[] = {
    fourcolours_H, fourcolours_T, fourcolours_P, fourcolours_F,
};

/*
 * Structure that describes how the colours in the above maps are
 * translated to output colours. This will vary with each kitemap our
 * coordinates pass through, in order to maintain consistency.
 */
typedef struct FourColourMap {
    unsigned char map[4];
} FourColourMap;

/*
 * Make an initial FourColourMap by choosing the initial permutation
 * of the three 'normal' hat colours randomly.
 */
static inline FourColourMap fourcolourmap_initial(random_state *rs)
{
    FourColourMap f;
    unsigned i;

    /* Start with the identity mapping */
    for (i = 0; i < 4; i++)
        f.map[i] = i;

    /* Randomly permute colours 0,1,2, leaving 3 as the distinguished
     * colour for reflected hats */
    shuffle(f.map, 3, sizeof(f.map[0]), rs);

    return f;
}

static inline FourColourMap fourcolourmap_update(
    FourColourMap prevm, HatCoords *prevc, HatCoords *currc, KiteStep step,
    HatCoordContext *ctx)
{
    size_t i, m1, m2;
    const int *f1, *f2;
    unsigned sum;
    int missing;
    FourColourMap newm;
    HatCoords *prev2c;

    /*
     * If prevc and currc are in the same kitemap anyway, that's the
     * easy case: the colour map for the new kitemap is the same as
     * for the old one, because they're the same kitemap.
     */
    ensure_coords(ctx, prevc, currc->nc);
    ensure_coords(ctx, currc, prevc->nc);
    for (i = 3; i < prevc->nc; i++)
        if (currc->c[i].index != prevc->c[i].index)
            goto mismatch;
    return prevm;
  mismatch:

    /*
     * The step_coords algorithm guarantees that the _new_ coordinate
     * currc is expected to be in a kitemap containing both this kite
     * and the previous one (because it first transformed the previous
     * coordinate until it _could_ take a step within the same
     * kitemap, and then did).
     *
     * So if we reverse the last step we took, we should get a second
     * HatCoords describing the same kite as prevc but showing its
     * position in the _new_ kitemap. This lets us figure out a pair
     * of corresponding metatile indices within the old and new
     * kitemaps (by looking at which metatile prevc and prev2c claim
     * to be in).
     *
     * That metatile will also always be a P or an F (because all
     * metatiles overlapping the next kitemap are of those types),
     * which means it will have two hats in it. And those hats will be
     * adjacent, so differently coloured. Hence, we have enough
     * information to decide how two of the new kitemap's three normal
     * colours map to the colours we were using in the old kitemap -
     * and then the third is determined by process of elimination.
     */
    prev2c = step_coords(
        ctx, currc, (step == KS_LEFT ? KS_RIGHT :
                     step == KS_RIGHT ? KS_LEFT :
                     step == KS_F_LEFT ? KS_F_RIGHT : KS_F_LEFT));

    /* Metatile indices within the old and new kitemaps */
    m1 = prevc->c[2].index;
    m2 = prev2c->c[2].index;

    /* The colourings of those metatiles' hats in our fixed fourcolours[] */
    f1 = fourcolours[prevc->c[3].type] + 4*m1;
    f2 = fourcolours[prev2c->c[3].type] + 4*m2;

    /*
     * Start making our new output map, filling in all three normal
     * colours to 255 = "don't know yet".
     */
    newm.map[3] = 3;
    newm.map[0] = newm.map[1] = newm.map[2] = 255;

    /*
     * Iterate over the tile colourings in fourcolours[] for these
     * metatiles, matching up our mappings.
     */
    for (i = 0; i < 4; i++) {
        /* They should be the same metatile, so have same number of hats! */
        assert((f1[i] == -1) == (f2[i] == -1));

        if (f1[i] != 255)
            newm.map[f2[i]] = prevm.map[f1[i]];
    }

    /*
     * We expect to have filled in exactly two of the three normal
     * colours. Find the missing index, and fill in its colour by
     * arithmetic (using the fact that the three colours add up to 3).
     */
    sum = 0;
    missing = -1;
    for (i = 0; i < 3; i++) {
        if (newm.map[i] == 255) {
            assert(missing == -1); /* shouldn't have two missing colours */
            missing = i;
        } else {
            sum += newm.map[i];
        }
    }
    assert(missing != -1);
    assert(0 < sum && sum <= 3);
    newm.map[missing] = 3 - sum;

    return newm;
}

static bool unit_tests(void)
{
    int fails = 0;
    HatCoordContext ctx[1];
    HatCoords *hc_in, *hc_out;

    ctx->rs = NULL;
    ctx->prototype = hc_construct(TT_KITE, 0, TT_HAT, 0, TT_H, -1);

    /* Simple steps within a hat */

    hc_in = hc_construct(TT_KITE, 6, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_out = step_coords(ctx, hc_in, KS_LEFT);
    EXPECT(hc_out, TT_KITE, 5, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    hc_in = hc_construct(TT_KITE, 6, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_out = step_coords(ctx, hc_in, KS_RIGHT);
    EXPECT(hc_out, TT_KITE, 7, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    hc_in = hc_construct(TT_KITE, 5, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_out = step_coords(ctx, hc_in, KS_F_LEFT);
    EXPECT(hc_out, TT_KITE, 2, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    hc_in = hc_construct(TT_KITE, 5, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_out = step_coords(ctx, hc_in, KS_F_RIGHT);
    EXPECT(hc_out, TT_KITE, 1, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    /* Step between hats in the same kitemap, which can change the
     * metatile type at layer 2 */

    hc_in = hc_construct(TT_KITE, 6, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_out = step_coords(ctx, hc_in, KS_F_LEFT);
    EXPECT(hc_out, TT_KITE, 3, TT_HAT, 0, TT_H, 0, TT_H, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    hc_in = hc_construct(TT_KITE, 7, TT_HAT, 2, TT_H, 1, TT_H, -1);
    hc_out = step_coords(ctx, hc_in, KS_F_RIGHT);
    EXPECT(hc_out, TT_KITE, 4, TT_HAT, 0, TT_T, 3, TT_H, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    /* Step off the edge of one kitemap, necessitating a metamap
     * rewrite of layers 2,3 to get into a different kitemap where
     * that step can be made */

    hc_in = hc_construct(TT_KITE, 6, TT_HAT, 0, TT_P, 2, TT_P, 3, TT_P, -1);
    hc_out = step_coords(ctx, hc_in, KS_F_RIGHT);
    /* Working:
     *     kite 6 . hat 0 . P 2 . P 3 . P ?
     *  -> kite 6 . hat 0 . P 6 . H 0 . P ?   (P metamap says 2.3 = 6.0)
     */
    EXPECT(hc_out, TT_KITE, 7, TT_HAT, 1, TT_H, 1, TT_H, 0, TT_P, -1);
    hc_free(hc_in);
    hc_free(hc_out);

    cleanup_coords(ctx);
    return fails == 0;
}

typedef struct pspoint {
    float x, y;
} pspoint;

typedef struct psbbox {
    bool started;
    pspoint bl, tr;
} psbbox;

static inline void psbbox_add(psbbox *bbox, pspoint p)
{
    if (!bbox->started || bbox->bl.x > p.x)
        bbox->bl.x = p.x;
    if (!bbox->started || bbox->tr.x < p.x)
        bbox->tr.x = p.x;
    if (!bbox->started || bbox->bl.y > p.y)
        bbox->bl.y = p.y;
    if (!bbox->started || bbox->tr.y < p.y)
        bbox->tr.y = p.y;
    bbox->started = true;
}

typedef enum OutFmt { OF_POSTSCRIPT, OF_PYTHON } OutFmt;
typedef enum ColourMode { CM_SEMANTIC, CM_FOURCOLOUR } ColourMode;

typedef struct drawctx {
    OutFmt outfmt;
    ColourMode colourmode;
    psbbox *bbox;
    KiteEnum *kiteenum;
    FourColourMap fourcolourmap[KE_NKEEP];
} drawctx;

static void bbox_add_hat(void *vctx, Kite kite0, HatCoords *hc, int *coords)
{
    drawctx *ctx = (drawctx *)vctx;
    pspoint p;
    size_t i;

    for (i = 0; i < 14; i++) {
        p.x = coords[2*i] * 1.5;
        p.y = coords[2*i+1] * sqrt(0.75);
        psbbox_add(ctx->bbox, p);
    }
}

static void header(drawctx *ctx)
{
    switch (ctx->outfmt) {
      case OF_POSTSCRIPT: {
        float xext = ctx->bbox->tr.x - ctx->bbox->bl.x;
        float yext = ctx->bbox->tr.y - ctx->bbox->bl.y;
        float ext = (xext > yext ? xext : yext);
        float scale = 500 / ext;
        float ox = 287 - scale * (ctx->bbox->bl.x + ctx->bbox->tr.x) / 2;
        float oy = 421 - scale * (ctx->bbox->bl.y + ctx->bbox->tr.y) / 2;

        printf("%%!PS-Adobe-2.0\n%%%%Creator: hat-test from Simon Tatham's "
               "Portable Puzzle Collection\n%%%%Pages: 1\n"
               "%%%%BoundingBox: %f %f %f %f\n"
               "%%%%EndComments\n%%%%Page: 1 1\n",
               ox + scale * ctx->bbox->bl.x - 20,
               oy + scale * ctx->bbox->bl.y - 20,
               ox + scale * ctx->bbox->tr.x + 20,
               oy + scale * ctx->bbox->tr.y + 20);

        printf("%f %f translate %f dup scale\n", ox, oy, scale);
        printf("%f setlinewidth\n", scale * 0.03);
        printf("0 setgray 1 setlinejoin 1 setlinecap\n");
        break;
      }
      default:
        break;
    }
}

static void draw_hat(void *vctx, Kite kite0, HatCoords *hc, int *coords)
{
    drawctx *ctx = (drawctx *)vctx;
    pspoint p;
    size_t i;
    int orientation;

    /*
     * Determine an index for the hat's orientation, based on the axis
     * of symmetry of its kite #0.
     */
    {
        int dx = kite0.outer.x - kite0.centre.x;
        int dy = kite0.outer.y - kite0.centre.y;
        orientation = 0;
        while (dx < 0 || dy < 0) {
            int newdx = dx + dy;
            int newdy = -dx;
            dx = newdx;
            dy = newdy;
            orientation++;
            assert(orientation < 6);
        }
    }

    switch (ctx->outfmt) {
      case OF_POSTSCRIPT: {
        const char *colour;

        printf("newpath");
        for (i = 0; i < 14; i++) {
            p.x = coords[2*i] * 1.5;
            p.y = coords[2*i+1] * sqrt(0.75);
            printf(" %f %f %s", p.x, p.y, i ? "lineto" : "moveto");
        }
        printf(" closepath gsave");

        switch (ctx->colourmode) {
          case CM_SEMANTIC:
            if (hc->c[2].type == TT_H) {
                colour = (hc->c[1].index == 3 ? "0 0.5 0.8 setrgbcolor" :
                          "0.6 0.8 1 setrgbcolor");
            } else if (hc->c[2].type == TT_F) {
                colour = "0.7 setgray";
            } else {
                colour = "1 setgray";
            }
            break;

          default /* case CM_FOURCOLOUR */: {
            /*
             * Determine the colour of this tile by translating the
             * fixed colour from fourcolours[] through our current
             * FourColourMap.
             */
            FourColourMap f = ctx->fourcolourmap[ctx->kiteenum->curr_index];
            const int *m = fourcolours[hc->c[3].type];
            static const char *const colours[] = {
                "1 0.7 0.7 setrgbcolor",
                "1 1 0.7 setrgbcolor",
                "0.7 1 0.7 setrgbcolor",
                "0.6 0.6 1 setrgbcolor",
            };
            colour = colours[f.map[m[hc->c[2].index * 4 + hc->c[1].index]]];
            break;
          }
        }
        printf(" %s fill grestore", colour);
        printf(" stroke\n");
        break;
      }
      case OF_PYTHON: {
        printf("hat('%c', %d, %d, [", "HTPF"[hc->c[2].type], hc->c[1].index,
               orientation);
        for (i = 0; i < 14; i++)
            printf("%s(%d,%d)", i ? ", " : "", coords[2*i], coords[2*i+1]);
        printf("])\n");
        break;
      }
    }
}

static void trailer(drawctx *dctx)
{
    switch (dctx->outfmt) {
      case OF_POSTSCRIPT: {
        printf("showpage\n");
        printf("%%%%Trailer\n");
        printf("%%%%EOF\n");
        break;
      }
      default:
        break;
    }
}

int main(int argc, char **argv)
{
    psbbox bbox[1];
    KiteEnum s[1];
    HatCoordContext ctx[1];
    HatCoords *coords[KE_NKEEP];
    random_state *rs;
    const char *random_seed = "12345";
    int w = 10, h = 10;
    int argpos = 0;
    size_t i;
    drawctx dctx[1];

    dctx->outfmt = OF_POSTSCRIPT;
    dctx->colourmode = CM_SEMANTIC;
    dctx->kiteenum = s;

    while (--argc > 0) {
        const char *arg = *++argv;
        if (!strcmp(arg, "--help")) {
            printf("  usage: hat-test [options] [<width>] [<height>]\n"
                   "options: --python   write a Python function call per hat\n"
                   "         --seed=STR vary the starting random seed\n"
                   "   also: hat-test --test\n");
            return 0;
        } else if (!strcmp(arg, "--test")) {
            return unit_tests() ? 0 : 1;
        } else if (!strcmp(arg, "--python")) {
            dctx->outfmt = OF_PYTHON;
        } else if (!strcmp(arg, "--fourcolour")) {
            dctx->colourmode = CM_FOURCOLOUR;
        } else if (!strncmp(arg, "--seed=", 7)) {
            random_seed = arg+7;
        } else if (arg[0] == '-') {
            fprintf(stderr, "unrecognised option '%s'\n", arg);
            return 1;
        } else {
            switch (argpos++) {
              case 0:
                w = atoi(arg);
                break;
              case 1:
                h = atoi(arg);
                break;
              default:
                fprintf(stderr, "unexpected extra argument '%s'\n", arg);
                return 1;
            }
        }
    }

    for (i = 0; i < lenof(coords); i++)
        coords[i] = NULL;

    rs = random_new(random_seed, strlen(random_seed));
    init_coords_random(ctx, rs);

    bbox->started = false;
    dctx->bbox = bbox;

    first_kite(s, w, h);
    coords[s->curr_index] = initial_coords(ctx);
    maybe_report_hat(w, h, *s->curr, coords[s->curr_index],
                     bbox_add_hat, dctx);
    while (next_kite(s)) {
        hc_free(coords[s->curr_index]);
        coords[s->curr_index] = step_coords(
            ctx, coords[s->last_index], s->last_step);
        maybe_report_hat(w, h, *s->curr, coords[s->curr_index],
                         bbox_add_hat, dctx);
    }
    for (i = 0; i < lenof(coords); i++) {
        hc_free(coords[i]);
        coords[i] = NULL;
    }

    header(dctx);

    first_kite(s, w, h);
    coords[s->curr_index] = initial_coords(ctx);
    dctx->fourcolourmap[s->curr_index] = fourcolourmap_initial(rs);
    maybe_report_hat(w, h, *s->curr, coords[s->curr_index],
                     draw_hat, dctx);
    while (next_kite(s)) {
        hc_free(coords[s->curr_index]);
        coords[s->curr_index] = step_coords(
            ctx, coords[s->last_index], s->last_step);
        dctx->fourcolourmap[s->curr_index] = fourcolourmap_update(
            dctx->fourcolourmap[s->last_index], coords[s->last_index],
            coords[s->curr_index], s->last_step, ctx);
        maybe_report_hat(w, h, *s->curr, coords[s->curr_index],
                         draw_hat, dctx);
    }
    for (i = 0; i < lenof(coords); i++) {
        hc_free(coords[i]);
        coords[i] = NULL;
    }

    trailer(dctx);

    cleanup_coords(ctx);

    return 0;
}
#endif