ref: 61e9c782487ea982498955b93d1b94921131059e
parent: 23d4322fec38f8ffac930bc541e08309e1d02f11
author: Simon Tatham <[email protected]>
date: Fri Jul 7 14:16:04 EDT 2023
grid.c: new and improved Penrose tiling generator. The new generator works on the same 'combinatorial coordinates' system as the more recently written Hats and Spectre generators. When I came up with that technique for the Hats tiling, I was already tempted to rewrite the Penrose generator on the same principle, because it has a lot of advantages over the previous technique of picking a randomly selected patch out of the subdivision of a huge starting tile. It generates the exact limiting distribution of possible tiling patches rather than an approximation based on how big a tile you subdivided; it doesn't use any overly large integers or overly precise floating point to suffer overflow or significance loss, because it never has to even _think_ about the coordinates of any point not in the actual output area. But I resisted the temptation to throw out our existing Penrose generator and move to the shiny new algorithm, for just one reason: backwards compatibility. There's no sensible way to translate existing Loopy game IDs into the new notation, so they would stop working, unless we kept the old generator around as well to interpret the previous system. And although _random seeds_ aren't guaranteed to generate the same result from one build of Puzzles to the next, I do try to keep existing descriptive game IDs working. So if we had to keep the old generator around anyway, it didn't seem worth writing a new one... ... at least, until we discovered that the old generator has a latent bug. The function that decides whether each tile is within the target area, and hence whether to make it part of the output grid, is based on floating-point calculation of the tile's vertices. So a change in FP rounding behaviour between two platforms could cause them to interpret the same grid description differently, invalidating a Loopy game on that grid. This is _unlikely_, as long as everyone uses IEEE 754 double, but the C standard doesn't actually guarantee that. We found this out when I started investigating a slight distortion in large instances of Penrose Loopy. For example, the Loopy random seed "40x40t11#12345", as of just before this commit, generates a game description beginning with the Penrose grid string "G-4944,5110,108", in which you can see a star shape of five darts a few tiles down the left edge, where two of the radii of the star don't properly line up with edges in the surrounding shell of kites that they should be collinear with. This turns out to be due to FP error: not because _double precision_ ran out, but because the definitions of COS54, SIN54, COS18 and SIN18 in penrose.c were specified to only 7 digits of precision. And when you expand a patch of tiling that large, you end up with integer combinations of those numbers with coefficients about 7 digits long, mostly cancelling out to leave a much smaller output value, and the inaccuracies in those constants start to be noticeable. But having noticed that, it then became clear that these badly computed values were also critical to _correctness_ of the grid. So they can't be fixed without invalidating old game IDs. _And_ then we realised that this also means existing platforms might disagree on a game ID's validity. So if we have to break compatibility anyway, we should go all the way, and instead of fixing the old system, introduce the shiny new system that gets all of this right. Therefore, the new penrose.c uses the more reliable combinatorial-coordinates system which doesn't deal in integers that large in the first place. Also, there's no longer any floating point at all in the calculation of tile vertex coordinates: the combinations of 1 and sqrt(5) are computed exactly by the n_times_root_k function. So now a large Penrose grid should never have obvious failures of alignment like that. The old system is kept for backwards compatibility. It's not fully reliable, because that bug hasn't been fixed - but any Penrose Loopy game ID that was working before on _some_ platform should still work on that one. However, new Penrose Loopy game IDs are based on combinatorial coordinates, computed in exact arithmetic, and should be properly stable. The new code looks suspiciously like the Spectre code (though considerably simpler - the Penrose coordinate maps are easy enough that I just hand-typed them rather than having an elaborate auxiliary data-generation tool). That's because they were similar enough in concept to make it possible to clone and hack. But there are enough divergences that it didn't seem like a good idea to try to merge them again afterwards - in particular, the fact that output Penrose tiles are generated by merging two triangular metatiles while Spectres are subdivisions of their metatiles.
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,7 +8,8 @@
add_library(core_obj OBJECT
combi.c divvy.c drawing.c dsf.c findloop.c grid.c latin.c
laydomino.c loopgen.c malloc.c matching.c midend.c misc.c penrose.c
- ps.c random.c sort.c tdq.c tree234.c version.c ${platform_common_sources})+ penrose-legacy.c ps.c random.c sort.c tdq.c tree234.c version.c
+ ${platform_common_sources})add_library(core $<TARGET_OBJECTS:core_obj>)
add_library(common $<TARGET_OBJECTS:core_obj> hat.c spectre.c)
--- a/auxiliary/CMakeLists.txt
+++ b/auxiliary/CMakeLists.txt
@@ -5,6 +5,7 @@
cliprogram(latin-test latin-test.c)
cliprogram(matching matching.c)
cliprogram(obfusc obfusc.c)
+cliprogram(penrose-legacy-test penrose-legacy-test.c)
cliprogram(penrose-test penrose-test.c)
cliprogram(sort-test sort-test.c)
cliprogram(spectre-gen spectre-gen.c spectre-help.c CORE_LIB)
--- /dev/null
+++ b/auxiliary/penrose-legacy-test.c
@@ -1,0 +1,98 @@
+#include <stdio.h>
+#include <string.h>
+
+#include "puzzles.h"
+#include "penrose-legacy.h"
+
+static int show_recursion = 0;
+static int ntiles, nfinal;
+
+static int test_cb(penrose_legacy_state *state, vector *vs, int n, int depth)
+{+ int i, xoff = 0, yoff = 0;
+ double l = penrose_legacy_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) ? "" : " ",+ penrose_legacy_vx(vs, i) + xoff,
+ penrose_legacy_vy(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\" />",+ penrose_legacy_vx(vs, 0) + xoff, penrose_legacy_vy(vs, 0) + yoff,
+ rball, rball, col);
+
+ return 0;
+}
+
+static void usage_exit(void)
+{+ fprintf(stderr, "Usage: penrose-legacy-test [--recursion] "
+ "P2|P3 SIZE DEPTH\n");
+ exit(1);
+}
+
+int main(int argc, char *argv[])
+{+ penrose_legacy_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_legacy(&ps, which, 0);
+ printf("</g>\n");+
+ printf("<!-- %d tiles and %d leaf tiles total -->\n",+ ntiles, nfinal);
+
+ printf("</svg>");+
+ return 0;
+}
--- a/auxiliary/penrose-test.c
+++ b/auxiliary/penrose-test.c
@@ -1,95 +1,46 @@
#include <stdio.h>
-#include <string.h>
+#include <math.h>
#include "puzzles.h"
#include "penrose.h"
+#include "penrose-internal.h"
-static int show_recursion = 0;
-static int ntiles, nfinal;
+struct testctx {+ double sqrt5, xunit, yunit;
+};
-static int test_cb(penrose_state *state, vector *vs, int n, int depth)
+static void tile(void *vctx, const int *coords)
{- int i, xoff = 0, yoff = 0;
- double l = penrose_side_length(state->start_size, depth);
- double rball = l / 10.0;
- const char *col;
+ struct testctx *tctx = (struct testctx *)vctx;
+ size_t i;
- ntiles++;
- if (state->max_depth == depth) {- col = n == 4 ? "black" : "green";
- nfinal++;
- } else {- if (!show_recursion)
- return 0;
- col = n == 4 ? "red" : "blue";
+ printf("newpath");+ for (i = 0; i < 4; i++) {+ printf(" %g %g %s",+ tctx->xunit * (coords[4*i+0] + tctx->sqrt5 * coords[4*i+1]),
+ tctx->yunit * (coords[4*i+2] + tctx->sqrt5 * coords[4*i+3]),
+ i == 0 ? "moveto" : "lineto");
}
- 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;
+ printf(" closepath gsave 0.7 setgray fill grestore stroke\n");}
-static void usage_exit(void)
+int main(void)
{- fprintf(stderr, "Usage: penrose-test [--recursion] P2|P3 SIZE DEPTH\n");
- exit(1);
-}
+ random_state *rs;
+ struct PenrosePatchParams params[1];
+ struct testctx tctx[1];
+ int w = 50, h = 40;
-int main(int argc, char *argv[])
-{- penrose_state ps;
- int which = 0;
+ tctx->sqrt5 = sqrt(5);
+ tctx->xunit = 25 * 0.25;
+ tctx->yunit = 25 * sin(atan2(0, -1)/5) / 2;
+ printf("newpath 0 0 moveto %g 0 rlineto 0 %g rlineto %g 0 rlineto "+ "closepath stroke\n", w * tctx->xunit, h * tctx->yunit,
+ -w * tctx->xunit);
- 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, 0);
- printf("</g>\n");-
- printf("<!-- %d tiles and %d leaf tiles total -->\n",- ntiles, nfinal);
-
- printf("</svg>");-
- return 0;
+ rs = random_new("12345", 5);+ penrose_tiling_randomise(params, PENROSE_P2, w, h, rs);
+ penrose_tiling_generate(params, w, h, tile, tctx);
+ sfree(params->coords);
+ random_free(rs);
}
--- a/grid.c
+++ b/grid.c
@@ -22,6 +22,7 @@
#include "puzzles.h"
#include "tree234.h"
#include "grid.h"
+#include "penrose-legacy.h"
#include "penrose.h"
#include "hat.h"
#include "spectre.h"
@@ -3024,13 +3025,48 @@
return g;
}
-typedef struct setface_ctx
+/*
+ * Penrose tilings. For historical reasons, we support two totally
+ * different generation algorithms: the legacy one is only supported
+ * by grid_new_penrose, for backwards compatibility with game
+ * descriptions generated before we switched. New grid generation uses
+ * only the new system.
+ */
+
+#define PENROSE_TILESIZE 100
+
+static const char *grid_validate_params_penrose(int width, int height)
{+ int l = PENROSE_TILESIZE;
+
+ if (width > INT_MAX / l || /* xextent */
+ height > INT_MAX / l || /* yextent */
+ width > INT_MAX / (3 * 3 * 4 * height)) /* max_dots */
+ return "Grid must not be unreasonably large";
+ return NULL;
+}
+
+static 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;
+}
+
+/*
+ * Legacy generation by selecting a patch of tiling from the expansion
+ * of a big triangle.
+ */
+
+typedef struct penrose_legacy_set_faces_ctx {int xmin, xmax, ymin, ymax;
grid *g;
tree234 *points;
-} setface_ctx;
+} penrose_legacy_set_faces_ctx;
static double round_int_nearest_away(double r)
{@@ -3037,9 +3073,11 @@
return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
}
-static int set_faces(penrose_state *state, vector *vs, int n, int depth)
+static int penrose_legacy_set_faces(penrose_legacy_state *state, vector *vs,
+ int n, int depth)
{- setface_ctx *sf_ctx = (setface_ctx *)state->ctx;
+ penrose_legacy_set_faces_ctx *sf_ctx =
+ (penrose_legacy_set_faces_ctx *)state->ctx;
int i;
int xs[4], ys[4];
@@ -3049,7 +3087,7 @@
#endif
for (i = 0; i < n; i++) {- double tx = v_x(vs, i), ty = v_y(vs, i);
+ double tx = penrose_legacy_vx(vs, i), ty = penrose_legacy_vy(vs, i);
xs[i] = (int)round_int_nearest_away(tx);
ys[i] = (int)round_int_nearest_away(ty);
@@ -3060,100 +3098,25 @@
grid_face_add_new(sf_ctx->g, n);
debug(("penrose: new face l=%f gen=%d...",- penrose_side_length(state->start_size, depth), depth));
+ penrose_legacy_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)));
+ d, d->x, d->y, penrose_legacy_vx(vs, i),
+ penrose_legacy_vy(vs, i)));
}
return 0;
}
-#define PENROSE_TILESIZE 100
+static grid *grid_new_penrose_legacy(int width, int height, int which,
+ const char *desc);
-static const char *grid_validate_params_penrose(int width, int height)
+static const char *grid_validate_desc_penrose_legacy(
+ grid_type type, int width, int height, const char *desc)
{- int l = PENROSE_TILESIZE;
-
- if (width > INT_MAX / l || /* xextent */
- height > INT_MAX / l || /* yextent */
- width > INT_MAX / (3 * 3 * 4 * height)) /* max_dots */
- return "Grid must not be unreasonably large";
- return NULL;
-}
-
-static 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 grid *grid_new_penrose(int width, int height, int which, const char *desc); /* forward reference */
-
-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);
- grid *g;
-
- while (1) {- /* 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);
-
- /*
- * Now test-generate our grid, to make sure it actually
- * produces something.
- */
- g = grid_new_penrose(width, height, which, gd);
- if (g) {- grid_free(g);
- break;
- }
- /* If not, go back to the top of this while loop and try again
- * with a different random offset. */
- }
-
- return dupstr(gd);
-}
-
-static const char *grid_validate_desc_penrose(grid_type type,
- int width, int height,
- const 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);
@@ -3162,8 +3125,8 @@
if (!desc)
return "Missing grid description string.";
- penrose_calculate_size(which, tilesize, width, height,
- &outer_radius, &startsz, &depth);
+ penrose_legacy_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)
@@ -3178,7 +3141,7 @@
* Test-generate to ensure these parameters don't end us up with
* no grid at all.
*/
- g = grid_new_penrose(width, height, which, desc);
+ g = grid_new_penrose_legacy(width, height, which, desc);
if (!g)
return "Patch coordinates do not identify a usable grid fragment";
grid_free(g);
@@ -3186,13 +3149,8 @@
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, const char *desc)
+static grid *grid_new_penrose_legacy(int width, int height, int which,
+ const char *desc)
{int tilesize = PENROSE_TILESIZE;
int xsz, ysz, xoff, yoff, aoff;
@@ -3201,16 +3159,16 @@
tree234 *points;
grid *g;
- penrose_state ps;
- setface_ctx sf_ctx;
+ penrose_legacy_state ps;
+ penrose_legacy_set_faces_ctx sf_ctx;
- penrose_calculate_size(which, tilesize, width, height,
- &rradius, &ps.start_size, &ps.max_depth);
+ penrose_legacy_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.new_tile = penrose_legacy_set_faces;
ps.ctx = &sf_ctx;
g = grid_empty();
@@ -3242,7 +3200,7 @@
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, aoff);
+ penrose_legacy(&ps, which, aoff);
freetree234(points);
@@ -3278,6 +3236,213 @@
g->highest_y = g->lowest_y + (sf_ctx.ymax - sf_ctx.ymin);
return g;
+}
+
+/*
+ * Combinatorial-coordinate generation.
+ *
+ * We receive coordinates from the generator in the form of x,y pairs
+ * each of which is an integer combination of 1 and sqrt(5), but those
+ * pairs have different scale units in the x and y directions. The
+ * scale units are 1/4 for x and sin(pi/5)/2 for y, which makes their
+ * ratio equal to 2 sin(pi/5) ~= 1.1756. We fudge that irrational
+ * aspect ratio into a rational approximation, by simply taking a pair
+ * of integer scale factors for the x and y dimensions; this distorts
+ * the output tiling slightly, but the distortion is consistent, and
+ * doesn't accumulate over a large patch of tiling, so it won't make
+ * anything end up totally out of place.
+ *
+ * (However, we compute the subsequent combination of 1 and sqrt(5)
+ * exactly, because using an approximation to sqrt(5) _could_ mean
+ * that in a sufficiently large patch of tiling two such combinations
+ * ended up misordered.)
+ *
+ * Adding to the confusion, we also flip round the x and y
+ * coordinates, because it's slightly nicer to have vertical edges in
+ * the tiling rather than horizontal ones. (Both for aesthetics, and
+ * also because if two P3 thin rhombs are separated by a horizontal
+ * line and both contain numeric clues then the clue numbers look a
+ * bit crowded, due to digits being taller than they are wide.)
+ *
+ * Finally, we have different base unit sizes for the two tiling
+ * types, because sensible sizes for the two are actually different.
+ * Each of P2 and P3 can be subdivided into the other, via dividing
+ * the larger triangle type in two, so that L large and S small become
+ * L+S large and L small. In the limit, this means that you expect the
+ * number of triangles (hence tiles) to grow by a factor of phi in
+ * each of those subdivisions (and hence by a factor of phi^2 in a
+ * full subdivision of P2 to a finer P2). So a sensible size ratio
+ * between the two tilings is one that makes them fit about the same
+ * number of tiles into the same area - and since tile area is
+ * proportional to _square_ of length, it follows that the P2 and P3
+ * length unit should differ by a factor of sqrt(phi).
+ */
+#define PENROSE_XUNIT_P2 37
+#define PENROSE_YUNIT_P2 44
+#define PENROSE_XUNIT_P3 30
+#define PENROSE_YUNIT_P3 35
+
+struct size { int w, h; };+static struct size api_size_penrose(int width, int height, int which)
+{+ int xunit = (which == PENROSE_P2 ? PENROSE_XUNIT_P2 : PENROSE_XUNIT_P3);
+ int yunit = (which == PENROSE_P2 ? PENROSE_YUNIT_P2 : PENROSE_YUNIT_P3);
+ struct size size = {+ width * PENROSE_TILESIZE / yunit,
+ height * PENROSE_TILESIZE / xunit,
+ };
+ return size;
+}
+
+static grid *grid_new_penrose(int width, int height, int which,
+ const char *desc); /* forward reference */
+
+static char *grid_new_desc_penrose(grid_type type, int width, int height,
+ random_state *rs)
+{+ char *buf;
+ struct PenrosePatchParams params;
+ int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
+ struct size size = api_size_penrose(width, height, which);
+
+ penrose_tiling_randomise(¶ms, which, size.h, size.w, rs);
+
+ buf = snewn(params.ncoords + 3, char);
+ buf[0] = '0' + params.orientation;
+ buf[1] = '0' + params.start_vertex;
+ memcpy(buf + 2, params.coords, params.ncoords);
+ buf[2 + params.ncoords] = '\0';
+
+ sfree(params.coords);
+ return buf;
+}
+
+/* Shared code between validating and reading grid descs.
+ * Always allocates params->coords, whether or not it returns an error. */
+static const char *grid_desc_to_penrose_params(
+ const char *desc, int which, struct PenrosePatchParams *params)
+{+ int i;
+
+ if (!*desc)
+ return "empty grid description";
+
+ params->ncoords = strlen(desc) - 2;
+ params->coords = snewn(params->ncoords, char);
+
+ {+ char c = desc[0];
+ if (isdigit((unsigned char)c))
+ params->orientation = c - '0';
+ else
+ return "expected digit at start of grid description";
+
+ c = desc[1];
+ if (c >= '0' && c < '3')
+ params->start_vertex = c - '0';
+ else
+ return "expected digit as second char of grid description";
+ }
+
+ for (i = 0; i < params->ncoords; i++) {+ char c = desc[i+2];
+ if (!penrose_valid_letter(c, which))
+ return "expected tile letter in grid description";
+ params->coords[i] = c;
+ }
+
+ return NULL;
+}
+
+static const char *grid_validate_desc_penrose(grid_type type,
+ int width, int height,
+ const char *desc)
+{+ struct PenrosePatchParams params;
+ const char *error = NULL;
+ int which = (type == GRID_PENROSE_P2 ? PENROSE_P2 : PENROSE_P3);
+
+ if (!desc)
+ return "Missing grid description string.";
+
+ if (*desc == 'G')
+ return grid_validate_desc_penrose_legacy(type, width, height, desc);
+
+ error = grid_desc_to_penrose_params(desc, which, ¶ms);
+ if (!error)
+ error = penrose_tiling_params_invalid(¶ms, which);
+
+ sfree(params.coords);
+ return error;
+}
+
+struct penrosecontext {+ grid *g;
+ tree234 *points;
+ int xunit, yunit;
+};
+
+static void grid_penrose_callback(void *vctx, const int *coords)
+{+ struct penrosecontext *ctx = (struct penrosecontext *)vctx;
+ size_t i;
+
+ grid_face_add_new(ctx->g, 4);
+ for (i = 0; i < 4; i++) {+ grid_dot *d = grid_get_dot(
+ ctx->g, ctx->points,
+ coords[4*i+2] * ctx->yunit + n_times_root_k(
+ coords[4*i+3] * ctx->yunit, 5),
+ coords[4*i+0] * ctx->xunit + n_times_root_k(
+ coords[4*i+1] * ctx->xunit, 5));
+ grid_face_set_dot(ctx->g, d, i);
+ }
+}
+
+static grid *grid_new_penrose(int width, int height, int which,
+ const char *desc)
+{+ struct PenrosePatchParams params;
+ const char *error = NULL;
+ struct penrosecontext ctx[1];
+ struct size size;
+
+ if (*desc == 'G')
+ return grid_new_penrose_legacy(width, height, which, desc);
+
+ error = grid_desc_to_penrose_params(desc, which, ¶ms);
+ assert(error == NULL && "grid_validate_desc_penrose should have failed");
+
+ ctx->g = grid_empty();
+ ctx->g->tilesize = PENROSE_TILESIZE;
+
+ ctx->points = newtree234(grid_point_cmp_fn);
+
+ ctx->xunit = (which == PENROSE_P2 ? PENROSE_XUNIT_P2 : PENROSE_XUNIT_P3);
+ ctx->yunit = (which == PENROSE_P2 ? PENROSE_YUNIT_P2 : PENROSE_YUNIT_P3);
+
+ size = api_size_penrose(width, height, which);
+ penrose_tiling_generate(¶ms, size.h, size.w,
+ grid_penrose_callback, ctx);
+
+ freetree234(ctx->points);
+ sfree(params.coords);
+
+ grid_trim_vigorously(ctx->g);
+ grid_make_consistent(ctx->g);
+
+ /*
+ * Centre the grid in its originally promised rectangle.
+ */
+ {+ int w = width * PENROSE_TILESIZE, h = height * PENROSE_TILESIZE;
+ ctx->g->lowest_x -= (w - (ctx->g->highest_x - ctx->g->lowest_x))/2;
+ ctx->g->lowest_y -= (h - (ctx->g->highest_y - ctx->g->lowest_y))/2;
+ ctx->g->highest_x = ctx->g->lowest_x + w;
+ ctx->g->highest_y = ctx->g->lowest_y + h;
+ }
+
+ return ctx->g;
}
static const char *grid_validate_params_penrose_p2_kite(int width, int height)
--- /dev/null
+++ b/penrose-internal.h
@@ -1,0 +1,289 @@
+#include "penrose.h"
+
+static inline unsigned num_subtriangles(char t)
+{+ return (t == 'A' || t == 'B' || t == 'X' || t == 'Y') ? 3 : 2;
+}
+
+static inline unsigned sibling_edge(char t)
+{+ switch (t) {+ case 'A': case 'U': return 2;
+ case 'B': case 'V': return 1;
+ default: return 0;
+ }
+}
+
+/*
+ * Coordinate system for tracking Penrose-tile half-triangles.
+ * PenroseCoords simply stores an array of triangle types.
+ */
+typedef struct PenroseCoords {+ char *c;
+ size_t nc, csize;
+} PenroseCoords;
+
+PenroseCoords *penrose_coords_new(void);
+void penrose_coords_free(PenroseCoords *pc);
+void penrose_coords_make_space(PenroseCoords *pc, size_t size);
+PenroseCoords *penrose_coords_copy(PenroseCoords *pc_in);
+
+/*
+ * Coordinate system for locating Penrose tiles in the plane.
+ *
+ * The 'Point' structure represents a single point by means of an
+ * integer linear combination of {1, t, t^2, t^3}, where t is the+ * complex number exp(i pi/5) representing 1/10 of a turn about the
+ * origin.
+ *
+ * The 'PenroseTriangle' structure represents a half-tile triangle,
+ * giving both the locations of its vertices and its combinatorial
+ * coordinates. It also contains a linked-list pointer and a boolean
+ * flag, used during breadth-first search to generate all the tiles in
+ * an area and report them exactly once.
+ */
+typedef struct Point {+ int coeffs[4];
+} Point;
+typedef struct PenroseTriangle PenroseTriangle;
+struct PenroseTriangle {+ Point vertices[3];
+ PenroseCoords *pc;
+ PenroseTriangle *next; /* used in breadth-first search */
+ bool reported;
+};
+
+/* Fill in all the coordinates of a triangle starting from any single edge.
+ * Requires tri->pc to have been filled in, so that we know which shape of
+ * triangle we're placing. */
+void penrose_place(PenroseTriangle *tri, Point u, Point v, int index_of_u);
+
+/* Free a PenroseHalf and its contained coordinates, or a whole PenroseTile */
+void penrose_free(PenroseTriangle *tri);
+
+/*
+ * A Point is really a complex number, so we can add, subtract and
+ * multiply them.
+ */
+static inline Point point_add(Point a, Point b)
+{+ Point r;
+ size_t i;
+ for (i = 0; i < 4; i++)
+ r.coeffs[i] = a.coeffs[i] + b.coeffs[i];
+ return r;
+}
+static inline Point point_sub(Point a, Point b)
+{+ Point r;
+ size_t i;
+ for (i = 0; i < 4; i++)
+ r.coeffs[i] = a.coeffs[i] - b.coeffs[i];
+ return r;
+}
+static inline Point point_mul_by_t(Point x)
+{+ Point r;
+ /* Multiply by t by using the identity t^4 - t^3 + t^2 - t + 1 = 0,
+ * so t^4 = t^3 - t^2 + t - 1 */
+ r.coeffs[0] = -x.coeffs[3];
+ r.coeffs[1] = x.coeffs[0] + x.coeffs[3];
+ r.coeffs[2] = x.coeffs[1] - x.coeffs[3];
+ r.coeffs[3] = x.coeffs[2] + x.coeffs[3];
+ return r;
+}
+static inline Point point_mul(Point a, Point b)
+{+ size_t i, j;
+ Point r;
+
+ /* Initialise r to be a, scaled by b's t^3 term */
+ for (j = 0; j < 4; j++)
+ r.coeffs[j] = a.coeffs[j] * b.coeffs[3];
+
+ /* Now iterate r = t*r + (next coefficient down), by Horner's rule */
+ for (i = 3; i-- > 0 ;) {+ r = point_mul_by_t(r);
+ for (j = 0; j < 4; j++)
+ r.coeffs[j] += a.coeffs[j] * b.coeffs[i];
+ }
+
+ return r;
+}
+static inline bool point_equal(Point a, Point b)
+{+ size_t i;
+ for (i = 0; i < 4; i++)
+ if (a.coeffs[i] != b.coeffs[i])
+ return false;
+ return true;
+}
+
+/*
+ * Return the Point corresponding to a rotation of s steps around the
+ * origin, i.e. a rotation by 36*s degrees or s*pi/5 radians.
+ */
+static inline Point point_rot(int s)
+{+ Point r = {{ 1, 0, 0, 0 }};+ Point tpower = {{ 0, 1, 0, 0 }};+
+ /* Reduce to a sensible range */
+ s = s % 10;
+ if (s < 0)
+ s += 10;
+
+ while (true) {+ if (s & 1)
+ r = point_mul(r, tpower);
+ s >>= 1;
+ if (!s)
+ break;
+ tpower = point_mul(tpower, tpower);
+ }
+
+ return r;
+}
+
+/*
+ * PenroseContext is the shared context of a whole run of the
+ * algorithm. Its 'prototype' PenroseCoords object represents the
+ * coordinates of the starting triangle, and is extended as necessary;
+ * any other PenroseCoord that needs extending will copy the
+ * higher-order values from ctx->prototype as needed, so that once
+ * each choice has been made, it remains consistent.
+ *
+ * When we're inventing a random piece of tiling in the first place,
+ * we append to ctx->prototype by choosing a random (but legal)
+ * higher-level metatile for the current topmost one to turn out to be
+ * part of. When we're replaying a generation whose parameters are
+ * already stored, we don't have a random_state, and we make fixed
+ * decisions if not enough coordinates were provided, as in the
+ * corresponding hat.c system.
+ *
+ * For a normal (non-testing) caller, penrosectx_generate() is the
+ * main useful function. It breadth-first searches a whole area to
+ * generate all the triangles in it, starting from a (typically
+ * central) one with the coordinates of ctx->prototype. It takes two
+ * callback function: one that checks whether a triangle is within the
+ * bounds of the target area (and therefore the search should continue
+ * exploring its neighbours), and another that reports a full Penrose
+ * tile once both of its halves have been found and determined to be
+ * in bounds.
+ */
+typedef struct PenroseContext {+ random_state *rs;
+ bool must_free_rs;
+ unsigned start_vertex; /* which vertex of 'prototype' is at the origin? */
+ int orientation; /* orientation to put in PenrosePatchParams */
+ PenroseCoords *prototype;
+} PenroseContext;
+
+void penrosectx_init_random(PenroseContext *ctx, random_state *rs, int which);
+void penrosectx_init_from_params(
+ PenroseContext *ctx, const struct PenrosePatchParams *ps);
+void penrosectx_cleanup(PenroseContext *ctx);
+PenroseCoords *penrosectx_initial_coords(PenroseContext *ctx);
+void penrosectx_extend_coords(PenroseContext *ctx, PenroseCoords *pc,
+ size_t n);
+void penrosectx_step(PenroseContext *ctx, PenroseCoords *pc,
+ unsigned edge, unsigned *outedge);
+void penrosectx_generate(
+ PenroseContext *ctx,
+ bool (*inbounds)(void *inboundsctx,
+ const PenroseTriangle *tri), void *inboundsctx,
+ void (*tile)(void *tilectx, const Point *vertices), void *tilectx);
+
+/* Subroutines that step around the tiling specified by a PenroseCtx,
+ * delivering both plane and combinatorial coordinates as they go */
+PenroseTriangle *penrose_initial(PenroseContext *ctx);
+PenroseTriangle *penrose_adjacent(PenroseContext *ctx,
+ const PenroseTriangle *src_spec,
+ unsigned src_edge, unsigned *dst_edge);
+
+/* For extracting the point coordinates */
+typedef struct Coord {+ int c1, cr5; /* coefficients of 1 and sqrt(5) respectively */
+} Coord;
+
+static inline Coord point_x(Point p)
+{+ Coord x = {+ 4 * p.coeffs[0] + p.coeffs[1] - p.coeffs[2] + p.coeffs[3],
+ p.coeffs[1] + p.coeffs[2] - p.coeffs[3],
+ };
+ return x;
+}
+
+static inline Coord point_y(Point p)
+{+ Coord y = {+ 2 * p.coeffs[1] + p.coeffs[2] + p.coeffs[3],
+ p.coeffs[2] + p.coeffs[3],
+ };
+ return y;
+}
+
+static inline int coord_sign(Coord x)
+{+ if (x.c1 == 0 && x.cr5 == 0)
+ return 0;
+ if (x.c1 >= 0 && x.cr5 >= 0)
+ return +1;
+ if (x.c1 <= 0 && x.cr5 <= 0)
+ return -1;
+
+ if (x.c1 * x.c1 > 5 * x.cr5 * x.cr5)
+ return x.c1 < 0 ? -1 : +1;
+ else
+ return x.cr5 < 0 ? -1 : +1;
+}
+
+static inline Coord coord_construct(int c1, int cr5)
+{+ Coord c = { c1, cr5 };+ return c;
+}
+
+static inline Coord coord_integer(int c1)
+{+ return coord_construct(c1, 0);
+}
+
+static inline Coord coord_add(Coord a, Coord b)
+{+ Coord sum;
+ sum.c1 = a.c1 + b.c1;
+ sum.cr5 = a.cr5 + b.cr5;
+ return sum;
+}
+
+static inline Coord coord_sub(Coord a, Coord b)
+{+ Coord diff;
+ diff.c1 = a.c1 - b.c1;
+ diff.cr5 = a.cr5 - b.cr5;
+ return diff;
+}
+
+static inline Coord coord_mul(Coord a, Coord b)
+{+ Coord prod;
+ prod.c1 = a.c1 * b.c1 + 5 * a.cr5 * b.cr5;
+ prod.cr5 = a.c1 * b.cr5 + a.cr5 * b.c1;
+ return prod;
+}
+
+static inline Coord coord_abs(Coord a)
+{+ int sign = coord_sign(a);
+ Coord abs;
+ abs.c1 = a.c1 * sign;
+ abs.cr5 = a.cr5 * sign;
+ return abs;
+}
+
+static inline int coord_cmp(Coord a, Coord b)
+{+ return coord_sign(coord_sub(a, b));
+}
--- /dev/null
+++ b/penrose-legacy.c
@@ -1,0 +1,506 @@
+/* penrose-legacy.c: legacy Penrose tile generator.
+ *
+ * Works by choosing a small patch from a recursively expanded large
+ * region of tiling, using one of the two algorithms described at
+ *
+ * https://www.chiark.greenend.org.uk/~sgtatham/quasiblog/aperiodic-tilings/
+ *
+ * This method of generating Penrose tiling fragments is superseded by
+ * the completely different algorithm in penrose.c, using the other
+ * algorithm in that article (the 'combinatorial coordinates' one). We
+ * keep the legacy algorithm around only for interpreting Loopy game
+ * IDs generated by older versions of the code.
+ */
+
+#include <assert.h>
+#include <string.h>
+#ifdef NO_TGMATH_H
+# include <math.h>
+#else
+# include <tgmath.h>
+#endif
+#include <stdio.h>
+
+#include "puzzles.h" /* for malloc routines, and PI */
+#include "penrose-legacy.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(void)
+{+ 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(void)
+{+ 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 penrose_legacy_vx(vector *vs, int i)
+{+ return (vs[i].a + vs[i].d) * COS54 +
+ (vs[i].b + vs[i].c) * COS18;
+}
+
+double penrose_legacy_vy(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.
+ */
+
+static 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_legacy_state *state, int depth, int flip,
+ vector v_orig, vector v_edge);
+
+static int penrose_p2_large(penrose_legacy_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_legacy_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_legacy_state *state, int depth, int flip,
+ vector v_orig, vector v_edge);
+
+static int penrose_p3_large(penrose_legacy_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_legacy_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_legacy_side_length(double start_size, int depth)
+{+ return start_size / pow(PHI, depth);
+}
+
+/*
+ * 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_legacy(penrose_legacy_state *state, int which, int angle)
+{+ 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;
+
+ vo = v_rotate(vo, angle);
+ vb = v_rotate(vb, angle);
+
+ 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_legacy_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;
+}
--- /dev/null
+++ b/penrose-legacy.h
@@ -1,0 +1,63 @@
+/* penrose-legacy.h: legacy 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.
+ *
+ * This method of generating Penrose tiling fragments is superseded by
+ * the completely different algorithm in penrose.c. We keep the legacy
+ * algorithm around only for interpreting Loopy game IDs generated by
+ * older versions of the code.
+ */
+
+#ifndef PUZZLES_PENROSE_LEGACY_H
+#define PUZZLES_PENROSE_LEGACY_H
+
+#ifndef PHI
+#define PHI 1.6180339887
+#endif
+
+typedef struct vector vector;
+
+double penrose_legacy_vx(vector *vs, int i);
+double penrose_legacy_vy(vector *vs, int i);
+
+typedef struct penrose_legacy_state penrose_legacy_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_legacy_state *, vector *, int, int);
+
+struct penrose_legacy_state {+ int start_size; /* initial side length */
+ int max_depth; /* Recursion depth */
+
+ tile_callback new_tile;
+ void *ctx; /* for callback */
+};
+
+#ifndef PENROSE_ENUM_DEFINED
+#define PENROSE_ENUM_DEFINED
+enum { PENROSE_P2, PENROSE_P3 };+#endif
+
+extern int penrose_legacy(penrose_legacy_state *state, int which, int angle);
+
+/* Returns the side-length of a penrose tile at recursion level
+ * gen, given a starting side length. */
+extern double penrose_legacy_side_length(double start_size, int depth);
+
+/* 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_legacy_calculate_size(
+ int which, int tilesize, int w, int h,
+ double *required_radius, int *start_size, int *depth);
+
+#endif
--- a/penrose.c
+++ b/penrose.c
@@ -1,501 +1,894 @@
-/* penrose.c
+/*
+ * Generate Penrose tilings via combinatorial coordinates.
*
- * Penrose tile generator.
- *
- * Works by choosing a small patch from a recursively expanded large
- * region of tiling, using one of the two algorithms described at
- *
+ * For general explanation of the algorithm:
* https://www.chiark.greenend.org.uk/~sgtatham/quasiblog/aperiodic-tilings/
+ *
+ * I use exactly the same indexing system here that's described in the
+ * article. For the P2 tiling, acute isosceles triangles (half-kites)
+ * are assigned letters A,B, and obtuse ones (half-darts) U,V; for P3,
+ * acute triangles (half of a thin rhomb) are C,D and obtuse ones
+ * (half a thick rhomb) are X,Y. Edges of all triangles are indexed
+ * anticlockwise around the triangle, with 0 being the base and 1,2
+ * being the two equal legs.
*/
#include <assert.h>
+#include <stddef.h>
#include <string.h>
-#ifdef NO_TGMATH_H
-# include <math.h>
-#else
-# include <tgmath.h>
-#endif
-#include <stdio.h>
-#include "puzzles.h" /* for malloc routines, and PI */
+#include "puzzles.h"
#include "penrose.h"
+#include "penrose-internal.h"
+#include "tree234.h"
-/* -------------------------------------------------------
- * 36-degree basis vector arithmetic routines.
- */
+bool penrose_valid_letter(char c, int which)
+{+ switch (c) {+ case 'A': case 'B': case 'U': case 'V':
+ return which == PENROSE_P2;
+ case 'C': case 'D': case 'X': case 'Y':
+ return which == PENROSE_P3;
+ default:
+ return false;
+ }
+}
-/* Imagine drawing a
- * ten-point 'clock face' like this:
+/*
+ * Result of attempting a transition within the coordinate system.
+ * INTERNAL means we've moved to a different child of the same parent,
+ * so the 'internal' substructure gives the type of the new triangle
+ * and which edge of it we came in through; EXTERNAL means we've moved
+ * out of the parent entirely, and the 'external' substructure tells
+ * us which edge of the parent triangle we left by, and if it's
+ * divided in two, which end of that edge (-1 for the left end or +1
+ * for the right end). If the parent edge is undivided, end == 0.
*
- * -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.)
+ * The type FAIL _shouldn't_ ever come up! It occurs if you try to
+ * compute an incoming transition with an illegal value of 'end' (i.e.
+ * having the wrong idea of whether the edge is divided), or if you
+ * refer to a child triangle type that doesn't exist in that parent.
+ * If it ever happens in the production code then an assertion will
+ * fail. But it might be useful to other users of the same code.
*/
+typedef struct TransitionResult {+ enum { INTERNAL, EXTERNAL, FAIL } type;+ union {+ struct {+ char new_child;
+ unsigned char new_edge;
+ } internal;
+ struct {+ unsigned char parent_edge;
+ signed char end;
+ } external;
+ } u;
+} TransitionResult;
-struct vector { int a, b, c, d; };+/* Construction function to make an INTERNAL-type TransitionResult */
+static inline TransitionResult internal(char new_child, unsigned new_edge)
+{+ TransitionResult tr;
+ tr.type = INTERNAL;
+ tr.u.internal.new_child = new_child;
+ tr.u.internal.new_edge = new_edge;
+ return tr;
+}
-static vector v_origin(void)
+/* Construction function to make an EXTERNAL-type TransitionResult */
+static inline TransitionResult external(unsigned parent_edge, int end)
{- vector v;
- v.a = v.b = v.c = v.d = 0;
- return v;
+ TransitionResult tr;
+ tr.type = EXTERNAL;
+ tr.u.external.parent_edge = parent_edge;
+ tr.u.external.end = end;
+ return tr;
}
-/* 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(void)
+/* Construction function to make a FAIL-type TransitionResult */
+static inline TransitionResult fail(void)
{- vector v;
+ TransitionResult tr;
+ tr.type = FAIL;
+ return tr;
+}
- v.b = 1;
- v.a = v.c = v.d = 0;
- return v;
+/*
+ * Compute a transition out of a triangle. Can return either INTERNAL
+ * or EXTERNAL types (or FAIL if it gets invalid data).
+ */
+static TransitionResult transition(char parent, char child, unsigned edge)
+{+ switch (parent) {+ case 'A':
+ switch (child) {+ case 'A':
+ switch (edge) {+ case 0: return external(2, -1);
+ case 1: return external(0, 0);
+ case 2: return internal('B', 1);+ }
+ case 'B':
+ switch (edge) {+ case 0: return internal('U', 1);+ case 1: return internal('A', 2);+ case 2: return external(1, +1);
+ }
+ case 'U':
+ switch (edge) {+ case 0: return external(2, +1);
+ case 1: return internal('B', 0);+ case 2: return external(1, -1);
+ }
+ default:
+ return fail();
+ }
+ case 'B':
+ switch (child) {+ case 'A':
+ switch (edge) {+ case 0: return internal('V', 2);+ case 1: return external(2, -1);
+ case 2: return internal('B', 1);+ }
+ case 'B':
+ switch (edge) {+ case 0: return external(1, +1);
+ case 1: return internal('A', 2);+ case 2: return external(0, 0);
+ }
+ case 'V':
+ switch (edge) {+ case 0: return external(1, -1);
+ case 1: return external(2, +1);
+ case 2: return internal('A', 0);+ }
+ default:
+ return fail();
+ }
+ case 'U':
+ switch (child) {+ case 'B':
+ switch (edge) {+ case 0: return internal('U', 1);+ case 1: return external(2, 0);
+ case 2: return external(0, +1);
+ }
+ case 'U':
+ switch (edge) {+ case 0: return external(1, 0);
+ case 1: return internal('B', 0);+ case 2: return external(0, -1);
+ }
+ default:
+ return fail();
+ }
+ case 'V':
+ switch (child) {+ case 'A':
+ switch (edge) {+ case 0: return internal('V', 2);+ case 1: return external(0, -1);
+ case 2: return external(1, 0);
+ }
+ case 'V':
+ switch (edge) {+ case 0: return external(2, 0);
+ case 1: return external(0, +1);
+ case 2: return internal('A', 0);+ }
+ default:
+ return fail();
+ }
+ case 'C':
+ switch (child) {+ case 'C':
+ switch (edge) {+ case 0: return external(1, +1);
+ case 1: return internal('Y', 1);+ case 2: return external(0, 0);
+ }
+ case 'Y':
+ switch (edge) {+ case 0: return external(2, 0);
+ case 1: return internal('C', 1);+ case 2: return external(1, -1);
+ }
+ default:
+ return fail();
+ }
+ case 'D':
+ switch (child) {+ case 'D':
+ switch (edge) {+ case 0: return external(2, -1);
+ case 1: return external(0, 0);
+ case 2: return internal('X', 2);+ }
+ case 'X':
+ switch (edge) {+ case 0: return external(1, 0);
+ case 1: return external(2, +1);
+ case 2: return internal('D', 2);+ }
+ default:
+ return fail();
+ }
+ case 'X':
+ switch (child) {+ case 'C':
+ switch (edge) {+ case 0: return external(2, +1);
+ case 1: return internal('Y', 1);+ case 2: return internal('X', 1);+ }
+ case 'X':
+ switch (edge) {+ case 0: return external(1, 0);
+ case 1: return internal('C', 2);+ case 2: return external(0, -1);
+ }
+ case 'Y':
+ switch (edge) {+ case 0: return external(0, +1);
+ case 1: return internal('C', 1);+ case 2: return external(2, -1);
+ }
+ default:
+ return fail();
+ }
+ case 'Y':
+ switch (child) {+ case 'D':
+ switch (edge) {+ case 0: return external(1, -1);
+ case 1: return internal('Y', 2);+ case 2: return internal('X', 2);+ }
+ case 'X':
+ switch (edge) {+ case 0: return external(0, -1);
+ case 1: return external(1, +1);
+ case 2: return internal('D', 2);+ }
+ case 'Y':
+ switch (edge) {+ case 0: return external(2, 0);
+ case 1: return external(0, +1);
+ case 2: return internal('D', 1);+ }
+ default:
+ return fail();
+ }
+ default:
+ return fail();
+ }
}
-#endif
+/*
+ * Compute a transition into a parent triangle, after the above
+ * function reported an EXTERNAL transition out of a neighbouring
+ * parent and we had to recurse. Because we're coming inwards, this
+ * should always return an INTERNAL TransitionResult (or FAIL if it
+ * gets invalid data).
+ */
+static TransitionResult transition_in(char parent, unsigned edge, int end)
+{+ #define EDGEEND(edge, end) (3 * (edge) + 1 + (end))
-#define COS54 0.5877852
-#define SIN54 0.8090169
-#define COS18 0.9510565
-#define SIN18 0.3090169
+ switch (parent) {+ case 'A':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, 0): return internal('A', 1);+ case EDGEEND(1, -1): return internal('B', 2);+ case EDGEEND(1, +1): return internal('U', 2);+ case EDGEEND(2, -1): return internal('U', 0);+ case EDGEEND(2, +1): return internal('A', 0);+ default:
+ return fail();
+ }
+ case 'B':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, 0): return internal('B', 2);+ case EDGEEND(1, -1): return internal('B', 0);+ case EDGEEND(1, +1): return internal('V', 0);+ case EDGEEND(2, -1): return internal('V', 1);+ case EDGEEND(2, +1): return internal('A', 1);+ default:
+ return fail();
+ }
+ case 'U':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, -1): return internal('B', 2);+ case EDGEEND(0, +1): return internal('U', 2);+ case EDGEEND(1, 0): return internal('U', 0);+ case EDGEEND(2, 0): return internal('B', 1);+ default:
+ return fail();
+ }
+ case 'V':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, -1): return internal('V', 1);+ case EDGEEND(0, +1): return internal('A', 1);+ case EDGEEND(1, 0): return internal('A', 2);+ case EDGEEND(2, 0): return internal('V', 0);+ default:
+ return fail();
+ }
+ case 'C':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, 0): return internal('C', 2);+ case EDGEEND(1, -1): return internal('C', 0);+ case EDGEEND(1, +1): return internal('Y', 2);+ case EDGEEND(2, 0): return internal('Y', 0);+ default:
+ return fail();
+ }
+ case 'D':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, 0): return internal('D', 1);+ case EDGEEND(1, 0): return internal('X', 0);+ case EDGEEND(2, -1): return internal('X', 1);+ case EDGEEND(2, +1): return internal('D', 0);+ default:
+ return fail();
+ }
+ case 'X':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, -1): return internal('Y', 0);+ case EDGEEND(0, +1): return internal('X', 2);+ case EDGEEND(1, 0): return internal('X', 0);+ case EDGEEND(2, -1): return internal('C', 0);+ case EDGEEND(2, +1): return internal('Y', 2);+ default:
+ return fail();
+ }
+ case 'Y':
+ switch (EDGEEND(edge, end)) {+ case EDGEEND(0, +1): return internal('X', 0);+ case EDGEEND(0, -1): return internal('Y', 1);+ case EDGEEND(1, -1): return internal('X', 1);+ case EDGEEND(1, +1): return internal('D', 0);+ case EDGEEND(2, 0): return internal('Y', 0);+ default:
+ return fail();
+ }
+ default:
+ return fail();
+ }
-/* 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;
+ #undef EDGEEND
}
-double v_y(vector *vs, int i)
+PenroseCoords *penrose_coords_new(void)
{- return (vs[i].a - vs[i].d) * SIN54 +
- (vs[i].b - vs[i].c) * SIN18;
+ PenroseCoords *pc = snew(PenroseCoords);
+ pc->nc = pc->csize = 0;
+ pc->c = NULL;
+ return pc;
+}
+void penrose_coords_free(PenroseCoords *pc)
+{+ if (pc) {+ sfree(pc->c);
+ sfree(pc);
+ }
}
-static vector v_trans(vector v, vector trans)
+void penrose_coords_make_space(PenroseCoords *pc, size_t size)
{- v.a += trans.a;
- v.b += trans.b;
- v.c += trans.c;
- v.d += trans.d;
- return v;
+ if (pc->csize < size) {+ pc->csize = pc->csize * 5 / 4 + 16;
+ if (pc->csize < size)
+ pc->csize = size;
+ pc->c = sresize(pc->c, pc->csize, char);
+ }
}
-static vector v_rotate_36(vector v)
+PenroseCoords *penrose_coords_copy(PenroseCoords *pc_in)
{- 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;
+ PenroseCoords *pc_out = penrose_coords_new();
+ penrose_coords_make_space(pc_out, pc_in->nc);
+ memcpy(pc_out->c, pc_in->c, pc_in->nc * sizeof(*pc_out->c));
+ pc_out->nc = pc_in->nc;
+ return pc_out;
}
-static vector v_rotate(vector v, int ang)
+/*
+ * The main recursive function for computing the next triangle's
+ * combinatorial coordinates.
+ */
+static void penrosectx_step_recurse(
+ PenroseContext *ctx, PenroseCoords *pc, size_t depth,
+ unsigned edge, unsigned *outedge)
{- int i;
+ TransitionResult tr;
- 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;
-}
+ penrosectx_extend_coords(ctx, pc, depth+2);
-#ifdef TEST_VECTORS
+ /* Look up the transition out of the starting triangle */
+ tr = transition(pc->c[depth+1], pc->c[depth], edge);
-static vector v_scale(vector v, int sc)
+ /* If we've left the parent triangle, recurse to find out what new
+ * triangle we've landed in at the next size up, and then call
+ * transition_in to find out which child of that parent we're
+ * going to */
+ if (tr.type == EXTERNAL) {+ unsigned parent_outedge;
+ penrosectx_step_recurse(
+ ctx, pc, depth+1, tr.u.external.parent_edge, &parent_outedge);
+ tr = transition_in(pc->c[depth+1], parent_outedge, tr.u.external.end);
+ }
+
+ /* Now we should definitely have ended up in a child of the
+ * (perhaps rewritten) parent triangle */
+ assert(tr.type == INTERNAL);
+ pc->c[depth] = tr.u.internal.new_child;
+ *outedge = tr.u.internal.new_edge;
+}
+
+void penrosectx_step(PenroseContext *ctx, PenroseCoords *pc,
+ unsigned edge, unsigned *outedge)
{- v.a *= sc;
- v.b *= sc;
- v.c *= sc;
- v.d *= sc;
- return v;
+ /* Allow outedge to be NULL harmlessly, just in case */
+ unsigned dummy_outedge;
+ if (!outedge)
+ outedge = &dummy_outedge;
+
+ penrosectx_step_recurse(ctx, pc, 0, edge, outedge);
}
-#endif
+static Point penrose_triangle_post_edge(char c, unsigned edge)
+{+ static const Point acute_post_edge[3] = {+ {{-1, 1, 0, 1}}, /* phi * t^3 */+ {{-1, 1, -1, 1}}, /* t^4 */+ {{-1, 1, 0, 0}}, /* 1/phi * t^3 */+ };
+ static const Point obtuse_post_edge[3] = {+ {{0, -1, 1, 0}}, /* 1/phi * t^4 */+ {{0, 0, 1, 0}}, /* t^2 */+ {{-1, 0, 0, 1}}, /* phi * t^4 */+ };
-static vector v_growphi(vector v)
+ switch (c) {+ case 'A': case 'B': case 'C': case 'D':
+ return acute_post_edge[edge];
+ default: /* case 'U': case 'V': case 'X': case 'Y': */
+ return obtuse_post_edge[edge];
+ }
+}
+
+void penrose_place(PenroseTriangle *tri, Point u, Point v, int index_of_u)
{- 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;
+ unsigned i;
+ Point here = u, delta = point_sub(v, u);
+
+ for (i = 0; i < 3; i++) {+ unsigned edge = (index_of_u + i) % 3;
+ tri->vertices[edge] = here;
+ here = point_add(here, delta);
+ delta = point_mul(delta, penrose_triangle_post_edge(
+ tri->pc->c[0], edge));
+ }
}
-static vector v_shrinkphi(vector v)
+void penrose_free(PenroseTriangle *tri)
{- 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;
+ penrose_coords_free(tri->pc);
+ sfree(tri);
}
-#ifdef TEST_VECTORS
-
-static const char *v_debug(vector v)
+static bool penrose_relative_probability(char c)
{- 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;
+ /* Penrose tile probability ratios are always phi, so we can
+ * approximate that very well using two consecutive Fibonacci
+ * numbers */
+ switch (c) {+ case 'A': case 'B': case 'X': case 'Y':
+ return 165580141;
+ case 'C': case 'D': case 'U': case 'V':
+ return 102334155;
+ default:
+ return 0;
+ }
}
-#endif
+static char penrose_choose_random(const char *possibilities, random_state *rs)
+{+ const char *p;
+ unsigned long value, limit = 0;
-/* -------------------------------------------------------
- * Tiling routines.
- */
+ for (p = possibilities; *p; p++)
+ limit += penrose_relative_probability(*p);
+ value = random_upto(rs, limit);
+ for (p = possibilities; *p; p++) {+ unsigned long curr = penrose_relative_probability(*p);
+ if (value < curr)
+ return *p;
+ value -= curr;
+ }
+ assert(false && "Probability overflow!");
+ return possibilities[0];
+}
-static vector xform_coord(vector vo, int shrink, vector vtrans, int ang)
+static const char *penrose_starting_tiles(int which)
{- if (shrink < 0)
- vo = v_shrinkphi(vo);
- else if (shrink > 0)
- vo = v_growphi(vo);
+ return which == PENROSE_P2 ? "ABUV" : "CDXY";
+}
- vo = v_rotate(vo, ang);
- vo = v_trans(vo, vtrans);
+static const char *penrose_valid_parents(char tile)
+{+ switch (tile) {+ case 'A': return "ABV";
+ case 'B': return "ABU";
+ case 'U': return "AU";
+ case 'V': return "BV";
+ case 'C': return "CX";
+ case 'D': return "DY";
+ case 'X': return "DXY";
+ case 'Y': return "CXY";
+ default: return NULL;
+ }
+}
- return vo;
+void penrosectx_init_random(PenroseContext *ctx, random_state *rs, int which)
+{+ ctx->rs = rs;
+ ctx->must_free_rs = false;
+ ctx->prototype = penrose_coords_new();
+ penrose_coords_make_space(ctx->prototype, 1);
+ ctx->prototype->c[0] = penrose_choose_random(
+ penrose_starting_tiles(which), rs);
+ ctx->prototype->nc = 1;
+ ctx->start_vertex = random_upto(rs, 3);
+ ctx->orientation = random_upto(rs, 10);
}
+void penrosectx_init_from_params(
+ PenroseContext *ctx, const struct PenrosePatchParams *ps)
+{+ size_t i;
-#define XFORM(n,o,s,a) vs[(n)] = xform_coord(v_edge, (s), vs[(o)], (a))
+ ctx->rs = NULL;
+ ctx->must_free_rs = false;
+ ctx->prototype = penrose_coords_new();
+ penrose_coords_make_space(ctx->prototype, ps->ncoords);
+ for (i = 0; i < ps->ncoords; i++)
+ ctx->prototype->c[i] = ps->coords[i];
+ ctx->prototype->nc = ps->ncoords;
+ ctx->start_vertex = ps->start_vertex;
+ ctx->orientation = ps->orientation;
+}
-static int penrose_p2_small(penrose_state *state, int depth, int flip,
- vector v_orig, vector v_edge);
+void penrosectx_cleanup(PenroseContext *ctx)
+{+ if (ctx->must_free_rs)
+ random_free(ctx->rs);
+ penrose_coords_free(ctx->prototype);
+}
-static int penrose_p2_large(penrose_state *state, int depth, int flip,
- vector v_orig, vector v_edge)
+PenroseCoords *penrosectx_initial_coords(PenroseContext *ctx)
{- vector vv_orig, vv_edge;
+ return penrose_coords_copy(ctx->prototype);
+}
-#ifdef DEBUG_PENROSE
- {- vector vs[3];
- vs[0] = v_orig;
- XFORM(1, 0, 0, 0);
- XFORM(2, 0, 0, -36*flip);
+void penrosectx_extend_coords(PenroseContext *ctx, PenroseCoords *pc,
+ size_t n)
+{+ if (ctx->prototype->nc < n) {+ penrose_coords_make_space(ctx->prototype, n);
+ while (ctx->prototype->nc < n) {+ if (!ctx->rs) {+ /*
+ * For safety, similarly to spectre.c, we respond to a
+ * lack of available random_state by making a
+ * deterministic one.
+ */
+ ctx->rs = random_new("dummy", 5);+ ctx->must_free_rs = true;
+ }
- state->new_tile(state, vs, 3, depth);
+ ctx->prototype->c[ctx->prototype->nc] = penrose_choose_random(
+ penrose_valid_parents(ctx->prototype->c[ctx->prototype->nc-1]),
+ ctx->rs);
+ ctx->prototype->nc++;
+ }
}
-#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);
+ penrose_coords_make_space(pc, n);
+ while (pc->nc < n) {+ pc->c[pc->nc] = ctx->prototype->c[pc->nc];
+ pc->nc++;
}
- 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)
+static Point penrose_triangle_edge_0_length(char c)
{- vector vv_orig;
+ static const Point one = {{ 1, 0, 0, 0 }};+ static const Point phi = {{ 1, 0, 1, -1 }};+ static const Point invphi = {{ 0, 0, 1, -1 }};-#ifdef DEBUG_PENROSE
- {- vector vs[3];
- vs[0] = v_orig;
- XFORM(1, 0, 0, 0);
- XFORM(2, 0, -1, -36*flip);
+ switch (c) {+ /* P2 tiling: unit-length edges are the long edges, i.e. edges
+ * 1,2 of AB and edge 0 of UV. So AB have edge 0 short. */
+ case 'A': case 'B':
+ return invphi;
+ case 'U': case 'V':
+ return one;
- state->new_tile(state, vs, 3, depth);
+ /* P3 tiling: unit-length edges are edges 1,2 of everything,
+ * so CD have edge 0 short and XY have it long. */
+ case 'C': case 'D':
+ return invphi;
+ default: /* case 'X': case 'Y': */
+ return phi;
}
-#endif
+}
- if (flip > 0) {- vector vs[4];
+PenroseTriangle *penrose_initial(PenroseContext *ctx)
+{+ char type = ctx->prototype->c[0];
+ Point origin = {{ 0, 0, 0, 0 }};+ Point edge0 = penrose_triangle_edge_0_length(type);
+ Point negoffset;
+ size_t i;
+ PenroseTriangle *tri = snew(PenroseTriangle);
- vs[0] = v_orig;
- XFORM(1, 0, 0, -72);
- XFORM(2, 0, -1, -36);
- XFORM(3, 0, 0, 0);
+ /* Orient the triangle by deciding what vector edge #0 should traverse */
+ edge0 = point_mul(edge0, point_rot(ctx->orientation));
- state->new_tile(state, vs, 4, depth);
- }
+ /* Place the triangle at an arbitrary position, in that orientation */
+ tri->pc = penrose_coords_copy(ctx->prototype);
+ penrose_place(tri, origin, edge0, 0);
- if (depth >= state->max_depth) return 0;
+ /* Now translate so that the appropriate vertex is at the origin */
+ negoffset = tri->vertices[ctx->start_vertex];
+ for (i = 0; i < 3; i++)
+ tri->vertices[i] = point_sub(tri->vertices[i], negoffset);
- vv_orig = v_trans(v_orig, v_edge);
+ return tri;
+}
- penrose_p2_large(state, depth+1, -flip,
- v_orig, v_shrinkphi(v_rotate(v_edge, -36*flip)));
+PenroseTriangle *penrose_adjacent(PenroseContext *ctx,
+ const PenroseTriangle *src_tri,
+ unsigned src_edge, unsigned *dst_edge_out)
+{+ unsigned dst_edge;
+ PenroseTriangle *dst_tri = snew(PenroseTriangle);
+ dst_tri->pc = penrose_coords_copy(src_tri->pc);
+ penrosectx_step(ctx, dst_tri->pc, src_edge, &dst_edge);
+ penrose_place(dst_tri, src_tri->vertices[(src_edge+1) % 3],
+ src_tri->vertices[src_edge], dst_edge);
+ if (dst_edge_out)
+ *dst_edge_out = dst_edge;
+ return dst_tri;
+}
- penrose_p2_small(state, depth+1, flip,
- vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
+static int penrose_cmp(void *av, void *bv)
+{+ PenroseTriangle *a = (PenroseTriangle *)av, *b = (PenroseTriangle *)bv;
+ size_t i, j;
+ /* We should only ever need to compare the first two vertices of
+ * any triangle, because those force the rest */
+ for (i = 0; i < 2; i++) {+ for (j = 0; j < 4; j++) {+ int ac = a->vertices[i].coeffs[j], bc = b->vertices[i].coeffs[j];
+ if (ac < bc)
+ return -1;
+ if (ac > bc)
+ return +1;
+ }
+ }
+
return 0;
}
-static int penrose_p3_small(penrose_state *state, int depth, int flip,
- vector v_orig, vector v_edge);
+static unsigned penrose_sibling_edge_index(char c)
+{+ switch (c) {+ case 'A': case 'U': return 2;
+ case 'B': case 'V': return 1;
+ default: return 0;
+ }
+}
-static int penrose_p3_large(penrose_state *state, int depth, int flip,
- vector v_orig, vector v_edge)
+void penrosectx_generate(
+ PenroseContext *ctx,
+ bool (*inbounds)(void *inboundsctx,
+ const PenroseTriangle *tri), void *inboundsctx,
+ void (*tile)(void *tilectx, const Point *vertices), void *tilectx)
{- vector vv_orig;
+ tree234 *placed = newtree234(penrose_cmp);
+ PenroseTriangle *qhead = NULL, *qtail = NULL;
-#ifdef DEBUG_PENROSE
{- vector vs[3];
- vs[0] = v_orig;
- XFORM(1, 0, 1, 0);
- XFORM(2, 0, 0, -36*flip);
+ PenroseTriangle *tri = penrose_initial(ctx);
- state->new_tile(state, vs, 3, depth);
+ add234(placed, tri);
+
+ tri->next = NULL;
+ tri->reported = false;
+
+ if (inbounds(inboundsctx, tri))
+ qhead = qtail = tri;
}
-#endif
- if (flip > 0) {- vector vs[4];
+ while (qhead) {+ PenroseTriangle *tri = qhead;
+ unsigned edge;
+ unsigned sibling_edge = penrose_sibling_edge_index(tri->pc->c[0]);
- vs[0] = v_orig;
- XFORM(1, 0, 0, -36);
- XFORM(2, 0, 1, 0);
- XFORM(3, 0, 0, 36);
+ for (edge = 0; edge < 3; edge++) {+ PenroseTriangle *new_tri, *found_tri;
- state->new_tile(state, vs, 4, depth);
- }
- if (depth >= state->max_depth) return 0;
+ new_tri = penrose_adjacent(ctx, tri, edge, NULL);
- vv_orig = v_trans(v_orig, v_edge);
+ if (!inbounds(inboundsctx, new_tri)) {+ penrose_free(new_tri);
+ continue;
+ }
- penrose_p3_large(state, depth+1, -flip,
- vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
+ found_tri = find234(placed, new_tri, NULL);
+ if (found_tri) {+ if (edge == sibling_edge && !tri->reported &&
+ !found_tri->reported) {+ /*
+ * found_tri and tri are opposite halves of the
+ * same tile; both are in the tree, and haven't
+ * yet been reported as a completed tile.
+ */
+ unsigned new_sibling_edge = penrose_sibling_edge_index(
+ found_tri->pc->c[0]);
+ Point tilevertices[4] = {+ tri->vertices[(sibling_edge + 1) % 3],
+ tri->vertices[(sibling_edge + 2) % 3],
+ found_tri->vertices[(new_sibling_edge + 1) % 3],
+ found_tri->vertices[(new_sibling_edge + 2) % 3],
+ };
+ tile(tilectx, tilevertices);
- penrose_p3_small(state, depth+1, flip,
- vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
+ tri->reported = true;
+ found_tri->reported = true;
+ }
- vv_orig = v_trans(v_orig, v_growphi(v_edge));
+ penrose_free(new_tri);
+ continue;
+ }
- penrose_p3_large(state, depth+1, flip,
- vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
+ add234(placed, new_tri);
+ qtail->next = new_tri;
+ qtail = new_tri;
+ new_tri->next = NULL;
+ new_tri->reported = false;
+ }
+ qhead = qhead->next;
+ }
- return 0;
+ {+ PenroseTriangle *tri;
+ while ((tri = delpos234(placed, 0)) != NULL)
+ penrose_free(tri);
+ freetree234(placed);
+ }
}
-static int penrose_p3_small(penrose_state *state, int depth, int flip,
- vector v_orig, vector v_edge)
+const char *penrose_tiling_params_invalid(
+ const struct PenrosePatchParams *params, int which)
{- vector vv_orig;
+ size_t i;
-#ifdef DEBUG_PENROSE
- {- vector vs[3];
- vs[0] = v_orig;
- XFORM(1, 0, 0, 0);
- XFORM(2, 0, 0, -36*flip);
+ if (params->ncoords == 0)
+ return "expected at least one coordinate";
- state->new_tile(state, vs, 3, depth);
+ for (i = 0; i < params->ncoords; i++) {+ if (!penrose_valid_letter(params->coords[i], which))
+ return "invalid coordinate letter";
+ if (i > 0 && !strchr(penrose_valid_parents(params->coords[i-1]),
+ params->coords[i]))
+ return "invalid pair of consecutive coordinates";
}
-#endif
- if (flip > 0) {- vector vs[4];
+ return NULL;
+}
- vs[0] = v_orig;
- XFORM(1, 0, 0, -36);
- XFORM(3, 0, 0, 0);
- XFORM(2, 3, 0, -36);
+struct PenroseOutputCtx {+ int xoff, yoff;
+ Coord xmin, xmax, ymin, ymax;
- state->new_tile(state, vs, 4, depth);
- }
- if (depth >= state->max_depth) return 0;
+ penrose_tile_callback_fn external_cb;
+ void *external_cbctx;
+};
- /* NB these two are identical to the first two of p3_large. */
- vv_orig = v_trans(v_orig, v_edge);
+static bool inbounds(void *vctx, const PenroseTriangle *tri)
+{+ struct PenroseOutputCtx *octx = (struct PenroseOutputCtx *)vctx;
+ size_t i;
- penrose_p3_large(state, depth+1, -flip,
- vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
+ for (i = 0; i < 3; i++) {+ Coord x = point_x(tri->vertices[i]);
+ Coord y = point_y(tri->vertices[i]);
- penrose_p3_small(state, depth+1, flip,
- vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
+ if (coord_cmp(x, octx->xmin) < 0 || coord_cmp(x, octx->xmax) > 0 ||
+ coord_cmp(y, octx->ymin) < 0 || coord_cmp(y, octx->ymax) > 0)
+ return false;
+ }
- return 0;
+ return true;
}
-/* -------------------------------------------------------
- * Utility routines.
- */
-
-double penrose_side_length(double start_size, int depth)
+static void null_output_tile(void *vctx, const Point *vertices)
{- return start_size / pow(PHI, depth);
}
-/*
- * 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, int angle)
+static void really_output_tile(void *vctx, const Point *vertices)
{- vector vo = v_origin();
- vector vb = v_origin();
+ struct PenroseOutputCtx *octx = (struct PenroseOutputCtx *)vctx;
+ size_t i;
+ int coords[16];
- vo.b = vo.c = -state->start_size;
- vo = v_shrinkphi(v_shrinkphi(vo));
+ for (i = 0; i < 4; i++) {+ Coord x = point_x(vertices[i]);
+ Coord y = point_y(vertices[i]);
- vb.b = state->start_size;
+ coords[4*i + 0] = x.c1 + octx->xoff;
+ coords[4*i + 1] = x.cr5;
+ coords[4*i + 2] = y.c1 + octx->yoff;
+ coords[4*i + 3] = y.cr5;
+ }
- vo = v_rotate(vo, angle);
- vb = v_rotate(vb, angle);
+ octx->external_cb(octx->external_cbctx, coords);
+}
- if (which == PENROSE_P2)
- return penrose_p2_large(state, 0, 1, vo, vb);
- else
- return penrose_p3_small(state, 0, 1, vo, vb);
+static void penrose_set_bounds(struct PenroseOutputCtx *octx, int w, int h)
+{+ octx->xoff = w/2;
+ octx->yoff = h/2;
+ octx->xmin.c1 = -octx->xoff;
+ octx->xmax.c1 = -octx->xoff + w;
+ octx->ymin.c1 = octx->yoff - h;
+ octx->ymax.c1 = octx->yoff;
+ octx->xmin.cr5 = 0;
+ octx->xmax.cr5 = 0;
+ octx->ymin.cr5 = 0;
+ octx->ymax.cr5 = 0;
}
-/*
- * 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)
+void penrose_tiling_randomise(struct PenrosePatchParams *params, int which,
+ int w, int h, random_state *rs)
{- double rradius, size;
- int n = 0;
+ PenroseContext ctx[1];
+ struct PenroseOutputCtx octx[1];
- /*
- * 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;
+ penrose_set_bounds(octx, w, h);
- rradius = tilesize * 3.11 * sqrt((double)(w*w + h*h));
- size = tilesize;
+ penrosectx_init_random(ctx, rs, which);
+ penrosectx_generate(ctx, inbounds, octx, null_output_tile, NULL);
- while ((size * INCIRCLE_RADIUS) < rradius) {- n++;
- size = size * PHI;
- }
+ params->orientation = ctx->orientation;
+ params->start_vertex = ctx->start_vertex;
+ params->ncoords = ctx->prototype->nc;
+ params->coords = snewn(params->ncoords, char);
+ memcpy(params->coords, ctx->prototype->c, params->ncoords);
- *start_size = (int)size;
- *depth = n;
- *required_radius = rradius;
+ penrosectx_cleanup(ctx);
+}
+
+void penrose_tiling_generate(
+ const struct PenrosePatchParams *params, int w, int h,
+ penrose_tile_callback_fn cb, void *cbctx)
+{+ PenroseContext ctx[1];
+ struct PenroseOutputCtx octx[1];
+
+ penrose_set_bounds(octx, w, h);
+ octx->external_cb = cb;
+ octx->external_cbctx = cbctx;
+
+ penrosectx_init_from_params(ctx, params);
+ penrosectx_generate(ctx, inbounds, octx, really_output_tile, octx);
+ penrosectx_cleanup(ctx);
}
--- a/penrose.h
+++ b/penrose.h
@@ -1,56 +1,86 @@
-/* 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 PUZZLES_PENROSE_H
#define PUZZLES_PENROSE_H
-#ifndef PHI
-#define PHI 1.6180339887
+struct PenrosePatchParams {+ /*
+ * A patch of Penrose tiling is identified by giving
+ *
+ * - the coordinates of the starting triangle, using a
+ * combinatorial coordinate system
+ *
+ * - which vertex of that triangle is at the centre point of the
+ * tiling
+ *
+ * - the orientation of the triangle's base edge, as a number
+ * from 0 to 9, measured in tenths of a turn
+ *
+ * Coordinates are a sequence of letters. For a P2 tiling all
+ * letters are from the set {A,B,U,V}; for P3, {C,D,X,Y}.+ */
+ unsigned start_vertex;
+ int orientation;
+ size_t ncoords;
+ char *coords;
+};
+
+#ifndef PENROSE_ENUM_DEFINED
+#define PENROSE_ENUM_DEFINED
+enum { PENROSE_P2, PENROSE_P3 };#endif
-typedef struct vector vector;
+bool penrose_valid_letter(char c, int which);
-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).
+/*
+ * Fill in PenrosePatchParams with a randomly selected set of
+ * coordinates, in enough detail to generate a patch of tiling filling
+ * a w x h area.
*
- * Parameters are state, vector array, npoints, depth.
- * ctx is inside state.
+ * Units of measurement: the tiling will be oriented such that
+ * horizontal tile edges are possible (and therefore vertical ones are
+ * not). Therefore, all x-coordinates will be rational linear
+ * combinations of 1 and sqrt(5), and all y-coordinates will be
+ * sin(pi/5) times such a rational linear combination. By scaling up
+ * appropriately we can turn those rational combinations into
+ * _integer_ combinations, so we do. Therefore, w is measured in units
+ * of 1/4, and h is measured in units of sin(pi/5)/2, on a scale where
+ * a length of 1 corresponds to the legs of the acute isosceles
+ * triangles in the tiling (hence, the long edges of P2 kites and
+ * darts, or all edges of P3 rhombs).
+ *
+ * (In case it's useful, the y scale factor sin(pi/5)/2 is an
+ * algebraic number. Its minimal polynomial is 256x^4 - 80x^2 + 5.)
+ *
+ * The 'coords' field of the structure will be filled in with a new
+ * dynamically allocated array. Any previous pointer in that field
+ * will be overwritten.
*/
-typedef int (*tile_callback)(penrose_state *, vector *, int, int);
+void penrose_tiling_randomise(struct PenrosePatchParams *params, int which,
+ int w, int h, random_state *rs);
-struct penrose_state {- int start_size; /* initial side length */
- int max_depth; /* Recursion depth */
+/*
+ * Validate a PenrosePatchParams to ensure it contains no illegal
+ * coordinates. Returns NULL if it's acceptable, or an error string if
+ * not.
+ */
+const char *penrose_tiling_params_invalid(
+ const struct PenrosePatchParams *params, int which);
- tile_callback new_tile;
- void *ctx; /* for callback */
-};
+/*
+ * Generate the actual set of Penrose tiles from a PenrosePatchParams,
+ * passing each one to a callback. The callback receives the vertices
+ * of each point, in the form of an array of 4*4 integers. Each vertex
+ * is represented by four consecutive integers in this array, with the
+ * first two giving the x coordinate and the last two the y
+ * coordinate. Each pair of integers a,b represent a single coordinate
+ * whose value is a + b*sqrt(5). The units of measurement for x and y
+ * are as described above.
+ */
+typedef void (*penrose_tile_callback_fn)(void *ctx, const int *coords);
-enum { PENROSE_P2, PENROSE_P3 };+#define PENROSE_NVERTICES 4
-extern int penrose(penrose_state *state, int which, int angle);
-
-/* 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);
-
-/* 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);
+void penrose_tiling_generate(
+ const struct PenrosePatchParams *params, int w, int h,
+ penrose_tile_callback_fn cb, void *cbctx);
#endif