From James Harvey (via a period of collaborative polishing), a patch

to add two kinds of Penrose tiling to the grid types supported by
Loopy.

This has involved a certain amount of infrastructure work, because of
course the whole point of Penrose tilings is that they don't have to
be the same every time: so now grid.c has grown the capacity to
describe its grids as strings, and reconstitute them from those string
descriptions. Hence a Penrose Loopy game description consists of a
string identifying a particular piece of Penrose tiling, followed by
the normal Loopy clue encoding.

All the existing grid types decline to provide a grid description
string, so their Loopy game descriptions have not changed encoding.

[originally from svn r9159]

grid.h

@@ -9,6 +9,8 @@
 #ifndef PUZZLES_GRID_H
 #define PUZZLES_GRID_H
 
+#include "puzzles.h" /* for random_state */
+
 /* Useful macros */
 #define SQ(x) ( (x) * (x) )
 
@@ -89,22 +91,41 @@ typedef struct grid {
   int refcount;
 } grid;
 
-grid *grid_new_square(int width, int height);
+/* Grids are specified by type: GRID_SQUARE, GRID_KITE, etc. */
-grid *grid_new_honeycomb(int width, int height);
+
-grid *grid_new_triangular(int width, int height);
+#define GRIDGEN_LIST(A) \
-grid *grid_new_snubsquare(int width, int height);
+  A(SQUARE,square) \
-grid *grid_new_cairo(int width, int height);
+  A(HONEYCOMB,honeycomb) \
-grid *grid_new_greathexagonal(int width, int height);
+  A(TRIANGULAR,triangular) \
-grid *grid_new_octagonal(int width, int height);
+  A(SNUBSQUARE,snubsquare) \
-grid *grid_new_kites(int width, int height);
+  A(CAIRO,cairo) \
-grid *grid_new_floret(int width, int height);
+  A(GREATHEXAGONAL,greathexagonal) \
-grid *grid_new_dodecagonal(int width, int height);
+  A(OCTAGONAL,octagonal) \
-grid *grid_new_greatdodecagonal(int width, int height);
+  A(KITE,kites) \
+  A(FLORET,floret) \
+  A(DODECAGONAL,dodecagonal) \
+  A(GREATDODECAGONAL,greatdodecagonal) \
+  A(PENROSE_P2,penrose_p2_kite) \
+  A(PENROSE_P3,penrose_p3_thick)
+
+#define ENUM(upper,lower) GRID_ ## upper,
+typedef enum grid_type { GRIDGEN_LIST(ENUM) GRID_TYPE_MAX } grid_type;
+#undef ENUM
+
+/* Free directly after use if non-NULL. Will never contain an underscore
+ * (so clients can safely use that as a separator). */
+char *grid_new_desc(grid_type type, int width, int height, random_state *rs);
+char *grid_validate_desc(grid_type type, int width, int height, char *desc);
+
+grid *grid_new(grid_type type, int width, int height, char *desc);
 
 void grid_free(grid *g);
 
 grid_edge *grid_nearest_edge(grid *g, int x, int y);
 
+void grid_compute_size(grid_type type, int width, int height,
+                       int *tilesize, int *xextent, int *yextent);
+
 void grid_find_incentre(grid_face *f);
 
 #endif /* PUZZLES_GRID_H */

loopy.R

@@ -1,6 +1,6 @@
 # -*- makefile -*-
 
-LOOPY_EXTRA = tree234 dsf grid
+LOOPY_EXTRA = tree234 dsf grid penrose
 
 loopy     : [X] GTK COMMON loopy LOOPY_EXTRA loopy-icon|no-icon
 
@@ -9,6 +9,13 @@ loopy     : [G] WINDOWS COMMON loopy LOOPY_EXTRA loopy.res|noicon.res
 loopysolver :   [U] loopy[STANDALONE_SOLVER] LOOPY_EXTRA STANDALONE m.lib
 loopysolver :   [C] loopy[STANDALONE_SOLVER] LOOPY_EXTRA STANDALONE
 
+#penrose :    [U] penrose[TEST_PENROSE] STANDALONE m.lib
+#penrose :    [C] penrose[TEST_PENROSE] STANDALONE
+
+#test-basis : [U] penrose[TEST_VECTORS] tree234 STANDALONE m.lib
+#test-basis : [C] penrose[TEST_VECTORS] tree234 STANDALONE
+
+
 ALL += loopy[COMBINED] LOOPY_EXTRA
 
 !begin gtk

loopy.c

@@ -107,7 +107,7 @@ enum {
 };
 
 struct game_state {
-    grid *game_grid;
+    grid *game_grid; /* ref-counted (internally) */
 
     /* Put -1 in a face that doesn't get a clue */
     signed char *clues;
@@ -207,10 +207,6 @@ struct game_params {
     int w, h;
     int diff;
     int type;
-
-    /* Grid generation is expensive, so keep a (ref-counted) reference to the
-     * grid for these parameters, and only generate when required. */
-    grid *game_grid;
 };
 
 /* line_drawstate is the same as line_state, but with the extra ERROR
@@ -247,29 +243,31 @@ static void check_caches(const solver_state* sstate);
 
 /* ------- List of grid generators ------- */
 #define GRIDLIST(A) \
-    A(Squares,grid_new_square,3,3) \
+    A(Squares,GRID_SQUARE,3,3) \
-    A(Triangular,grid_new_triangular,3,3) \
+    A(Triangular,GRID_TRIANGULAR,3,3) \
-    A(Honeycomb,grid_new_honeycomb,3,3) \
+    A(Honeycomb,GRID_HONEYCOMB,3,3) \
-    A(Snub-Square,grid_new_snubsquare,3,3) \
+    A(Snub-Square,GRID_SNUBSQUARE,3,3) \
-    A(Cairo,grid_new_cairo,3,4) \
+    A(Cairo,GRID_CAIRO,3,4) \
-    A(Great-Hexagonal,grid_new_greathexagonal,3,3) \
+    A(Great-Hexagonal,GRID_GREATHEXAGONAL,3,3) \
-    A(Octagonal,grid_new_octagonal,3,3) \
+    A(Octagonal,GRID_OCTAGONAL,3,3) \
-    A(Kites,grid_new_kites,3,3) \
+    A(Kites,GRID_KITE,3,3) \
-    A(Floret,grid_new_floret,1,2) \
+    A(Floret,GRID_FLORET,1,2) \
-    A(Dodecagonal,grid_new_dodecagonal,2,2) \
+    A(Dodecagonal,GRID_DODECAGONAL,2,2) \
-    A(Great-Dodecagonal,grid_new_greatdodecagonal,2,2)
+    A(Great-Dodecagonal,GRID_GREATDODECAGONAL,2,2) \
+    A(Penrose (kite/dart),GRID_PENROSE_P2,3,3) \
+    A(Penrose (rhombs),GRID_PENROSE_P3,3,3)
 
-#define GRID_NAME(title,fn,amin,omin) #title,
+#define GRID_NAME(title,type,amin,omin) #title,
-#define GRID_CONFIG(title,fn,amin,omin) ":" #title
+#define GRID_CONFIG(title,type,amin,omin) ":" #title
-#define GRID_FN(title,fn,amin,omin) &fn,
+#define GRID_TYPE(title,type,amin,omin) type,
-#define GRID_SIZES(title,fn,amin,omin) \
+#define GRID_SIZES(title,type,amin,omin) \
     {amin, omin, \
      "Width and height for this grid type must both be at least " #amin, \
      "At least one of width and height for this grid type must be at least " #omin,},
 static char const *const gridnames[] = { GRIDLIST(GRID_NAME) };
 #define GRID_CONFIGS GRIDLIST(GRID_CONFIG)
-static grid * (*(grid_fns[]))(int w, int h) = { GRIDLIST(GRID_FN) };
+static grid_type grid_types[] = { GRIDLIST(GRID_TYPE) };
-#define NUM_GRID_TYPES (sizeof(grid_fns) / sizeof(grid_fns[0]))
+#define NUM_GRID_TYPES (sizeof(grid_types) / sizeof(grid_types[0]))
 static const struct {
     int amin, omin;
     char *aerr, *oerr;
@@ -277,13 +275,10 @@ static const struct {
 
 /* Generates a (dynamically allocated) new grid, according to the
  * type and size requested in params.  Does nothing if the grid is already
- * generated.  The allocated grid is owned by the params object, and will be
+ * generated. */
- * freed in free_params(). */
+static grid *loopy_generate_grid(game_params *params, char *grid_desc)
-static void params_generate_grid(game_params *params)
 {
-    if (!params->game_grid) {
+    return grid_new(grid_types[params->type], params->w, params->h, grid_desc);
-        params->game_grid = grid_fns[params->type](params->w, params->h);
-    }
 }
 
 /* ----------------------------------------------------------------------
@@ -480,8 +475,6 @@ static game_params *default_params(void)
     ret->diff = DIFF_EASY;
     ret->type = 0;
 
-    ret->game_grid = NULL;
-
     return ret;
 }
 
@@ -490,44 +483,45 @@ static game_params *dup_params(game_params *params)
     game_params *ret = snew(game_params);
 
     *ret = *params;                       /* structure copy */
-    if (ret->game_grid) {
-        ret->game_grid->refcount++;
-    }
     return ret;
 }
 
 static const game_params presets[] = {
 #ifdef SMALL_SCREEN
-    {  7,  7, DIFF_EASY, 0, NULL },
+    {  7,  7, DIFF_EASY, 0 },
-    {  7,  7, DIFF_NORMAL, 0, NULL },
+    {  7,  7, DIFF_NORMAL, 0 },
-    {  7,  7, DIFF_HARD, 0, NULL },
+    {  7,  7, DIFF_HARD, 0 },
-    {  7,  7, DIFF_HARD, 1, NULL },
+    {  7,  7, DIFF_HARD, 1 },
-    {  7,  7, DIFF_HARD, 2, NULL },
+    {  7,  7, DIFF_HARD, 2 },
-    {  5,  5, DIFF_HARD, 3, NULL },
+    {  5,  5, DIFF_HARD, 3 },
-    {  7,  7, DIFF_HARD, 4, NULL },
+    {  7,  7, DIFF_HARD, 4 },
-    {  5,  4, DIFF_HARD, 5, NULL },
+    {  5,  4, DIFF_HARD, 5 },
-    {  5,  5, DIFF_HARD, 6, NULL },
+    {  5,  5, DIFF_HARD, 6 },
-    {  5,  5, DIFF_HARD, 7, NULL },
+    {  5,  5, DIFF_HARD, 7 },
-    {  3,  3, DIFF_HARD, 8, NULL },
+    {  3,  3, DIFF_HARD, 8 },
-    {  3,  3, DIFF_HARD, 9, NULL },
+    {  3,  3, DIFF_HARD, 9 },
-    {  3,  3, DIFF_HARD, 10, NULL },
+    {  3,  3, DIFF_HARD, 10 },
+    {  6,  6, DIFF_HARD, 11 },
+    {  6,  6, DIFF_HARD, 12 },
 #else
-    {  7,  7, DIFF_EASY, 0, NULL },
+    {  7,  7, DIFF_EASY, 0 },
-    {  10,  10, DIFF_EASY, 0, NULL },
+    {  10,  10, DIFF_EASY, 0 },
-    {  7,  7, DIFF_NORMAL, 0, NULL },
+    {  7,  7, DIFF_NORMAL, 0 },
-    {  10,  10, DIFF_NORMAL, 0, NULL },
+    {  10,  10, DIFF_NORMAL, 0 },
-    {  7,  7, DIFF_HARD, 0, NULL },
+    {  7,  7, DIFF_HARD, 0 },
-    {  10,  10, DIFF_HARD, 0, NULL },
+    {  10,  10, DIFF_HARD, 0 },
-    {  10,  10, DIFF_HARD, 1, NULL },
+    {  10,  10, DIFF_HARD, 1 },
-    {  12,  10, DIFF_HARD, 2, NULL },
+    {  12,  10, DIFF_HARD, 2 },
-    {  7,  7, DIFF_HARD, 3, NULL },
+    {  7,  7, DIFF_HARD, 3 },
-    {  9,  9, DIFF_HARD, 4, NULL },
+    {  9,  9, DIFF_HARD, 4 },
-    {  5,  4, DIFF_HARD, 5, NULL },
+    {  5,  4, DIFF_HARD, 5 },
-    {  7,  7, DIFF_HARD, 6, NULL },
+    {  7,  7, DIFF_HARD, 6 },
-    {  5,  5, DIFF_HARD, 7, NULL },
+    {  5,  5, DIFF_HARD, 7 },
-    {  5,  5, DIFF_HARD, 8, NULL },
+    {  5,  5, DIFF_HARD, 8 },
-    {  5,  4, DIFF_HARD, 9, NULL },
+    {  5,  4, DIFF_HARD, 9 },
-    {  5,  4, DIFF_HARD, 10, NULL },
+    {  5,  4, DIFF_HARD, 10 },
+    {  10, 10, DIFF_HARD, 11 },
+    {  10, 10, DIFF_HARD, 12 }
 #endif
 };
 
@@ -551,18 +545,11 @@ static int game_fetch_preset(int i, char **name, game_params **params)
 
 static void free_params(game_params *params)
 {
-    if (params->game_grid) {
-        grid_free(params->game_grid);
-    }
     sfree(params);
 }
 
 static void decode_params(game_params *params, char const *string)
 {
-    if (params->game_grid) {
-        grid_free(params->game_grid);
-        params->game_grid = NULL;
-    }
     params->h = params->w = atoi(string);
     params->diff = DIFF_EASY;
     while (*string && isdigit((unsigned char)*string)) string++;
@@ -641,7 +628,6 @@ static game_params *custom_params(config_item *cfg)
     ret->type = cfg[2].ival;
     ret->diff = cfg[3].ival;
 
-    ret->game_grid = NULL;
     return ret;
 }
 
@@ -702,14 +688,44 @@ static char *state_to_text(const game_state *state)
     return retval;
 }
 
+#define GRID_DESC_SEP '_'
+
+/* Splits up a (optional) grid_desc from the game desc. Returns the
+ * grid_desc (which needs freeing) and updates the desc pointer to
+ * start of real desc, or returns NULL if no desc. */
+static char *extract_grid_desc(char **desc)
+{
+    char *sep = strchr(*desc, GRID_DESC_SEP), *gd;
+    int gd_len;
+
+    if (!sep) return NULL;
+
+    gd_len = sep - (*desc);
+    gd = snewn(gd_len+1, char);
+    memcpy(gd, *desc, gd_len);
+    gd[gd_len] = '\0';
+
+    *desc = sep+1;
+
+    return gd;
+}
+
 /* We require that the params pass the test in validate_params and that the
  * description fills the entire game area */
 static char *validate_desc(game_params *params, char *desc)
 {
     int count = 0;
     grid *g;
-    params_generate_grid(params);
+    char *grid_desc, *ret;
-    g = params->game_grid;
+
+    /* It's pretty inefficient to do this just for validation. All we need to
+     * know is the precise number of faces. */
+    grid_desc = extract_grid_desc(&desc);
+    ret = grid_validate_desc(grid_types[params->type], params->w, params->h, grid_desc);
+    if (ret) return ret;
+
+    g = loopy_generate_grid(params, grid_desc);
+    if (grid_desc) sfree(grid_desc);
 
     for (; *desc; ++desc) {
         if ((*desc >= '0' && *desc <= '9') || (*desc >= 'A' && *desc <= 'Z')) {
@@ -728,6 +744,8 @@ static char *validate_desc(game_params *params, char *desc)
     if (count > g->num_faces)
         return "Description too long for board size";
 
+    grid_free(g);
+
     return NULL;
 }
 
@@ -809,16 +827,15 @@ static void game_changed_state(game_ui *ui, game_state *oldstate,
 static void game_compute_size(game_params *params, int tilesize,
                               int *x, int *y)
 {
-    grid *g;
     int grid_width, grid_height, rendered_width, rendered_height;
+    int g_tilesize;
+
+    grid_compute_size(grid_types[params->type], params->w, params->h,
+                      &g_tilesize, &grid_width, &grid_height);
 
-    params_generate_grid(params);
-    g = params->game_grid;
-    grid_width = g->highest_x - g->lowest_x;
-    grid_height = g->highest_y - g->lowest_y;
     /* multiply first to minimise rounding error on integer division */
-    rendered_width = grid_width * tilesize / g->tilesize;
+    rendered_width = grid_width * tilesize / g_tilesize;
-    rendered_height = grid_height * tilesize / g->tilesize;
+    rendered_height = grid_height * tilesize / g_tilesize;
     *x = rendered_width + 2 * BORDER(tilesize) + 1;
     *y = rendered_height + 2 * BORDER(tilesize) + 1;
 }
@@ -1836,13 +1853,14 @@ static char *new_game_desc(game_params *params, random_state *rs,
                            char **aux, int interactive)
 {
     /* solution and description both use run-length encoding in obvious ways */
-    char *retval;
+    char *retval, *game_desc, *grid_desc;
     grid *g;
     game_state *state = snew(game_state);
     game_state *state_new;
-    params_generate_grid(params);
+
-    state->game_grid = g = params->game_grid;
+    grid_desc = grid_new_desc(grid_types[params->type], params->w, params->h, rs);
-    g->refcount++;
+    state->game_grid = g = loopy_generate_grid(params, grid_desc);
+
     state->clues = snewn(g->num_faces, signed char);
     state->lines = snewn(g->num_edges, char);
     state->line_errors = snewn(g->num_edges, unsigned char);
@@ -1875,10 +1893,19 @@ static char *new_game_desc(game_params *params, random_state *rs,
         goto newboard_please;
     }
 
-    retval = state_to_text(state);
+    game_desc = state_to_text(state);
 
     free_game(state);
 
+    if (grid_desc) {
+        retval = snewn(strlen(grid_desc) + 1 + strlen(game_desc) + 1, char);
+        sprintf(retval, "%s%c%s", grid_desc, GRID_DESC_SEP, game_desc);
+        sfree(grid_desc);
+        sfree(game_desc);
+    } else {
+        retval = game_desc;
+    }
+
     assert(!validate_desc(params, retval));
 
     return retval;
@@ -1890,13 +1917,17 @@ static game_state *new_game(midend *me, game_params *params, char *desc)
     game_state *state = snew(game_state);
     int empties_to_make = 0;
     int n,n2;
-    const char *dp = desc;
+    const char *dp;
+    char *grid_desc;
     grid *g;
     int num_faces, num_edges;
 
-    params_generate_grid(params);
+    grid_desc = extract_grid_desc(&desc);
-    state->game_grid = g = params->game_grid;
+    state->game_grid = g = loopy_generate_grid(params, grid_desc);
-    g->refcount++;
+    if (grid_desc) sfree(grid_desc);
+
+    dp = desc;
+
     num_faces = g->num_faces;
     num_edges = g->num_edges;
 
@@ -4017,3 +4048,5 @@ int main(int argc, char **argv)
 }
 
 #endif
+
+/* vim: set shiftwidth=4 tabstop=8: */

penrose.c

@@ -0,0 +1,626 @@
+/* penrose.c
+ *
+ * Penrose tile generator.
+ *
+ * Uses half-tile technique outlined on:
+ *
+ * http://tartarus.org/simon/20110412-penrose/penrose.xhtml
+ */
+
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "puzzles.h" /* for malloc routines, and PI */
+#include "penrose.h"
+
+/* -------------------------------------------------------
+ * 36-degree basis vector arithmetic routines.
+ */
+
+/* Imagine drawing a
+ * ten-point 'clock face' like this:
+ *
+ *                     -E
+ *                 -D   |   A
+ *                   \  |  /
+ *              -C.   \ | /   ,B
+ *                 `-._\|/_,-'
+ *                 ,-' /|\ `-.
+ *              -B'   / | \   `C
+ *                   /  |  \
+ *                 -A   |   D
+ *                      E
+ *
+ * In case the ASCII art isn't clear, those are supposed to be ten
+ * vectors of length 1, all sticking out from the origin at equal
+ * angular spacing (hence 36 degrees). Our basis vectors are A,B,C,D (I
+ * choose them to be symmetric about the x-axis so that the final
+ * translation into 2d coordinates will also be symmetric, which I
+ * think will avoid minor rounding uglinesses), so our vector
+ * representation sets
+ *
+ *   A = (1,0,0,0)
+ *   B = (0,1,0,0)
+ *   C = (0,0,1,0)
+ *   D = (0,0,0,1)
+ *
+ * The fifth vector E looks at first glance as if it needs to be
+ * another basis vector, but in fact it doesn't, because it can be
+ * represented in terms of the other four. Imagine starting from the
+ * origin and following the path -A, +B, -C, +D: you'll find you've
+ * traced four sides of a pentagram, and ended up one E-vector away
+ * from the origin. So we have
+ *
+ *   E = (-1,1,-1,1)
+ *
+ * This tells us that we can rotate any vector in this system by 36
+ * degrees: if we start with a*A + b*B + c*C + d*D, we want to end up
+ * with a*B + b*C + c*D + d*E, and we substitute our identity for E to
+ * turn that into a*B + b*C + c*D + d*(-A+B-C+D). In other words,
+ *
+ *   rotate_one_notch_clockwise(a,b,c,d) = (-d, d+a, -d+b, d+c)
+ *
+ * and you can verify for yourself that applying that operation
+ * repeatedly starting with (1,0,0,0) cycles round ten vectors and
+ * comes back to where it started.
+ *
+ * The other operation that may be required is to construct vectors
+ * with lengths that are multiples of phi. That can be done by
+ * observing that the vector C-B is parallel to E and has length 1/phi,
+ * and the vector D-A is parallel to E and has length phi. So this
+ * tells us that given any vector, we can construct one which points in
+ * the same direction and is 1/phi or phi times its length, like this:
+ *
+ *   divide_by_phi(vector) = rotate(vector, 2) - rotate(vector, 3)
+ *   multiply_by_phi(vector) = rotate(vector, 1) - rotate(vector, 4)
+ *
+ * where rotate(vector, n) means applying the above
+ * rotate_one_notch_clockwise primitive n times. Expanding out the
+ * applications of rotate gives the following direct representation in
+ * terms of the vector coordinates:
+ *
+ *   divide_by_phi(a,b,c,d) = (b-d, c+d-b, a+b-c, c-a)
+ *   multiply_by_phi(a,b,c,d) = (a+b-d, c+d, a+b, c+d-a)
+ *
+ * and you can verify for yourself that those two operations are
+ * inverses of each other (as you'd hope!).
+ *
+ * Having done all of this, testing for equality between two vectors is
+ * a trivial matter of comparing the four integer coordinates. (Which
+ * it _wouldn't_ have been if we'd kept E as a fifth basis vector,
+ * because then (-1,1,-1,1,0) and (0,0,0,0,1) would have had to be
+ * considered identical. So leaving E out is vital.)
+ */
+
+struct vector { int a, b, c, d; };
+
+static vector v_origin()
+{
+    vector v;
+    v.a = v.b = v.c = v.d = 0;
+    return v;
+}
+
+/* We start with a unit vector of B: this means we can easily
+ * draw an isoceles triangle centred on the X axis. */
+#ifdef TEST_VECTORS
+
+static vector v_unit()
+{
+    vector v;
+
+    v.b = 1;
+    v.a = v.c = v.d = 0;
+    return v;
+}
+
+#endif
+
+#define COS54 0.5877852
+#define SIN54 0.8090169
+#define COS18 0.9510565
+#define SIN18 0.3090169
+
+/* These two are a bit rough-and-ready for now. Note that B/C are
+ * 18 degrees from the x-axis, and A/D are 54 degrees. */
+double v_x(vector *vs, int i)
+{
+    return (vs[i].a + vs[i].d) * COS54 +
+           (vs[i].b + vs[i].c) * COS18;
+}
+
+double v_y(vector *vs, int i)
+{
+    return (vs[i].a - vs[i].d) * SIN54 +
+           (vs[i].b - vs[i].c) * SIN18;
+
+}
+
+static vector v_trans(vector v, vector trans)
+{
+    v.a += trans.a;
+    v.b += trans.b;
+    v.c += trans.c;
+    v.d += trans.d;
+    return v;
+}
+
+static vector v_rotate_36(vector v)
+{
+    vector vv;
+    vv.a = -v.d;
+    vv.b =  v.d + v.a;
+    vv.c = -v.d + v.b;
+    vv.d =  v.d + v.c;
+    return vv;
+}
+
+static vector v_rotate(vector v, int ang)
+{
+    int i;
+
+    assert((ang % 36) == 0);
+    while (ang < 0) ang += 360;
+    ang = 360-ang;
+    for (i = 0; i < (ang/36); i++)
+        v = v_rotate_36(v);
+    return v;
+}
+
+#ifdef TEST_VECTORS
+
+static vector v_scale(vector v, int sc)
+{
+    v.a *= sc;
+    v.b *= sc;
+    v.c *= sc;
+    v.d *= sc;
+    return v;
+}
+
+#endif
+
+static vector v_growphi(vector v)
+{
+    vector vv;
+    vv.a = v.a + v.b - v.d;
+    vv.b = v.c + v.d;
+    vv.c = v.a + v.b;
+    vv.d = v.c + v.d - v.a;
+    return vv;
+}
+
+static vector v_shrinkphi(vector v)
+{
+    vector vv;
+    vv.a = v.b - v.d;
+    vv.b = v.c + v.d - v.b;
+    vv.c = v.a + v.b - v.c;
+    vv.d = v.c - v.a;
+    return vv;
+}
+
+#ifdef TEST_VECTORS
+
+static const char *v_debug(vector v)
+{
+    static char buf[255];
+    sprintf(buf,
+             "(%d,%d,%d,%d)[%2.2f,%2.2f]",
+             v.a, v.b, v.c, v.d, v_x(&v,0), v_y(&v,0));
+    return buf;
+}
+
+#endif
+
+/* -------------------------------------------------------
+ * Tiling routines.
+ */
+
+vector xform_coord(vector vo, int shrink, vector vtrans, int ang)
+{
+    if (shrink < 0)
+        vo = v_shrinkphi(vo);
+    else if (shrink > 0)
+        vo = v_growphi(vo);
+
+    vo = v_rotate(vo, ang);
+    vo = v_trans(vo, vtrans);
+
+    return vo;
+}
+
+
+#define XFORM(n,o,s,a) vs[(n)] = xform_coord(v_edge, (s), vs[(o)], (a))
+
+static int penrose_p2_small(penrose_state *state, int depth, int flip,
+                            vector v_orig, vector v_edge);
+
+static int penrose_p2_large(penrose_state *state, int depth, int flip,
+                            vector v_orig, vector v_edge)
+{
+    vector vv_orig, vv_edge;
+
+#ifdef DEBUG_PENROSE
+    {
+        vector vs[3];
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, 0);
+        XFORM(2, 0, 0, -36*flip);
+
+        state->new_tile(state, vs, 3, depth);
+    }
+#endif
+
+    if (flip > 0) {
+        vector vs[4];
+
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, -36);
+        XFORM(2, 0, 0, 0);
+        XFORM(3, 0, 0, 36);
+
+        state->new_tile(state, vs, 4, depth);
+    }
+    if (depth >= state->max_depth) return 0;
+
+    vv_orig = v_trans(v_orig, v_rotate(v_edge, -36*flip));
+    vv_edge = v_rotate(v_edge, 108*flip);
+
+    penrose_p2_small(state, depth+1, flip,
+                     v_orig, v_shrinkphi(v_edge));
+
+    penrose_p2_large(state, depth+1, flip,
+                     vv_orig, v_shrinkphi(vv_edge));
+    penrose_p2_large(state, depth+1, -flip,
+                     vv_orig, v_shrinkphi(vv_edge));
+
+    return 0;
+}
+
+static int penrose_p2_small(penrose_state *state, int depth, int flip,
+                            vector v_orig, vector v_edge)
+{
+    vector vv_orig;
+
+#ifdef DEBUG_PENROSE
+    {
+        vector vs[3];
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, 0);
+        XFORM(2, 0, -1, -36*flip);
+
+        state->new_tile(state, vs, 3, depth);
+    }
+#endif
+
+    if (flip > 0) {
+        vector vs[4];
+
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, -72);
+        XFORM(2, 0, -1, -36);
+        XFORM(3, 0, 0, 0);
+
+        state->new_tile(state, vs, 4, depth);
+    }
+
+    if (depth >= state->max_depth) return 0;
+
+    vv_orig = v_trans(v_orig, v_edge);
+
+    penrose_p2_large(state, depth+1, -flip,
+                     v_orig, v_shrinkphi(v_rotate(v_edge, -36*flip)));
+
+    penrose_p2_small(state, depth+1, flip,
+                     vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
+
+    return 0;
+}
+
+static int penrose_p3_small(penrose_state *state, int depth, int flip,
+                            vector v_orig, vector v_edge);
+
+static int penrose_p3_large(penrose_state *state, int depth, int flip,
+                            vector v_orig, vector v_edge)
+{
+    vector vv_orig;
+
+#ifdef DEBUG_PENROSE
+    {
+        vector vs[3];
+        vs[0] = v_orig;
+        XFORM(1, 0, 1, 0);
+        XFORM(2, 0, 0, -36*flip);
+
+        state->new_tile(state, vs, 3, depth);
+    }
+#endif
+
+    if (flip > 0) {
+        vector vs[4];
+
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, -36);
+        XFORM(2, 0, 1, 0);
+        XFORM(3, 0, 0, 36);
+
+        state->new_tile(state, vs, 4, depth);
+    }
+    if (depth >= state->max_depth) return 0;
+
+    vv_orig = v_trans(v_orig, v_edge);
+
+    penrose_p3_large(state, depth+1, -flip,
+                     vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
+
+    penrose_p3_small(state, depth+1, flip,
+                     vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
+
+    vv_orig = v_trans(v_orig, v_growphi(v_edge));
+
+    penrose_p3_large(state, depth+1, flip,
+                     vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
+
+
+    return 0;
+}
+
+static int penrose_p3_small(penrose_state *state, int depth, int flip,
+                            vector v_orig, vector v_edge)
+{
+    vector vv_orig;
+
+#ifdef DEBUG_PENROSE
+    {
+        vector vs[3];
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, 0);
+        XFORM(2, 0, 0, -36*flip);
+
+        state->new_tile(state, vs, 3, depth);
+    }
+#endif
+
+    if (flip > 0) {
+        vector vs[4];
+
+        vs[0] = v_orig;
+        XFORM(1, 0, 0, -36);
+        XFORM(3, 0, 0, 0);
+        XFORM(2, 3, 0, -36);
+
+        state->new_tile(state, vs, 4, depth);
+    }
+    if (depth >= state->max_depth) return 0;
+
+    /* NB these two are identical to the first two of p3_large. */
+    vv_orig = v_trans(v_orig, v_edge);
+
+    penrose_p3_large(state, depth+1, -flip,
+                     vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
+
+    penrose_p3_small(state, depth+1, flip,
+                     vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
+
+    return 0;
+}
+
+/* -------------------------------------------------------
+ * Utility routines.
+ */
+
+double penrose_side_length(double start_size, int depth)
+{
+  return start_size / pow(PHI, depth);
+}
+
+void penrose_count_tiles(int depth, int *nlarge, int *nsmall)
+{
+  /* Steal sgt's fibonacci thingummy. */
+}
+
+/*
+ * It turns out that an acute isosceles triangle with sides in ratio 1:phi:phi
+ * has an incentre which is conveniently 2*phi^-2 of the way from the apex to
+ * the base. Why's that convenient? Because: if we situate the incentre of the
+ * triangle at the origin, then we can place the apex at phi^-2 * (B+C), and
+ * the other two vertices at apex-B and apex-C respectively. So that's an acute
+ * triangle with its long sides of unit length, covering a circle about the
+ * origin of radius 1-(2*phi^-2), which is conveniently enough phi^-3.
+ *
+ * (later mail: this is an overestimate by about 5%)
+ */
+
+int penrose(penrose_state *state, int which)
+{
+    vector vo = v_origin();
+    vector vb = v_origin();
+
+    vo.b = vo.c = -state->start_size;
+    vo = v_shrinkphi(v_shrinkphi(vo));
+
+    vb.b = state->start_size;
+
+    if (which == PENROSE_P2)
+        return penrose_p2_large(state, 0, 1, vo, vb);
+    else
+        return penrose_p3_small(state, 0, 1, vo, vb);
+}
+
+/*
+ * We're asked for a MxN grid, which just means a tiling fitting into roughly
+ * an MxN space in some kind of reasonable unit - say, the side length of the
+ * two-arrow edges of the tiles. By some reasoning in a previous email, that
+ * means we want to pick some subarea of a circle of radius 3.11*sqrt(M^2+N^2).
+ * To cover that circle, we need to subdivide a triangle large enough that it
+ * contains a circle of that radius.
+ *
+ * Hence: start with those three vectors marking triangle vertices, scale them
+ * all up by phi repeatedly until the radius of the inscribed circle gets
+ * bigger than the target, and then recurse into that triangle with the same
+ * recursion depth as the number of times you scaled up. That will give you
+ * tiles of unit side length, covering a circle big enough that if you randomly
+ * choose an orientation and coordinates within the circle, you'll be able to
+ * get any valid piece of Penrose tiling of size MxN.
+ */
+#define INCIRCLE_RADIUS 0.22426 /* phi^-3 less 5%: see above */
+
+void penrose_calculate_size(int which, int tilesize, int w, int h,
+                            double *required_radius, int *start_size, int *depth)
+{
+    double rradius, size;
+    int n = 0;
+
+    /*
+     * Fudge factor to scale P2 and P3 tilings differently. This
+     * doesn't seem to have much relevance to questions like the
+     * average number of tiles per unit area; it's just aesthetic.
+     */
+    if (which == PENROSE_P2)
+        tilesize = tilesize * 3 / 2;
+    else
+        tilesize = tilesize * 5 / 4;
+
+    rradius = tilesize * 3.11 * sqrt((double)(w*w + h*h));
+    size = tilesize;
+
+    while ((size * INCIRCLE_RADIUS) < rradius) {
+        n++;
+        size = size * PHI;
+    }
+
+    *start_size = (int)size;
+    *depth = n;
+    *required_radius = rradius;
+}
+
+/* -------------------------------------------------------
+ * Test code.
+ */
+
+#ifdef TEST_PENROSE
+
+#include <stdio.h>
+#include <string.h>
+
+int show_recursion = 0;
+int ntiles, nfinal;
+
+int test_cb(penrose_state *state, vector *vs, int n, int depth)
+{
+    int i, xoff = 0, yoff = 0;
+    double l = penrose_side_length(state->start_size, depth);
+    double rball = l / 10.0;
+    const char *col;
+
+    ntiles++;
+    if (state->max_depth == depth) {
+        col = n == 4 ? "black" : "green";
+        nfinal++;
+    } else {
+        if (!show_recursion)
+            return 0;
+        col = n == 4 ? "red" : "blue";
+    }
+    if (n != 4) yoff = state->start_size;
+
+    printf("<polygon points=\"");
+    for (i = 0; i < n; i++) {
+        printf("%s%f,%f", (i == 0) ? "" : " ",
+               v_x(vs, i) + xoff, v_y(vs, i) + yoff);
+    }
+    printf("\" style=\"fill: %s; fill-opacity: 0.2; stroke: %s\" />\n", col, col);
+    printf("<ellipse cx=\"%f\" cy=\"%f\" rx=\"%f\" ry=\"%f\" fill=\"%s\" />",
+           v_x(vs, 0) + xoff, v_y(vs, 0) + yoff, rball, rball, col);
+
+    return 0;
+}
+
+void usage_exit()
+{
+    fprintf(stderr, "Usage: penrose-test [--recursion] P2|P3 SIZE DEPTH\n");
+    exit(1);
+}
+
+int main(int argc, char *argv[])
+{
+    penrose_state ps;
+    int which = 0;
+
+    while (--argc > 0) {
+        char *p = *++argv;
+        if (!strcmp(p, "-h") || !strcmp(p, "--help")) {
+            usage_exit();
+        } else if (!strcmp(p, "--recursion")) {
+            show_recursion = 1;
+        } else if (*p == '-') {
+            fprintf(stderr, "Unrecognised option '%s'\n", p);
+            exit(1);
+        } else {
+            break;
+        }
+    }
+
+    if (argc < 3) usage_exit();
+
+    if (strcmp(argv[0], "P2") == 0) which = PENROSE_P2;
+    else if (strcmp(argv[0], "P3") == 0) which = PENROSE_P3;
+    else usage_exit();
+
+    ps.start_size = atoi(argv[1]);
+    ps.max_depth = atoi(argv[2]);
+    ps.new_tile = test_cb;
+
+    ntiles = nfinal = 0;
+
+    printf("\
+<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n\
+<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 20010904//EN\"\n\
+\"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\
+\n\
+<svg xmlns=\"http://www.w3.org/2000/svg\"\n\
+xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n\n");
+
+    printf("<g>\n");
+    penrose(&ps, which);
+    printf("</g>\n");
+
+    printf("<!-- %d tiles and %d leaf tiles total -->\n",
+           ntiles, nfinal);
+
+    printf("</svg>");
+
+    return 0;
+}
+
+#endif
+
+#ifdef TEST_VECTORS
+
+static void dbgv(const char *msg, vector v)
+{
+    printf("%s: %s\n", msg, v_debug(v));
+}
+
+int main(int argc, const char *argv[])
+{
+   vector v = v_unit();
+
+   dbgv("unit vector", v);
+   v = v_rotate(v, 36);
+   dbgv("rotated 36", v);
+   v = v_scale(v, 2);
+   dbgv("scaled x2", v);
+   v = v_shrinkphi(v);
+   dbgv("shrunk phi", v);
+   v = v_rotate(v, -36);
+   dbgv("rotated -36", v);
+
+   return 0;
+}
+
+#endif
+/* vim: set shiftwidth=4 tabstop=8: */

penrose.h

@@ -0,0 +1,59 @@
+/* penrose.h
+ *
+ * Penrose tiling functions.
+ *
+ * Provides an interface with which to generate Penrose tilings
+ * by recursive subdivision of an initial tile of choice (one of the
+ * four sets of two pairs kite/dart, or thin/thick rhombus).
+ *
+ * You supply a callback function and a context pointer, which is
+ * called with each tile in turn: you choose how many times to recurse.
+ */
+
+#ifndef _PENROSE_H
+#define _PENROSE_H
+
+#ifndef PHI
+#define PHI 1.6180339887
+#endif
+
+typedef struct vector vector;
+
+double v_x(vector *vs, int i);
+double v_y(vector *vs, int i);
+
+typedef struct penrose_state penrose_state;
+
+/* Return non-zero to clip the tree here (i.e. not recurse
+ * below this tile).
+ *
+ * Parameters are state, vector array, npoints, depth.
+ * ctx is inside state.
+ */
+typedef int (*tile_callback)(penrose_state *, vector *, int, int);
+
+struct penrose_state {
+    int start_size;  /* initial side length */
+    int max_depth;      /* Recursion depth */
+
+    tile_callback new_tile;
+    void *ctx;          /* for callback */
+};
+
+enum { PENROSE_P2, PENROSE_P3 };
+
+extern int penrose(penrose_state *state, int which);
+
+/* Returns the side-length of a penrose tile at recursion level
+ * gen, given a starting side length. */
+extern double penrose_side_length(double start_size, int depth);
+
+/* Returns the count of each type of tile at a given recursion depth. */
+extern void penrose_count_tiles(int gen, int *nlarge, int *nsmall);
+
+/* Calculate start size and recursion depth required to produce a
+ * width-by-height sized patch of penrose tiles with the given tilesize */
+extern void penrose_calculate_size(int which, int tilesize, int w, int h,
+                                   double *required_radius, int *start_size, int *depth);
+
+#endif