diff --git a/CMakeLists.txt b/CMakeLists.txt index d3015a2..4153efc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,8 @@ include(cmake/setup.cmake) 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 $) add_library(common $ hat.c spectre.c) diff --git a/auxiliary/CMakeLists.txt b/auxiliary/CMakeLists.txt index 88428a7..c97db18 100644 --- a/auxiliary/CMakeLists.txt +++ b/auxiliary/CMakeLists.txt @@ -5,6 +5,7 @@ cliprogram(hat-test hat-test.c) 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) diff --git a/auxiliary/penrose-legacy-test.c b/auxiliary/penrose-legacy-test.c new file mode 100644 index 0000000..db018db --- /dev/null +++ b/auxiliary/penrose-legacy-test.c @@ -0,0 +1,98 @@ +#include +#include + +#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("\n", col, col); + printf("", + 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("\ +\n\ +\n\ +\n\ +\n\n"); + + printf("\n"); + penrose_legacy(&ps, which, 0); + printf("\n"); + + printf("\n", + ntiles, nfinal); + + printf(""); + + return 0; +} diff --git a/auxiliary/penrose-test.c b/auxiliary/penrose-test.c index 468997d..b2092e5 100644 --- a/auxiliary/penrose-test.c +++ b/auxiliary/penrose-test.c @@ -1,95 +1,46 @@ #include -#include +#include #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("\n", col, col); - printf("", - 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); -} - -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("\ -\n\ -\n\ -\n\ -\n\n"); - - printf("\n"); - penrose(&ps, which, 0); - printf("\n"); - - printf("\n", - ntiles, nfinal); - - printf(""); - - return 0; + random_state *rs; + struct PenrosePatchParams params[1]; + struct testctx tctx[1]; + int w = 50, h = 40; + + 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); + + 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); } diff --git a/grid.c b/grid.c index 2b21103..90a9278 100644 --- 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,53 +3025,13 @@ static grid *grid_new_compassdodecagonal(int width, int height, const char *desc return g; } -typedef struct setface_ctx -{ - int xmin, xmax, ymin, ymax; - - grid *g; - tree234 *points; -} setface_ctx; - -static double round_int_nearest_away(double r) -{ - 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) -{ - setface_ctx *sf_ctx = (setface_ctx *)state->ctx; - int i; - int xs[4], ys[4]; - - 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_int_nearest_away(tx); - ys[i] = (int)round_int_nearest_away(ty); - - 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; -} +/* + * 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 @@ -3086,7 +3047,7 @@ static const char *grid_validate_params_penrose(int width, int height) } static void grid_size_penrose(int width, int height, - int *tilesize, int *xextent, int *yextent) + int *tilesize, int *xextent, int *yextent) { int l = PENROSE_TILESIZE; @@ -3095,64 +3056,66 @@ static void grid_size_penrose(int width, int height, *yextent = l * height; } -static grid *grid_new_penrose(int width, int height, int which, const char *desc); /* forward reference */ +/* + * 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; -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; + tree234 *points; +} penrose_legacy_set_faces_ctx; - 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 double round_int_nearest_away(double r) +{ + return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5); } -static const char *grid_validate_desc_penrose(grid_type type, - int width, int height, - const char *desc) +static int penrose_legacy_set_faces(penrose_legacy_state *state, vector *vs, + int n, int depth) +{ + penrose_legacy_set_faces_ctx *sf_ctx = + (penrose_legacy_set_faces_ctx *)state->ctx; + int i; + int xs[4], ys[4]; + + 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 = 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); + + 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_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, penrose_legacy_vx(vs, i), + penrose_legacy_vy(vs, i))); + } + + return 0; +} + +static grid *grid_new_penrose_legacy(int width, int height, int which, + const char *desc); + +static const char *grid_validate_desc_penrose_legacy( + grid_type type, int width, int height, const char *desc) { int tilesize = PENROSE_TILESIZE, startsz, depth, xoff, yoff, aoff, inner_radius; double outer_radius; @@ -3162,8 +3125,8 @@ static const char *grid_validate_desc_penrose(grid_type type, 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 @@ static const char *grid_validate_desc_penrose(grid_type type, * 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 @@ static const char *grid_validate_desc_penrose(grid_type type, 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 @@ static grid *grid_new_penrose(int width, int height, int which, const char *desc 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 @@ static grid *grid_new_penrose(int width, int height, int which, const char *desc 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); @@ -3280,6 +3238,213 @@ static grid *grid_new_penrose(int width, int height, int which, const char *desc 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) { return grid_validate_params_penrose(width, height); diff --git a/penrose-internal.h b/penrose-internal.h new file mode 100644 index 0000000..984cd35 --- /dev/null +++ b/penrose-internal.h @@ -0,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)); +} diff --git a/penrose-legacy.c b/penrose-legacy.c new file mode 100644 index 0000000..709d68d --- /dev/null +++ b/penrose-legacy.c @@ -0,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 +#include +#ifdef NO_TGMATH_H +# include +#else +# include +#endif +#include + +#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; +} diff --git a/penrose-legacy.h b/penrose-legacy.h new file mode 100644 index 0000000..514fd8b --- /dev/null +++ b/penrose-legacy.h @@ -0,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 diff --git a/penrose.c b/penrose.c index 1c45ddb..4d7dcc4 100644 --- a/penrose.c +++ b/penrose.c @@ -1,501 +1,894 @@ -/* penrose.c - * - * 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 +/* + * Generate Penrose tilings via combinatorial coordinates. * + * 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 +#include #include -#ifdef NO_TGMATH_H -# include -#else -# include -#endif -#include -#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. - */ - -/* 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) +bool penrose_valid_letter(char c, int which) { - 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 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. - */ - -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_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); + 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; } -#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); } /* - * 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. + * 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. * - * (later mail: this is an overestimate by about 5%) + * 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; -int penrose(penrose_state *state, int which, int angle) +/* Construction function to make an INTERNAL-type TransitionResult */ +static inline TransitionResult internal(char new_child, unsigned new_edge) { - vector vo = v_origin(); - vector vb = v_origin(); + TransitionResult tr; + tr.type = INTERNAL; + tr.u.internal.new_child = new_child; + tr.u.internal.new_edge = new_edge; + return tr; +} - vo.b = vo.c = -state->start_size; - vo = v_shrinkphi(v_shrinkphi(vo)); +/* Construction function to make an EXTERNAL-type TransitionResult */ +static inline TransitionResult external(unsigned parent_edge, int end) +{ + TransitionResult tr; + tr.type = EXTERNAL; + tr.u.external.parent_edge = parent_edge; + tr.u.external.end = end; + return tr; +} - 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); +/* Construction function to make a FAIL-type TransitionResult */ +static inline TransitionResult fail(void) +{ + TransitionResult tr; + tr.type = FAIL; + return tr; } /* - * 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. + * Compute a transition out of a triangle. Can return either INTERNAL + * or EXTERNAL types (or FAIL if it gets invalid data). */ -#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) +static TransitionResult transition(char parent, char child, unsigned edge) { - double rradius, size; - int n = 0; + 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(); + } +} - /* - * 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; +/* + * 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)) - rradius = tilesize * 3.11 * sqrt((double)(w*w + h*h)); - size = tilesize; - - while ((size * INCIRCLE_RADIUS) < rradius) { - n++; - size = size * PHI; + 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(); } - *start_size = (int)size; - *depth = n; - *required_radius = rradius; + #undef EDGEEND +} + +PenroseCoords *penrose_coords_new(void) +{ + 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); + } +} + +void penrose_coords_make_space(PenroseCoords *pc, size_t size) +{ + 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); + } +} + +PenroseCoords *penrose_coords_copy(PenroseCoords *pc_in) +{ + 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; +} + +/* + * 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) +{ + TransitionResult tr; + + penrosectx_extend_coords(ctx, pc, depth+2); + + /* Look up the transition out of the starting triangle */ + tr = transition(pc->c[depth+1], pc->c[depth], edge); + + /* 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) +{ + /* 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); +} + +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 */ + }; + + 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) +{ + 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)); + } +} + +void penrose_free(PenroseTriangle *tri) +{ + penrose_coords_free(tri->pc); + sfree(tri); +} + +static bool penrose_relative_probability(char c) +{ + /* 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; + } +} + +static char penrose_choose_random(const char *possibilities, random_state *rs) +{ + const char *p; + unsigned long value, limit = 0; + + 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 const char *penrose_starting_tiles(int which) +{ + return which == PENROSE_P2 ? "ABUV" : "CDXY"; +} + +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; + } +} + +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; + + 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; +} + +void penrosectx_cleanup(PenroseContext *ctx) +{ + if (ctx->must_free_rs) + random_free(ctx->rs); + penrose_coords_free(ctx->prototype); +} + +PenroseCoords *penrosectx_initial_coords(PenroseContext *ctx) +{ + return penrose_coords_copy(ctx->prototype); +} + +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; + } + + ctx->prototype->c[ctx->prototype->nc] = penrose_choose_random( + penrose_valid_parents(ctx->prototype->c[ctx->prototype->nc-1]), + ctx->rs); + ctx->prototype->nc++; + } + } + + penrose_coords_make_space(pc, n); + while (pc->nc < n) { + pc->c[pc->nc] = ctx->prototype->c[pc->nc]; + pc->nc++; + } +} + +static Point penrose_triangle_edge_0_length(char c) +{ + static const Point one = {{ 1, 0, 0, 0 }}; + static const Point phi = {{ 1, 0, 1, -1 }}; + static const Point invphi = {{ 0, 0, 1, -1 }}; + + 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; + + /* 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; + } +} + +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); + + /* Orient the triangle by deciding what vector edge #0 should traverse */ + edge0 = point_mul(edge0, point_rot(ctx->orientation)); + + /* Place the triangle at an arbitrary position, in that orientation */ + tri->pc = penrose_coords_copy(ctx->prototype); + penrose_place(tri, origin, edge0, 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); + + return tri; +} + +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; +} + +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 unsigned penrose_sibling_edge_index(char c) +{ + switch (c) { + case 'A': case 'U': return 2; + case 'B': case 'V': return 1; + default: return 0; + } +} + +void penrosectx_generate( + PenroseContext *ctx, + bool (*inbounds)(void *inboundsctx, + const PenroseTriangle *tri), void *inboundsctx, + void (*tile)(void *tilectx, const Point *vertices), void *tilectx) +{ + tree234 *placed = newtree234(penrose_cmp); + PenroseTriangle *qhead = NULL, *qtail = NULL; + + { + PenroseTriangle *tri = penrose_initial(ctx); + + add234(placed, tri); + + tri->next = NULL; + tri->reported = false; + + if (inbounds(inboundsctx, tri)) + qhead = qtail = tri; + } + + while (qhead) { + PenroseTriangle *tri = qhead; + unsigned edge; + unsigned sibling_edge = penrose_sibling_edge_index(tri->pc->c[0]); + + for (edge = 0; edge < 3; edge++) { + PenroseTriangle *new_tri, *found_tri; + + new_tri = penrose_adjacent(ctx, tri, edge, NULL); + + if (!inbounds(inboundsctx, new_tri)) { + penrose_free(new_tri); + continue; + } + + 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); + + tri->reported = true; + found_tri->reported = true; + } + + penrose_free(new_tri); + continue; + } + + add234(placed, new_tri); + qtail->next = new_tri; + qtail = new_tri; + new_tri->next = NULL; + new_tri->reported = false; + } + + qhead = qhead->next; + } + + { + PenroseTriangle *tri; + while ((tri = delpos234(placed, 0)) != NULL) + penrose_free(tri); + freetree234(placed); + } +} + +const char *penrose_tiling_params_invalid( + const struct PenrosePatchParams *params, int which) +{ + size_t i; + + if (params->ncoords == 0) + return "expected at least one coordinate"; + + 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"; + } + + return NULL; +} + +struct PenroseOutputCtx { + int xoff, yoff; + Coord xmin, xmax, ymin, ymax; + + penrose_tile_callback_fn external_cb; + void *external_cbctx; +}; + +static bool inbounds(void *vctx, const PenroseTriangle *tri) +{ + struct PenroseOutputCtx *octx = (struct PenroseOutputCtx *)vctx; + size_t i; + + for (i = 0; i < 3; i++) { + Coord x = point_x(tri->vertices[i]); + Coord y = point_y(tri->vertices[i]); + + 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 true; +} + +static void null_output_tile(void *vctx, const Point *vertices) +{ +} + +static void really_output_tile(void *vctx, const Point *vertices) +{ + struct PenroseOutputCtx *octx = (struct PenroseOutputCtx *)vctx; + size_t i; + int coords[16]; + + for (i = 0; i < 4; i++) { + Coord x = point_x(vertices[i]); + Coord y = point_y(vertices[i]); + + 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; + } + + octx->external_cb(octx->external_cbctx, coords); +} + +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; +} + +void penrose_tiling_randomise(struct PenrosePatchParams *params, int which, + int w, int h, random_state *rs) +{ + PenroseContext ctx[1]; + struct PenroseOutputCtx octx[1]; + + penrose_set_bounds(octx, w, h); + + penrosectx_init_random(ctx, rs, which); + penrosectx_generate(ctx, inbounds, octx, null_output_tile, NULL); + + 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); + + 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); } diff --git a/penrose.h b/penrose.h index 838364d..e6e927d 100644 --- 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 -#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 */ +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 -extern int penrose(penrose_state *state, int which, int angle); +bool penrose_valid_letter(char c, 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); +/* + * 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. + * + * 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. + */ +void penrose_tiling_randomise(struct PenrosePatchParams *params, int which, + int w, int h, random_state *rs); -/* 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); +/* + * 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); + +/* + * 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); + +#define PENROSE_NVERTICES 4 + +void penrose_tiling_generate( + const struct PenrosePatchParams *params, int w, int h, + penrose_tile_callback_fn cb, void *cbctx); #endif