From James Harvey (via a period of collaborative polishing), a patch
authorsimon <simon@cda61777-01e9-0310-a592-d414129be87e>
Sun, 24 Apr 2011 09:10:52 +0000 (09:10 +0000)
committersimon <simon@cda61777-01e9-0310-a592-d414129be87e>
Sun, 24 Apr 2011 09:10:52 +0000 (09:10 +0000)
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.

git-svn-id: svn://svn.tartarus.org/sgt/puzzles@9159 cda61777-01e9-0310-a592-d414129be87e

grid.c
grid.h
loopy.R
loopy.c
penrose.c [new file with mode: 0644]
penrose.h [new file with mode: 0644]

diff --git a/grid.c b/grid.c
index 9bcbce7..f7c6a5c 100644 (file)
--- a/grid.c
+++ b/grid.c
 #include <assert.h>
 #include <ctype.h>
 #include <math.h>
+#include <errno.h>
 
 #include "puzzles.h"
 #include "tree234.h"
 #include "grid.h"
+#include "penrose.h"
 
 /* Debugging options */
 
@@ -50,7 +52,7 @@ void grid_free(grid *g)
 
 /* Used by the other grid generators.  Create a brand new grid with nothing
  * initialised (all lists are NULL) */
-static grid *grid_new(void)
+static grid *grid_empty()
 {
     grid *g = snew(grid);
     g->faces = NULL;
@@ -158,13 +160,93 @@ grid_edge *grid_nearest_edge(grid *g, int x, int y)
  * Grid generation
  */
 
-#ifdef DEBUG_GRID
+#ifdef SVG_GRID
+
+#define SVG_DOTS  1
+#define SVG_EDGES 2
+#define SVG_FACES 4
+
+#define FACE_COLOUR "red"
+#define EDGE_COLOUR "blue"
+#define DOT_COLOUR "black"
+
+static void grid_output_svg(FILE *fp, grid *g, int which)
+{
+    int i, j;
+
+    fprintf(fp,"\
+<?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");
+
+    if (which & SVG_FACES) {
+        fprintf(fp, "<g>\n");
+        for (i = 0; i < g->num_faces; i++) {
+            grid_face *f = g->faces + i;
+            fprintf(fp, "<polygon points=\"");
+            for (j = 0; j < f->order; j++) {
+                grid_dot *d = f->dots[j];
+                fprintf(fp, "%s%d,%d", (j == 0) ? "" : " ",
+                        d->x, d->y);
+            }
+            fprintf(fp, "\" style=\"fill: %s; fill-opacity: 0.2; stroke: %s\" />\n",
+                    FACE_COLOUR, FACE_COLOUR);
+        }
+        fprintf(fp, "</g>\n");
+    }
+    if (which & SVG_EDGES) {
+        fprintf(fp, "<g>\n");
+        for (i = 0; i < g->num_edges; i++) {
+            grid_edge *e = g->edges + i;
+            grid_dot *d1 = e->dot1, *d2 = e->dot2;
+
+            fprintf(fp, "<line x1=\"%d\" y1=\"%d\" x2=\"%d\" y2=\"%d\" "
+                        "style=\"stroke: %s\" />\n",
+                        d1->x, d1->y, d2->x, d2->y, EDGE_COLOUR);
+        }
+        fprintf(fp, "</g>\n");
+    }
+
+    if (which & SVG_DOTS) {
+        fprintf(fp, "<g>\n");
+        for (i = 0; i < g->num_dots; i++) {
+            grid_dot *d = g->dots + i;
+            fprintf(fp, "<ellipse cx=\"%d\" cy=\"%d\" rx=\"%d\" ry=\"%d\" fill=\"%s\" />",
+                    d->x, d->y, g->tilesize/20, g->tilesize/20, DOT_COLOUR);
+        }
+        fprintf(fp, "</g>\n");
+    }
+
+    fprintf(fp, "</svg>\n");
+}
+#endif
+
+#ifdef SVG_GRID
+static void grid_try_svg(grid *g, int which)
+{
+    char *svg = getenv("PUZZLES_SVG_GRID");
+    if (svg) {
+        FILE *svgf = fopen(svg, "w");
+        if (svgf) {
+            grid_output_svg(svgf, g, which);
+            fclose(svgf);
+        } else {
+            fprintf(stderr, "Unable to open file `%s': %s", svg, strerror(errno));
+        }
+    }
+}
+#endif
+
 /* Show the basic grid information, before doing grid_make_consistent */
-static void grid_print_basic(grid *g)
+static void grid_debug_basic(grid *g)
 {
     /* TODO: Maybe we should generate an SVG image of the dots and lines
      * of the grid here, before grid_make_consistent.
      * Would help with debugging grid generation. */
+#ifdef DEBUG_GRID
     int i;
     printf("--- Basic Grid Data ---\n");
     for (i = 0; i < g->num_faces; i++) {
@@ -177,10 +259,16 @@ static void grid_print_basic(grid *g)
         }
         printf("]\n");
     }
+#endif
+#ifdef SVG_GRID
+    grid_try_svg(g, SVG_FACES);
+#endif
 }
+
 /* Show the derived grid information, computed by grid_make_consistent */
-static void grid_print_derived(grid *g)
+static void grid_debug_derived(grid *g)
 {
+#ifdef DEBUG_GRID
     /* edges */
     int i;
     printf("--- Derived Grid Data ---\n");
@@ -220,8 +308,11 @@ static void grid_print_derived(grid *g)
         }
         printf("]\n");
     }
+#endif
+#ifdef SVG_GRID
+    grid_try_svg(g, SVG_DOTS | SVG_EDGES | SVG_FACES);
+#endif
 }
-#endif /* DEBUG_GRID */
 
 /* Helper function for building incomplete-edges list in
  * grid_make_consistent() */
@@ -250,6 +341,147 @@ static int grid_edge_bydots_cmpfn(void *v1, void *v2)
     return 0;
 }
 
+/*
+ * 'Vigorously trim' a grid, by which I mean deleting any isolated or
+ * uninteresting faces. By which, in turn, I mean: ensure that the
+ * grid is composed solely of faces adjacent to at least one
+ * 'landlocked' dot (i.e. one not in contact with the infinite
+ * exterior face), and that all those dots are in a single connected
+ * component.
+ *
+ * This function operates on, and returns, a grid satisfying the
+ * preconditions to grid_make_consistent() rather than the
+ * postconditions. (So call it first.)
+ */
+static void grid_trim_vigorously(grid *g)
+{
+    int *dotpairs, *faces, *dots;
+    int *dsf;
+    int i, j, k, size, newfaces, newdots;
+
+    /*
+     * First construct a matrix in which each ordered pair of dots is
+     * mapped to the index of the face in which those dots occur in
+     * that order.
+     */
+    dotpairs = snewn(g->num_dots * g->num_dots, int);
+    for (i = 0; i < g->num_dots; i++)
+        for (j = 0; j < g->num_dots; j++)
+            dotpairs[i*g->num_dots+j] = -1;
+    for (i = 0; i < g->num_faces; i++) {
+        grid_face *f = g->faces + i;
+        int dot0 = f->dots[f->order-1] - g->dots;
+        for (j = 0; j < f->order; j++) {
+            int dot1 = f->dots[j] - g->dots;
+            dotpairs[dot0 * g->num_dots + dot1] = i;
+            dot0 = dot1;
+        }
+    }
+
+    /*
+     * Now we can identify landlocked dots: they're the ones all of
+     * whose edges have a mirror-image counterpart in this matrix.
+     */
+    dots = snewn(g->num_dots, int);
+    for (i = 0; i < g->num_dots; i++) {
+        dots[i] = TRUE;
+        for (j = 0; j < g->num_dots; j++) {
+            if ((dotpairs[i*g->num_dots+j] >= 0) ^
+                (dotpairs[j*g->num_dots+i] >= 0))
+                dots[i] = FALSE;    /* non-duplicated edge: coastal dot */
+        }
+    }
+
+    /*
+     * Now identify connected pairs of landlocked dots, and form a dsf
+     * unifying them.
+     */
+    dsf = snew_dsf(g->num_dots);
+    for (i = 0; i < g->num_dots; i++)
+        for (j = 0; j < i; j++)
+            if (dots[i] && dots[j] &&
+                dotpairs[i*g->num_dots+j] >= 0 &&
+                dotpairs[j*g->num_dots+i] >= 0)
+                dsf_merge(dsf, i, j);
+
+    /*
+     * Now look for the largest component.
+     */
+    size = 0;
+    j = -1;
+    for (i = 0; i < g->num_dots; i++) {
+        int newsize;
+        if (dots[i] && dsf_canonify(dsf, i) == i &&
+            (newsize = dsf_size(dsf, i)) > size) {
+            j = i;
+            size = newsize;
+        }
+    }
+
+    /*
+     * Work out which faces we're going to keep (precisely those with
+     * at least one dot in the same connected component as j) and
+     * which dots (those required by any face we're keeping).
+     *
+     * At this point we reuse the 'dots' array to indicate the dots
+     * we're keeping, rather than the ones that are landlocked.
+     */
+    faces = snewn(g->num_faces, int);
+    for (i = 0; i < g->num_faces; i++)
+        faces[i] = 0;
+    for (i = 0; i < g->num_dots; i++)
+        dots[i] = 0;
+    for (i = 0; i < g->num_faces; i++) {
+        grid_face *f = g->faces + i;
+        int keep = FALSE;
+        for (k = 0; k < f->order; k++)
+            if (dsf_canonify(dsf, f->dots[k] - g->dots) == j)
+                keep = TRUE;
+        if (keep) {
+            faces[i] = TRUE;
+            for (k = 0; k < f->order; k++)
+                dots[f->dots[k]-g->dots] = TRUE;
+        }
+    }
+
+    /*
+     * Work out the new indices of those faces and dots, when we
+     * compact the arrays containing them.
+     */
+    for (i = newfaces = 0; i < g->num_faces; i++)
+        faces[i] = (faces[i] ? newfaces++ : -1);
+    for (i = newdots = 0; i < g->num_dots; i++)
+        dots[i] = (dots[i] ? newdots++ : -1);
+
+    /*
+     * Go through and compact the arrays.
+     */
+    for (i = 0; i < g->num_dots; i++)
+        if (dots[i] >= 0) {
+            grid_dot *dnew = g->dots + dots[i], *dold = g->dots + i;
+            *dnew = *dold;             /* structure copy */
+        }
+    for (i = 0; i < g->num_faces; i++)
+        if (faces[i] >= 0) {
+            grid_face *fnew = g->faces + faces[i], *fold = g->faces + i;
+            *fnew = *fold;             /* structure copy */
+            for (j = 0; j < fnew->order; j++) {
+                /*
+                 * Reindex the dots in this face.
+                 */
+                k = fnew->dots[j] - g->dots;
+                fnew->dots[j] = g->dots + dots[k];
+            }
+        }
+    g->num_faces = newfaces;
+    g->num_dots = newdots;
+
+    sfree(dotpairs);
+    sfree(dsf);
+    sfree(dots);
+    sfree(faces);
+}
+
 /* Input: grid has its dots and faces initialised:
  * - dots have (optionally) x and y coordinates, but no edges or faces
  * (pointers are NULL).
@@ -265,9 +497,7 @@ static void grid_make_consistent(grid *g)
     tree234 *incomplete_edges;
     grid_edge *next_new_edge; /* Where new edge will go into g->edges */
 
-#ifdef DEBUG_GRID
-    grid_print_basic(g);
-#endif
+    grid_debug_basic(g);
 
     /* ====== Stage 1 ======
      * Generate edges
@@ -545,10 +775,8 @@ static void grid_make_consistent(grid *g)
             g->highest_y = max(g->highest_y, d->y);
         }
     }
-    
-#ifdef DEBUG_GRID
-    grid_print_derived(g);
-#endif
+
+    grid_debug_derived(g);
 }
 
 /* Helpers for making grid-generation easier.  These functions are only
@@ -1139,11 +1367,23 @@ void grid_find_incentre(grid_face *f)
  * arithmetic here!
  */
 
-grid *grid_new_square(int width, int height)
+#define SQUARE_TILESIZE 20
+
+void grid_size_square(int width, int height,
+                      int *tilesize, int *xextent, int *yextent)
+{
+    int a = SQUARE_TILESIZE;
+
+    *tilesize = a;
+    *xextent = width * a;
+    *yextent = height * a;
+}
+
+grid *grid_new_square(int width, int height, char *desc)
 {
     int x, y;
     /* Side length */
-    int a = 20;
+    int a = SQUARE_TILESIZE;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = width * height;
@@ -1151,7 +1391,7 @@ grid *grid_new_square(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
+    grid *g = grid_empty();
     g->tilesize = a;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
@@ -1186,21 +1426,36 @@ grid *grid_new_square(int width, int height)
     return g;
 }
 
-grid *grid_new_honeycomb(int width, int height)
+#define HONEY_TILESIZE 45
+/* Vector for side of hexagon - ratio is close to sqrt(3) */
+#define HONEY_A 15
+#define HONEY_B 26
+
+void grid_size_honeycomb(int width, int height,
+                         int *tilesize, int *xextent, int *yextent)
+{
+    int a = HONEY_A;
+    int b = HONEY_B;
+
+    *tilesize = HONEY_TILESIZE;
+    *xextent = (3 * a * (width-1)) + 4*a;
+    *yextent = (2 * b * (height-1)) + 3*b;
+}
+
+grid *grid_new_honeycomb(int width, int height, char *desc)
 {
     int x, y;
-    /* Vector for side of hexagon - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = HONEY_A;
+    int b = HONEY_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = width * height;
     int max_dots = 2 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 3 * a;
+    grid *g = grid_empty();
+    g->tilesize = HONEY_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1240,23 +1495,39 @@ grid *grid_new_honeycomb(int width, int height)
     return g;
 }
 
+#define TRIANGLE_TILESIZE 18
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define TRIANGLE_VEC_X 15
+#define TRIANGLE_VEC_Y 26
+
+void grid_size_triangular(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int vec_x = TRIANGLE_VEC_X;
+    int vec_y = TRIANGLE_VEC_Y;
+
+    *tilesize = TRIANGLE_TILESIZE;
+    *xextent = width * 2 * vec_x + vec_x;
+    *yextent = height * vec_y;
+}
+
 /* Doesn't use the previous method of generation, it pre-dates it!
  * A triangular grid is just about simple enough to do by "brute force" */
-grid *grid_new_triangular(int width, int height)
+grid *grid_new_triangular(int width, int height, char *desc)
 {
     int x,y;
     
     /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int vec_x = 15;
-    int vec_y = 26;
+    int vec_x = TRIANGLE_VEC_X;
+    int vec_y = TRIANGLE_VEC_Y;
     
     int index;
 
     /* convenient alias */
     int w = width + 1;
 
-    grid *g = grid_new();
-    g->tilesize = 18; /* adjust to your taste */
+    grid *g = grid_empty();
+    g->tilesize = TRIANGLE_TILESIZE;
 
     g->num_faces = width * height * 2;
     g->num_dots = (width + 1) * (height + 1);
@@ -1317,21 +1588,36 @@ grid *grid_new_triangular(int width, int height)
     return g;
 }
 
-grid *grid_new_snubsquare(int width, int height)
+#define SNUBSQUARE_TILESIZE 18
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define SNUBSQUARE_A 15
+#define SNUBSQUARE_B 26
+
+void grid_size_snubsquare(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = SNUBSQUARE_A;
+    int b = SNUBSQUARE_B;
+
+    *tilesize = SNUBSQUARE_TILESIZE;
+    *xextent = (a+b) * (width-1) + a + b;
+    *yextent = (a+b) * (height-1) + a + b;
+}
+
+grid *grid_new_snubsquare(int width, int height, char *desc)
 {
     int x, y;
-    /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = SNUBSQUARE_A;
+    int b = SNUBSQUARE_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 3 * width * height;
     int max_dots = 2 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 18;
+    grid *g = grid_empty();
+    g->tilesize = SNUBSQUARE_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1416,21 +1702,35 @@ grid *grid_new_snubsquare(int width, int height)
     return g;
 }
 
-grid *grid_new_cairo(int width, int height)
+#define CAIRO_TILESIZE 40
+/* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
+#define CAIRO_A 14
+#define CAIRO_B 31
+
+void grid_size_cairo(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int b = CAIRO_B; /* a unused in determining grid size. */
+
+    *tilesize = CAIRO_TILESIZE;
+    *xextent = 2*b*(width-1) + 2*b;
+    *yextent = 2*b*(height-1) + 2*b;
+}
+
+grid *grid_new_cairo(int width, int height, char *desc)
 {
     int x, y;
-    /* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
-    int a = 14;
-    int b = 31;
+    int a = CAIRO_A;
+    int b = CAIRO_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 2 * width * height;
     int max_dots = 3 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 40;
+    grid *g = grid_empty();
+    g->tilesize = CAIRO_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1508,12 +1808,27 @@ grid *grid_new_cairo(int width, int height)
     return g;
 }
 
-grid *grid_new_greathexagonal(int width, int height)
+#define GREATHEX_TILESIZE 18
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define GREATHEX_A 15
+#define GREATHEX_B 26
+
+void grid_size_greathexagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = GREATHEX_A;
+    int b = GREATHEX_B;
+
+    *tilesize = GREATHEX_TILESIZE;
+    *xextent = (3*a + b) * (width-1) + 4*a;
+    *yextent = (2*a + 2*b) * (height-1) + 3*b + a;
+}
+
+grid *grid_new_greathexagonal(int width, int height, char *desc)
 {
     int x, y;
-    /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = GREATHEX_A;
+    int b = GREATHEX_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 6 * (width + 1) * (height + 1);
@@ -1521,8 +1836,8 @@ grid *grid_new_greathexagonal(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 18;
+    grid *g = grid_empty();
+    g->tilesize = GREATHEX_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1623,12 +1938,27 @@ grid *grid_new_greathexagonal(int width, int height)
     return g;
 }
 
-grid *grid_new_octagonal(int width, int height)
+#define OCTAGONAL_TILESIZE 40
+/* b/a approx sqrt(2) */
+#define OCTAGONAL_A 29
+#define OCTAGONAL_B 41
+
+void grid_size_octagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = OCTAGONAL_A;
+    int b = OCTAGONAL_B;
+
+    *tilesize = OCTAGONAL_TILESIZE;
+    *xextent = (2*a + b) * width;
+    *yextent = (2*a + b) * height;
+}
+
+grid *grid_new_octagonal(int width, int height, char *desc)
 {
     int x, y;
-    /* b/a approx sqrt(2) */
-    int a = 29;
-    int b = 41;
+    int a = OCTAGONAL_A;
+    int b = OCTAGONAL_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 2 * width * height;
@@ -1636,8 +1966,8 @@ grid *grid_new_octagonal(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 40;
+    grid *g = grid_empty();
+    g->tilesize = OCTAGONAL_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1691,12 +2021,27 @@ grid *grid_new_octagonal(int width, int height)
     return g;
 }
 
-grid *grid_new_kites(int width, int height)
+#define KITE_TILESIZE 40
+/* b/a approx sqrt(3) */
+#define KITE_A 15
+#define KITE_B 26
+
+void grid_size_kites(int width, int height,
+                     int *tilesize, int *xextent, int *yextent)
+{
+    int a = KITE_A;
+    int b = KITE_B;
+
+    *tilesize = KITE_TILESIZE;
+    *xextent = 4*b * width + 2*b;
+    *yextent = 6*a * (height-1) + 8*a;
+}
+
+grid *grid_new_kites(int width, int height, char *desc)
 {
     int x, y;
-    /* b/a approx sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = KITE_A;
+    int b = KITE_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 6 * width * height;
@@ -1704,8 +2049,8 @@ grid *grid_new_kites(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 40;
+    grid *g = grid_empty();
+    g->tilesize = KITE_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1796,25 +2141,45 @@ grid *grid_new_kites(int width, int height)
     return g;
 }
 
-grid *grid_new_floret(int width, int height)
+#define FLORET_TILESIZE 150
+/* -py/px is close to tan(30 - atan(sqrt(3)/9))
+ * using py=26 makes everything lean to the left, rather than right
+ */
+#define FLORET_PX 75
+#define FLORET_PY -26
+
+void grid_size_floret(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
+    int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
+    int ry = qy-py;
+    /* rx unused in determining grid size. */
+
+    *tilesize = FLORET_TILESIZE;
+    *xextent = (6*px+3*qx)/2 * (width-1) + 4*qx + 2*px;
+    *yextent = (5*qy-4*py) * (height-1) + 4*qy + 2*ry;
+}
+
+grid *grid_new_floret(int width, int height, char *desc)
 {
     int x, y;
     /* Vectors for sides; weird numbers needed to keep puzzle aligned with window
      * -py/px is close to tan(30 - atan(sqrt(3)/9))
      * using py=26 makes everything lean to the left, rather than right
      */
-    int px = 75, py = -26;       /* |( 75, -26)| = 79.43 */
-    int qx = 4*px/5, qy = -py*2; /* |( 60,  52)| = 79.40 */
-    int rx = qx-px, ry = qy-py;  /* |(-15,  78)| = 79.38 */
+    int px = FLORET_PX, py = FLORET_PY;         /* |( 75, -26)| = 79.43 */
+    int qx = 4*px/5, qy = -py*2;                /* |( 60,  52)| = 79.40 */
+    int rx = qx-px, ry = qy-py;                 /* |(-15,  78)| = 79.38 */
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 6 * width * height;
     int max_dots = 9 * (width + 1) * (height + 1);
-    
+
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = 2 * px;
+    grid *g = grid_empty();
+    g->tilesize = FLORET_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1884,12 +2249,28 @@ grid *grid_new_floret(int width, int height)
     return g;
 }
 
-grid *grid_new_dodecagonal(int width, int height)
+/* DODEC_* are used for dodecagonal and great-dodecagonal grids. */
+#define DODEC_TILESIZE 26
+/* Vector for side of triangle - ratio is close to sqrt(3) */
+#define DODEC_A 15
+#define DODEC_B 26
+
+void grid_size_dodecagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = DODEC_A;
+    int b = DODEC_B;
+
+    *tilesize = DODEC_TILESIZE;
+    *xextent = (4*a + 2*b) * (width-1) + 3*(2*a + b);
+    *yextent = (3*a + 2*b) * (height-1) + 2*(2*a + b);
+}
+
+grid *grid_new_dodecagonal(int width, int height, char *desc)
 {
     int x, y;
-    /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = DODEC_A;
+    int b = DODEC_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 3 * width * height;
@@ -1897,8 +2278,8 @@ grid *grid_new_dodecagonal(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = b;
+    grid *g = grid_empty();
+    g->tilesize = DODEC_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -1954,12 +2335,23 @@ grid *grid_new_dodecagonal(int width, int height)
     return g;
 }
 
-grid *grid_new_greatdodecagonal(int width, int height)
+void grid_size_greatdodecagonal(int width, int height,
+                          int *tilesize, int *xextent, int *yextent)
+{
+    int a = DODEC_A;
+    int b = DODEC_B;
+
+    *tilesize = DODEC_TILESIZE;
+    *xextent = (6*a + 2*b) * (width-1) + 2*(2*a + b) + 3*a + b;
+    *yextent = (3*a + 3*b) * (height-1) + 2*(2*a + b);
+}
+
+grid *grid_new_greatdodecagonal(int width, int height, char *desc)
 {
     int x, y;
     /* Vector for side of triangle - ratio is close to sqrt(3) */
-    int a = 15;
-    int b = 26;
+    int a = DODEC_A;
+    int b = DODEC_B;
 
     /* Upper bounds - don't have to be exact */
     int max_faces = 30 * width * height;
@@ -1967,8 +2359,8 @@ grid *grid_new_greatdodecagonal(int width, int height)
 
     tree234 *points;
 
-    grid *g = grid_new();
-    g->tilesize = b;
+    grid *g = grid_empty();
+    g->tilesize = DODEC_TILESIZE;
     g->faces = snewn(max_faces, grid_face);
     g->dots = snewn(max_dots, grid_dot);
 
@@ -2057,4 +2449,279 @@ grid *grid_new_greatdodecagonal(int width, int height)
     return g;
 }
 
+typedef struct setface_ctx
+{
+    int xmin, xmax, ymin, ymax;
+    int aoff;
+
+    grid *g;
+    tree234 *points;
+} setface_ctx;
+
+double round(double r)
+{
+    return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
+}
+
+int set_faces(penrose_state *state, vector *vs, int n, int depth)
+{
+    setface_ctx *sf_ctx = (setface_ctx *)state->ctx;
+    int i;
+    int xs[4], ys[4];
+    double cosa = cos(sf_ctx->aoff * PI / 180.0);
+    double sina = sin(sf_ctx->aoff * PI / 180.0);
+
+    if (depth < state->max_depth) return 0;
+#ifdef DEBUG_PENROSE
+    if (n != 4) return 0; /* triangles are sent as debugging. */
+#endif
+
+    for (i = 0; i < n; i++) {
+        double tx = v_x(vs, i), ty = v_y(vs, i);
+
+        xs[i] = (int)round( tx*cosa + ty*sina);
+        ys[i] = (int)round(-tx*sina + ty*cosa);
+
+        if (xs[i] < sf_ctx->xmin || xs[i] > sf_ctx->xmax) return 0;
+        if (ys[i] < sf_ctx->ymin || ys[i] > sf_ctx->ymax) return 0;
+    }
+
+    grid_face_add_new(sf_ctx->g, n);
+    debug(("penrose: new face l=%f gen=%d...",
+           penrose_side_length(state->start_size, depth), depth));
+    for (i = 0; i < n; i++) {
+        grid_dot *d = grid_get_dot(sf_ctx->g, sf_ctx->points,
+                                   xs[i], ys[i]);
+        grid_face_set_dot(sf_ctx->g, d, i);
+        debug((" ... dot 0x%x (%d,%d) (was %2.2f,%2.2f)",
+               d, d->x, d->y, v_x(vs, i), v_y(vs, i)));
+    }
+
+    return 0;
+}
+
+#define PENROSE_TILESIZE 100
+
+void grid_size_penrose(int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    int l = PENROSE_TILESIZE;
+
+    *tilesize = l;
+    *xextent = l * width;
+    *yextent = l * height;
+}
+
+static char *grid_new_desc_penrose(grid_type type, int width, int height, random_state *rs)
+{
+    int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff;
+    double outer_radius;
+    int inner_radius;
+    char gd[255];
+    int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
+
+    /* We want to produce a random bit of penrose tiling, so we calculate
+     * a random offset (within the patch that penrose.c calculates for us)
+     * and an angle (multiple of 36) to rotate the patch. */
+
+    penrose_calculate_size(which, tilesize, width, height,
+                           &outer_radius, &startsz, &depth);
+
+    /* Calculate radius of (circumcircle of) patch, subtract from
+     * radius calculated. */
+    inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
+
+    /* Pick a random offset (the easy way: choose within outer square,
+     * discarding while it's outside the circle) */
+    do {
+        xoff = random_upto(rs, 2*inner_radius) - inner_radius;
+        yoff = random_upto(rs, 2*inner_radius) - inner_radius;
+    } while (sqrt(xoff*xoff+yoff*yoff) > inner_radius);
+
+    aoff = random_upto(rs, 360/36) * 36;
+
+    debug(("grid_desc: ts %d, %dx%d patch, orad %2.2f irad %d",
+           tilesize, width, height, outer_radius, inner_radius));
+    debug(("    -> xoff %d yoff %d aoff %d", xoff, yoff, aoff));
+
+    sprintf(gd, "G%d,%d,%d", xoff, yoff, aoff);
+
+    return dupstr(gd);
+}
+
+static char *grid_validate_desc_penrose(grid_type type, int width, int height, char *desc)
+{
+    int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff, inner_radius;
+    double outer_radius;
+    int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
+
+    if (!desc)
+        return "Missing grid description string.";
+
+    penrose_calculate_size(which, tilesize, width, height,
+                           &outer_radius, &startsz, &depth);
+    inner_radius = (int)(outer_radius - sqrt(width*width + height*height));
+
+    if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &aoff) != 3)
+        return "Invalid format grid description string.";
+
+    if (sqrt(xoff*xoff + yoff*yoff) > inner_radius)
+        return "Patch offset out of bounds.";
+    if ((aoff % 36) != 0 || aoff < 0 || aoff >= 360)
+        return "Angle offset out of bounds.";
+
+    return NULL;
+}
+
+/*
+ * We're asked for a grid of a particular size, and we generate enough
+ * of the tiling so we can be sure to have enough random grid from which
+ * to pick.
+ */
+
+static grid *grid_new_penrose(int width, int height, int which, char *desc)
+{
+    int max_faces, max_dots, tilesize = PENROSE_TILESIZE;
+    int xsz, ysz, xoff, yoff;
+    double rradius;
+
+    tree234 *points;
+    grid *g;
+
+    penrose_state ps;
+    setface_ctx sf_ctx;
+
+    penrose_calculate_size(which, tilesize, width, height,
+                           &rradius, &ps.start_size, &ps.max_depth);
+
+    debug(("penrose: w%d h%d, tile size %d, start size %d, depth %d",
+           width, height, tilesize, ps.start_size, ps.max_depth));
+
+    ps.new_tile = set_faces;
+    ps.ctx = &sf_ctx;
+
+    max_faces = (width*3) * (height*3); /* somewhat paranoid... */
+    max_dots = max_faces * 4; /* ditto... */
+
+    g = grid_empty();
+    g->tilesize = tilesize;
+    g->faces = snewn(max_faces, grid_face);
+    g->dots = snewn(max_dots, grid_dot);
+
+    points = newtree234(grid_point_cmp_fn);
+
+    memset(&sf_ctx, 0, sizeof(sf_ctx));
+    sf_ctx.g = g;
+    sf_ctx.points = points;
+
+    if (desc != NULL) {
+        if (sscanf(desc, "G%d,%d,%d", &xoff, &yoff, &sf_ctx.aoff) != 3)
+            assert(!"Invalid grid description.");
+    } else {
+        xoff = yoff = 0;
+    }
+
+    xsz = width * tilesize;
+    ysz = height * tilesize;
+
+    sf_ctx.xmin = xoff - xsz/2;
+    sf_ctx.xmax = xoff + xsz/2;
+    sf_ctx.ymin = yoff - ysz/2;
+    sf_ctx.ymax = yoff + ysz/2;
+
+    debug(("penrose: centre (%f, %f) xsz %f ysz %f",
+           0.0, 0.0, xsz, ysz));
+    debug(("penrose: x range (%f --> %f), y range (%f --> %f)",
+           sf_ctx.xmin, sf_ctx.xmax, sf_ctx.ymin, sf_ctx.ymax));
+
+    penrose(&ps, which);
+
+    freetree234(points);
+    assert(g->num_faces <= max_faces);
+    assert(g->num_dots <= max_dots);
+
+    debug(("penrose: %d faces total (equivalent to %d wide by %d high)",
+           g->num_faces, g->num_faces/height, g->num_faces/width));
+
+    grid_trim_vigorously(g);
+    grid_make_consistent(g);
+
+    /*
+     * Centre the grid in its originally promised rectangle.
+     */
+    g->lowest_x -= ((sf_ctx.xmax - sf_ctx.xmin) -
+                    (g->highest_x - g->lowest_x)) / 2;
+    g->highest_x = g->lowest_x + (sf_ctx.xmax - sf_ctx.xmin);
+    g->lowest_y -= ((sf_ctx.ymax - sf_ctx.ymin) -
+                    (g->highest_y - g->lowest_y)) / 2;
+    g->highest_y = g->lowest_y + (sf_ctx.ymax - sf_ctx.ymin);
+
+    return g;
+}
+
+void grid_size_penrose_p2_kite(int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    grid_size_penrose(width, height, tilesize, xextent, yextent);
+}
+
+void grid_size_penrose_p3_thick(int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    grid_size_penrose(width, height, tilesize, xextent, yextent);
+}
+
+grid *grid_new_penrose_p2_kite(int width, int height, char *desc)
+{
+    return grid_new_penrose(width, height, PENROSE_P2, desc);
+}
+
+grid *grid_new_penrose_p3_thick(int width, int height, char *desc)
+{
+    return grid_new_penrose(width, height, PENROSE_P3, desc);
+}
+
 /* ----------- End of grid generators ------------- */
+
+#define FNNEW(upper,lower) &grid_new_ ## lower,
+#define FNSZ(upper,lower) &grid_size_ ## lower,
+
+static grid *(*(grid_news[]))(int, int, char*) = { GRIDGEN_LIST(FNNEW) };
+static void(*(grid_sizes[]))(int, int, int*, int*, int*) = { GRIDGEN_LIST(FNSZ) };
+
+char *grid_new_desc(grid_type type, int width, int height, random_state *rs)
+{
+    if (type != GRID_PENROSE_P2 && type != GRID_PENROSE_P3)
+        return NULL;
+
+    return grid_new_desc_penrose(type, width, height, rs);
+}
+
+char *grid_validate_desc(grid_type type, int width, int height, char *desc)
+{
+    if (type != GRID_PENROSE_P2 && type != GRID_PENROSE_P3) {
+        if (desc != NULL)
+            return "Grid description strings not used with this grid type";
+        return NULL;
+    }
+
+    return grid_validate_desc_penrose(type, width, height, desc);
+}
+
+grid *grid_new(grid_type type, int width, int height, char *desc)
+{
+    char *err = grid_validate_desc(type, width, height, desc);
+    if (err) assert(!"Invalid grid description.");
+
+    return grid_news[type](width, height, desc);
+}
+
+void grid_compute_size(grid_type type, int width, int height,
+                       int *tilesize, int *xextent, int *yextent)
+{
+    grid_sizes[type](width, height, tilesize, xextent, yextent);
+}
+
+/* ----------- End of grid helpers ------------- */
+
+/* vim: set shiftwidth=4 tabstop=8: */
diff --git a/grid.h b/grid.h
index 65ced86..d1c260e 100644 (file)
--- a/grid.h
+++ b/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);
-grid *grid_new_honeycomb(int width, int height);
-grid *grid_new_triangular(int width, int height);
-grid *grid_new_snubsquare(int width, int height);
-grid *grid_new_cairo(int width, int height);
-grid *grid_new_greathexagonal(int width, int height);
-grid *grid_new_octagonal(int width, int height);
-grid *grid_new_kites(int width, int height);
-grid *grid_new_floret(int width, int height);
-grid *grid_new_dodecagonal(int width, int height);
-grid *grid_new_greatdodecagonal(int width, int height);
+/* Grids are specified by type: GRID_SQUARE, GRID_KITE, etc. */
+
+#define GRIDGEN_LIST(A) \
+  A(SQUARE,square) \
+  A(HONEYCOMB,honeycomb) \
+  A(TRIANGULAR,triangular) \
+  A(SNUBSQUARE,snubsquare) \
+  A(CAIRO,cairo) \
+  A(GREATHEXAGONAL,greathexagonal) \
+  A(OCTAGONAL,octagonal) \
+  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 */
diff --git a/loopy.R b/loopy.R
index e77f9c7..c15e39c 100644 (file)
--- a/loopy.R
+++ b/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
diff --git a/loopy.c b/loopy.c
index 1b32635..6e97acd 100644 (file)
--- a/loopy.c
+++ b/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(Triangular,grid_new_triangular,3,3) \
-    A(Honeycomb,grid_new_honeycomb,3,3) \
-    A(Snub-Square,grid_new_snubsquare,3,3) \
-    A(Cairo,grid_new_cairo,3,4) \
-    A(Great-Hexagonal,grid_new_greathexagonal,3,3) \
-    A(Octagonal,grid_new_octagonal,3,3) \
-    A(Kites,grid_new_kites,3,3) \
-    A(Floret,grid_new_floret,1,2) \
-    A(Dodecagonal,grid_new_dodecagonal,2,2) \
-    A(Great-Dodecagonal,grid_new_greatdodecagonal,2,2)
-
-#define GRID_NAME(title,fn,amin,omin) #title,
-#define GRID_CONFIG(title,fn,amin,omin) ":" #title
-#define GRID_FN(title,fn,amin,omin) &fn,
-#define GRID_SIZES(title,fn,amin,omin) \
+    A(Squares,GRID_SQUARE,3,3) \
+    A(Triangular,GRID_TRIANGULAR,3,3) \
+    A(Honeycomb,GRID_HONEYCOMB,3,3) \
+    A(Snub-Square,GRID_SNUBSQUARE,3,3) \
+    A(Cairo,GRID_CAIRO,3,4) \
+    A(Great-Hexagonal,GRID_GREATHEXAGONAL,3,3) \
+    A(Octagonal,GRID_OCTAGONAL,3,3) \
+    A(Kites,GRID_KITE,3,3) \
+    A(Floret,GRID_FLORET,1,2) \
+    A(Dodecagonal,GRID_DODECAGONAL,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,type,amin,omin) #title,
+#define GRID_CONFIG(title,type,amin,omin) ":" #title
+#define GRID_TYPE(title,type,amin,omin) type,
+#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) };
-#define NUM_GRID_TYPES (sizeof(grid_fns) / sizeof(grid_fns[0]))
+static grid_type grid_types[] = { GRIDLIST(GRID_TYPE) };
+#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
- * freed in free_params(). */
-static void params_generate_grid(game_params *params)
+ * generated. */
+static grid *loopy_generate_grid(game_params *params, char *grid_desc)
 {
-    if (!params->game_grid) {
-        params->game_grid = grid_fns[params->type](params->w, params->h);
-    }
+    return grid_new(grid_types[params->type], params->w, params->h, grid_desc);
 }
 
 /* ----------------------------------------------------------------------
@@ -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_NORMAL, 0, NULL },
-    {  7,  7, DIFF_HARD, 0, NULL },
-    {  7,  7, DIFF_HARD, 1, NULL },
-    {  7,  7, DIFF_HARD, 2, NULL },
-    {  5,  5, DIFF_HARD, 3, NULL },
-    {  7,  7, DIFF_HARD, 4, NULL },
-    {  5,  4, DIFF_HARD, 5, NULL },
-    {  5,  5, DIFF_HARD, 6, NULL },
-    {  5,  5, DIFF_HARD, 7, NULL },
-    {  3,  3, DIFF_HARD, 8, NULL },
-    {  3,  3, DIFF_HARD, 9, NULL },
-    {  3,  3, DIFF_HARD, 10, NULL },
+    {  7,  7, DIFF_EASY, 0 },
+    {  7,  7, DIFF_NORMAL, 0 },
+    {  7,  7, DIFF_HARD, 0 },
+    {  7,  7, DIFF_HARD, 1 },
+    {  7,  7, DIFF_HARD, 2 },
+    {  5,  5, DIFF_HARD, 3 },
+    {  7,  7, DIFF_HARD, 4 },
+    {  5,  4, DIFF_HARD, 5 },
+    {  5,  5, DIFF_HARD, 6 },
+    {  5,  5, DIFF_HARD, 7 },
+    {  3,  3, DIFF_HARD, 8 },
+    {  3,  3, DIFF_HARD, 9 },
+    {  3,  3, DIFF_HARD, 10 },
+    {  6,  6, DIFF_HARD, 11 },
+    {  6,  6, DIFF_HARD, 12 },
 #else
-    {  7,  7, DIFF_EASY, 0, NULL },
-    {  10,  10, DIFF_EASY, 0, NULL },
-    {  7,  7, DIFF_NORMAL, 0, NULL },
-    {  10,  10, DIFF_NORMAL, 0, NULL },
-    {  7,  7, DIFF_HARD, 0, NULL },
-    {  10,  10, DIFF_HARD, 0, NULL },
-    {  10,  10, DIFF_HARD, 1, NULL },
-    {  12,  10, DIFF_HARD, 2, NULL },
-    {  7,  7, DIFF_HARD, 3, NULL },
-    {  9,  9, DIFF_HARD, 4, NULL },
-    {  5,  4, DIFF_HARD, 5, NULL },
-    {  7,  7, DIFF_HARD, 6, NULL },
-    {  5,  5, DIFF_HARD, 7, NULL },
-    {  5,  5, DIFF_HARD, 8, NULL },
-    {  5,  4, DIFF_HARD, 9, NULL },
-    {  5,  4, DIFF_HARD, 10, NULL },
+    {  7,  7, DIFF_EASY, 0 },
+    {  10,  10, DIFF_EASY, 0 },
+    {  7,  7, DIFF_NORMAL, 0 },
+    {  10,  10, DIFF_NORMAL, 0 },
+    {  7,  7, DIFF_HARD, 0 },
+    {  10,  10, DIFF_HARD, 0 },
+    {  10,  10, DIFF_HARD, 1 },
+    {  12,  10, DIFF_HARD, 2 },
+    {  7,  7, DIFF_HARD, 3 },
+    {  9,  9, DIFF_HARD, 4 },
+    {  5,  4, DIFF_HARD, 5 },
+    {  7,  7, DIFF_HARD, 6 },
+    {  5,  5, DIFF_HARD, 7 },
+    {  5,  5, DIFF_HARD, 8 },
+    {  5,  4, DIFF_HARD, 9 },
+    {  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);
-    g = params->game_grid;
+    char *grid_desc, *ret;
+
+    /* 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_height = grid_height * tilesize / g->tilesize;
+    rendered_width = grid_width * 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;
-    g->refcount++;
+
+    grid_desc = grid_new_desc(grid_types[params->type], params->w, params->h, rs);
+    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);
-    state->game_grid = g = params->game_grid;
-    g->refcount++;
+    grid_desc = extract_grid_desc(&desc);
+    state->game_grid = g = loopy_generate_grid(params, grid_desc);
+    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: */
diff --git a/penrose.c b/penrose.c
new file mode 100644 (file)
index 0000000..18c4258
--- /dev/null
+++ b/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: */
diff --git a/penrose.h b/penrose.h
new file mode 100644 (file)
index 0000000..f7b23e2
--- /dev/null
+++ b/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