3 * Penrose tile generator.
5 * Uses half-tile technique outlined on:
7 * http://tartarus.org/simon/20110412-penrose/penrose.xhtml
15 #include "puzzles.h" /* for malloc routines, and PI */
18 /* -------------------------------------------------------
19 * 36-degree basis vector arithmetic routines.
23 * ten-point 'clock face' like this:
36 * In case the ASCII art isn't clear, those are supposed to be ten
37 * vectors of length 1, all sticking out from the origin at equal
38 * angular spacing (hence 36 degrees). Our basis vectors are A,B,C,D (I
39 * choose them to be symmetric about the x-axis so that the final
40 * translation into 2d coordinates will also be symmetric, which I
41 * think will avoid minor rounding uglinesses), so our vector
49 * The fifth vector E looks at first glance as if it needs to be
50 * another basis vector, but in fact it doesn't, because it can be
51 * represented in terms of the other four. Imagine starting from the
52 * origin and following the path -A, +B, -C, +D: you'll find you've
53 * traced four sides of a pentagram, and ended up one E-vector away
54 * from the origin. So we have
58 * This tells us that we can rotate any vector in this system by 36
59 * degrees: if we start with a*A + b*B + c*C + d*D, we want to end up
60 * with a*B + b*C + c*D + d*E, and we substitute our identity for E to
61 * turn that into a*B + b*C + c*D + d*(-A+B-C+D). In other words,
63 * rotate_one_notch_clockwise(a,b,c,d) = (-d, d+a, -d+b, d+c)
65 * and you can verify for yourself that applying that operation
66 * repeatedly starting with (1,0,0,0) cycles round ten vectors and
67 * comes back to where it started.
69 * The other operation that may be required is to construct vectors
70 * with lengths that are multiples of phi. That can be done by
71 * observing that the vector C-B is parallel to E and has length 1/phi,
72 * and the vector D-A is parallel to E and has length phi. So this
73 * tells us that given any vector, we can construct one which points in
74 * the same direction and is 1/phi or phi times its length, like this:
76 * divide_by_phi(vector) = rotate(vector, 2) - rotate(vector, 3)
77 * multiply_by_phi(vector) = rotate(vector, 1) - rotate(vector, 4)
79 * where rotate(vector, n) means applying the above
80 * rotate_one_notch_clockwise primitive n times. Expanding out the
81 * applications of rotate gives the following direct representation in
82 * terms of the vector coordinates:
84 * divide_by_phi(a,b,c,d) = (b-d, c+d-b, a+b-c, c-a)
85 * multiply_by_phi(a,b,c,d) = (a+b-d, c+d, a+b, c+d-a)
87 * and you can verify for yourself that those two operations are
88 * inverses of each other (as you'd hope!).
90 * Having done all of this, testing for equality between two vectors is
91 * a trivial matter of comparing the four integer coordinates. (Which
92 * it _wouldn't_ have been if we'd kept E as a fifth basis vector,
93 * because then (-1,1,-1,1,0) and (0,0,0,0,1) would have had to be
94 * considered identical. So leaving E out is vital.)
97 struct vector { int a, b, c, d; };
99 static vector v_origin(void)
102 v.a = v.b = v.c = v.d = 0;
106 /* We start with a unit vector of B: this means we can easily
107 * draw an isoceles triangle centred on the X axis. */
110 static vector v_unit(void)
121 #define COS54 0.5877852
122 #define SIN54 0.8090169
123 #define COS18 0.9510565
124 #define SIN18 0.3090169
126 /* These two are a bit rough-and-ready for now. Note that B/C are
127 * 18 degrees from the x-axis, and A/D are 54 degrees. */
128 double v_x(vector *vs, int i)
130 return (vs[i].a + vs[i].d) * COS54 +
131 (vs[i].b + vs[i].c) * COS18;
134 double v_y(vector *vs, int i)
136 return (vs[i].a - vs[i].d) * SIN54 +
137 (vs[i].b - vs[i].c) * SIN18;
141 static vector v_trans(vector v, vector trans)
150 static vector v_rotate_36(vector v)
160 static vector v_rotate(vector v, int ang)
164 assert((ang % 36) == 0);
165 while (ang < 0) ang += 360;
167 for (i = 0; i < (ang/36); i++)
174 static vector v_scale(vector v, int sc)
185 static vector v_growphi(vector v)
188 vv.a = v.a + v.b - v.d;
191 vv.d = v.c + v.d - v.a;
195 static vector v_shrinkphi(vector v)
199 vv.b = v.c + v.d - v.b;
200 vv.c = v.a + v.b - v.c;
207 static const char *v_debug(vector v)
209 static char buf[255];
211 "(%d,%d,%d,%d)[%2.2f,%2.2f]",
212 v.a, v.b, v.c, v.d, v_x(&v,0), v_y(&v,0));
218 /* -------------------------------------------------------
222 static vector xform_coord(vector vo, int shrink, vector vtrans, int ang)
225 vo = v_shrinkphi(vo);
229 vo = v_rotate(vo, ang);
230 vo = v_trans(vo, vtrans);
236 #define XFORM(n,o,s,a) vs[(n)] = xform_coord(v_edge, (s), vs[(o)], (a))
238 static int penrose_p2_small(penrose_state *state, int depth, int flip,
239 vector v_orig, vector v_edge);
241 static int penrose_p2_large(penrose_state *state, int depth, int flip,
242 vector v_orig, vector v_edge)
244 vector vv_orig, vv_edge;
251 XFORM(2, 0, 0, -36*flip);
253 state->new_tile(state, vs, 3, depth);
265 state->new_tile(state, vs, 4, depth);
267 if (depth >= state->max_depth) return 0;
269 vv_orig = v_trans(v_orig, v_rotate(v_edge, -36*flip));
270 vv_edge = v_rotate(v_edge, 108*flip);
272 penrose_p2_small(state, depth+1, flip,
273 v_orig, v_shrinkphi(v_edge));
275 penrose_p2_large(state, depth+1, flip,
276 vv_orig, v_shrinkphi(vv_edge));
277 penrose_p2_large(state, depth+1, -flip,
278 vv_orig, v_shrinkphi(vv_edge));
283 static int penrose_p2_small(penrose_state *state, int depth, int flip,
284 vector v_orig, vector v_edge)
293 XFORM(2, 0, -1, -36*flip);
295 state->new_tile(state, vs, 3, depth);
304 XFORM(2, 0, -1, -36);
307 state->new_tile(state, vs, 4, depth);
310 if (depth >= state->max_depth) return 0;
312 vv_orig = v_trans(v_orig, v_edge);
314 penrose_p2_large(state, depth+1, -flip,
315 v_orig, v_shrinkphi(v_rotate(v_edge, -36*flip)));
317 penrose_p2_small(state, depth+1, flip,
318 vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
323 static int penrose_p3_small(penrose_state *state, int depth, int flip,
324 vector v_orig, vector v_edge);
326 static int penrose_p3_large(penrose_state *state, int depth, int flip,
327 vector v_orig, vector v_edge)
336 XFORM(2, 0, 0, -36*flip);
338 state->new_tile(state, vs, 3, depth);
350 state->new_tile(state, vs, 4, depth);
352 if (depth >= state->max_depth) return 0;
354 vv_orig = v_trans(v_orig, v_edge);
356 penrose_p3_large(state, depth+1, -flip,
357 vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
359 penrose_p3_small(state, depth+1, flip,
360 vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
362 vv_orig = v_trans(v_orig, v_growphi(v_edge));
364 penrose_p3_large(state, depth+1, flip,
365 vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
371 static int penrose_p3_small(penrose_state *state, int depth, int flip,
372 vector v_orig, vector v_edge)
381 XFORM(2, 0, 0, -36*flip);
383 state->new_tile(state, vs, 3, depth);
395 state->new_tile(state, vs, 4, depth);
397 if (depth >= state->max_depth) return 0;
399 /* NB these two are identical to the first two of p3_large. */
400 vv_orig = v_trans(v_orig, v_edge);
402 penrose_p3_large(state, depth+1, -flip,
403 vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
405 penrose_p3_small(state, depth+1, flip,
406 vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
411 /* -------------------------------------------------------
415 double penrose_side_length(double start_size, int depth)
417 return start_size / pow(PHI, depth);
420 void penrose_count_tiles(int depth, int *nlarge, int *nsmall)
422 /* Steal sgt's fibonacci thingummy. */
426 * It turns out that an acute isosceles triangle with sides in ratio 1:phi:phi
427 * has an incentre which is conveniently 2*phi^-2 of the way from the apex to
428 * the base. Why's that convenient? Because: if we situate the incentre of the
429 * triangle at the origin, then we can place the apex at phi^-2 * (B+C), and
430 * the other two vertices at apex-B and apex-C respectively. So that's an acute
431 * triangle with its long sides of unit length, covering a circle about the
432 * origin of radius 1-(2*phi^-2), which is conveniently enough phi^-3.
434 * (later mail: this is an overestimate by about 5%)
437 int penrose(penrose_state *state, int which, int angle)
439 vector vo = v_origin();
440 vector vb = v_origin();
442 vo.b = vo.c = -state->start_size;
443 vo = v_shrinkphi(v_shrinkphi(vo));
445 vb.b = state->start_size;
447 vo = v_rotate(vo, angle);
448 vb = v_rotate(vb, angle);
450 if (which == PENROSE_P2)
451 return penrose_p2_large(state, 0, 1, vo, vb);
453 return penrose_p3_small(state, 0, 1, vo, vb);
457 * We're asked for a MxN grid, which just means a tiling fitting into roughly
458 * an MxN space in some kind of reasonable unit - say, the side length of the
459 * two-arrow edges of the tiles. By some reasoning in a previous email, that
460 * means we want to pick some subarea of a circle of radius 3.11*sqrt(M^2+N^2).
461 * To cover that circle, we need to subdivide a triangle large enough that it
462 * contains a circle of that radius.
464 * Hence: start with those three vectors marking triangle vertices, scale them
465 * all up by phi repeatedly until the radius of the inscribed circle gets
466 * bigger than the target, and then recurse into that triangle with the same
467 * recursion depth as the number of times you scaled up. That will give you
468 * tiles of unit side length, covering a circle big enough that if you randomly
469 * choose an orientation and coordinates within the circle, you'll be able to
470 * get any valid piece of Penrose tiling of size MxN.
472 #define INCIRCLE_RADIUS 0.22426 /* phi^-3 less 5%: see above */
474 void penrose_calculate_size(int which, int tilesize, int w, int h,
475 double *required_radius, int *start_size, int *depth)
477 double rradius, size;
481 * Fudge factor to scale P2 and P3 tilings differently. This
482 * doesn't seem to have much relevance to questions like the
483 * average number of tiles per unit area; it's just aesthetic.
485 if (which == PENROSE_P2)
486 tilesize = tilesize * 3 / 2;
488 tilesize = tilesize * 5 / 4;
490 rradius = tilesize * 3.11 * sqrt((double)(w*w + h*h));
493 while ((size * INCIRCLE_RADIUS) < rradius) {
498 *start_size = (int)size;
500 *required_radius = rradius;
503 /* -------------------------------------------------------
512 int show_recursion = 0;
515 int test_cb(penrose_state *state, vector *vs, int n, int depth)
517 int i, xoff = 0, yoff = 0;
518 double l = penrose_side_length(state->start_size, depth);
519 double rball = l / 10.0;
523 if (state->max_depth == depth) {
524 col = n == 4 ? "black" : "green";
529 col = n == 4 ? "red" : "blue";
531 if (n != 4) yoff = state->start_size;
533 printf("<polygon points=\"");
534 for (i = 0; i < n; i++) {
535 printf("%s%f,%f", (i == 0) ? "" : " ",
536 v_x(vs, i) + xoff, v_y(vs, i) + yoff);
538 printf("\" style=\"fill: %s; fill-opacity: 0.2; stroke: %s\" />\n", col, col);
539 printf("<ellipse cx=\"%f\" cy=\"%f\" rx=\"%f\" ry=\"%f\" fill=\"%s\" />",
540 v_x(vs, 0) + xoff, v_y(vs, 0) + yoff, rball, rball, col);
545 void usage_exit(void)
547 fprintf(stderr, "Usage: penrose-test [--recursion] P2|P3 SIZE DEPTH\n");
551 int main(int argc, char *argv[])
558 if (!strcmp(p, "-h") || !strcmp(p, "--help")) {
560 } else if (!strcmp(p, "--recursion")) {
562 } else if (*p == '-') {
563 fprintf(stderr, "Unrecognised option '%s'\n", p);
570 if (argc < 3) usage_exit();
572 if (strcmp(argv[0], "P2") == 0) which = PENROSE_P2;
573 else if (strcmp(argv[0], "P3") == 0) which = PENROSE_P3;
576 ps.start_size = atoi(argv[1]);
577 ps.max_depth = atoi(argv[2]);
578 ps.new_tile = test_cb;
583 <?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n\
584 <!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 20010904//EN\"\n\
585 \"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n\
587 <svg xmlns=\"http://www.w3.org/2000/svg\"\n\
588 xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n\n");
594 printf("<!-- %d tiles and %d leaf tiles total -->\n",
606 static void dbgv(const char *msg, vector v)
608 printf("%s: %s\n", msg, v_debug(v));
611 int main(int argc, const char *argv[])
615 dbgv("unit vector", v);
617 dbgv("rotated 36", v);
619 dbgv("scaled x2", v);
621 dbgv("shrunk phi", v);
622 v = v_rotate(v, -36);
623 dbgv("rotated -36", v);
629 /* vim: set shiftwidth=4 tabstop=8: */