chiark / gitweb /
Rework the preset menu system to permit submenus.
[sgt-puzzles.git] / penrose.c
1 /* penrose.c
2  *
3  * Penrose tile generator.
4  *
5  * Uses half-tile technique outlined on:
6  *
7  * http://tartarus.org/simon/20110412-penrose/penrose.xhtml
8  */
9
10 #include <assert.h>
11 #include <string.h>
12 #include <math.h>
13 #include <stdio.h>
14
15 #include "puzzles.h" /* for malloc routines, and PI */
16 #include "penrose.h"
17
18 /* -------------------------------------------------------
19  * 36-degree basis vector arithmetic routines.
20  */
21
22 /* Imagine drawing a
23  * ten-point 'clock face' like this:
24  *
25  *                     -E
26  *                 -D   |   A
27  *                   \  |  /
28  *              -C.   \ | /   ,B
29  *                 `-._\|/_,-'
30  *                 ,-' /|\ `-.
31  *              -B'   / | \   `C
32  *                   /  |  \
33  *                 -A   |   D
34  *                      E
35  *
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
42  * representation sets
43  *
44  *   A = (1,0,0,0)
45  *   B = (0,1,0,0)
46  *   C = (0,0,1,0)
47  *   D = (0,0,0,1)
48  *
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
55  *
56  *   E = (-1,1,-1,1)
57  *
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,
62  *
63  *   rotate_one_notch_clockwise(a,b,c,d) = (-d, d+a, -d+b, d+c)
64  *
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.
68  *
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:
75  *
76  *   divide_by_phi(vector) = rotate(vector, 2) - rotate(vector, 3)
77  *   multiply_by_phi(vector) = rotate(vector, 1) - rotate(vector, 4)
78  *
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:
83  *
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)
86  *
87  * and you can verify for yourself that those two operations are
88  * inverses of each other (as you'd hope!).
89  *
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.)
95  */
96
97 struct vector { int a, b, c, d; };
98
99 static vector v_origin(void)
100 {
101     vector v;
102     v.a = v.b = v.c = v.d = 0;
103     return v;
104 }
105
106 /* We start with a unit vector of B: this means we can easily
107  * draw an isoceles triangle centred on the X axis. */
108 #ifdef TEST_VECTORS
109
110 static vector v_unit(void)
111 {
112     vector v;
113
114     v.b = 1;
115     v.a = v.c = v.d = 0;
116     return v;
117 }
118
119 #endif
120
121 #define COS54 0.5877852
122 #define SIN54 0.8090169
123 #define COS18 0.9510565
124 #define SIN18 0.3090169
125
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)
129 {
130     return (vs[i].a + vs[i].d) * COS54 +
131            (vs[i].b + vs[i].c) * COS18;
132 }
133
134 double v_y(vector *vs, int i)
135 {
136     return (vs[i].a - vs[i].d) * SIN54 +
137            (vs[i].b - vs[i].c) * SIN18;
138
139 }
140
141 static vector v_trans(vector v, vector trans)
142 {
143     v.a += trans.a;
144     v.b += trans.b;
145     v.c += trans.c;
146     v.d += trans.d;
147     return v;
148 }
149
150 static vector v_rotate_36(vector v)
151 {
152     vector vv;
153     vv.a = -v.d;
154     vv.b =  v.d + v.a;
155     vv.c = -v.d + v.b;
156     vv.d =  v.d + v.c;
157     return vv;
158 }
159
160 static vector v_rotate(vector v, int ang)
161 {
162     int i;
163
164     assert((ang % 36) == 0);
165     while (ang < 0) ang += 360;
166     ang = 360-ang;
167     for (i = 0; i < (ang/36); i++)
168         v = v_rotate_36(v);
169     return v;
170 }
171
172 #ifdef TEST_VECTORS
173
174 static vector v_scale(vector v, int sc)
175 {
176     v.a *= sc;
177     v.b *= sc;
178     v.c *= sc;
179     v.d *= sc;
180     return v;
181 }
182
183 #endif
184
185 static vector v_growphi(vector v)
186 {
187     vector vv;
188     vv.a = v.a + v.b - v.d;
189     vv.b = v.c + v.d;
190     vv.c = v.a + v.b;
191     vv.d = v.c + v.d - v.a;
192     return vv;
193 }
194
195 static vector v_shrinkphi(vector v)
196 {
197     vector vv;
198     vv.a = v.b - v.d;
199     vv.b = v.c + v.d - v.b;
200     vv.c = v.a + v.b - v.c;
201     vv.d = v.c - v.a;
202     return vv;
203 }
204
205 #ifdef TEST_VECTORS
206
207 static const char *v_debug(vector v)
208 {
209     static char buf[255];
210     sprintf(buf,
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));
213     return buf;
214 }
215
216 #endif
217
218 /* -------------------------------------------------------
219  * Tiling routines.
220  */
221
222 static vector xform_coord(vector vo, int shrink, vector vtrans, int ang)
223 {
224     if (shrink < 0)
225         vo = v_shrinkphi(vo);
226     else if (shrink > 0)
227         vo = v_growphi(vo);
228
229     vo = v_rotate(vo, ang);
230     vo = v_trans(vo, vtrans);
231
232     return vo;
233 }
234
235
236 #define XFORM(n,o,s,a) vs[(n)] = xform_coord(v_edge, (s), vs[(o)], (a))
237
238 static int penrose_p2_small(penrose_state *state, int depth, int flip,
239                             vector v_orig, vector v_edge);
240
241 static int penrose_p2_large(penrose_state *state, int depth, int flip,
242                             vector v_orig, vector v_edge)
243 {
244     vector vv_orig, vv_edge;
245
246 #ifdef DEBUG_PENROSE
247     {
248         vector vs[3];
249         vs[0] = v_orig;
250         XFORM(1, 0, 0, 0);
251         XFORM(2, 0, 0, -36*flip);
252
253         state->new_tile(state, vs, 3, depth);
254     }
255 #endif
256
257     if (flip > 0) {
258         vector vs[4];
259
260         vs[0] = v_orig;
261         XFORM(1, 0, 0, -36);
262         XFORM(2, 0, 0, 0);
263         XFORM(3, 0, 0, 36);
264
265         state->new_tile(state, vs, 4, depth);
266     }
267     if (depth >= state->max_depth) return 0;
268
269     vv_orig = v_trans(v_orig, v_rotate(v_edge, -36*flip));
270     vv_edge = v_rotate(v_edge, 108*flip);
271
272     penrose_p2_small(state, depth+1, flip,
273                      v_orig, v_shrinkphi(v_edge));
274
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));
279
280     return 0;
281 }
282
283 static int penrose_p2_small(penrose_state *state, int depth, int flip,
284                             vector v_orig, vector v_edge)
285 {
286     vector vv_orig;
287
288 #ifdef DEBUG_PENROSE
289     {
290         vector vs[3];
291         vs[0] = v_orig;
292         XFORM(1, 0, 0, 0);
293         XFORM(2, 0, -1, -36*flip);
294
295         state->new_tile(state, vs, 3, depth);
296     }
297 #endif
298
299     if (flip > 0) {
300         vector vs[4];
301
302         vs[0] = v_orig;
303         XFORM(1, 0, 0, -72);
304         XFORM(2, 0, -1, -36);
305         XFORM(3, 0, 0, 0);
306
307         state->new_tile(state, vs, 4, depth);
308     }
309
310     if (depth >= state->max_depth) return 0;
311
312     vv_orig = v_trans(v_orig, v_edge);
313
314     penrose_p2_large(state, depth+1, -flip,
315                      v_orig, v_shrinkphi(v_rotate(v_edge, -36*flip)));
316
317     penrose_p2_small(state, depth+1, flip,
318                      vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
319
320     return 0;
321 }
322
323 static int penrose_p3_small(penrose_state *state, int depth, int flip,
324                             vector v_orig, vector v_edge);
325
326 static int penrose_p3_large(penrose_state *state, int depth, int flip,
327                             vector v_orig, vector v_edge)
328 {
329     vector vv_orig;
330
331 #ifdef DEBUG_PENROSE
332     {
333         vector vs[3];
334         vs[0] = v_orig;
335         XFORM(1, 0, 1, 0);
336         XFORM(2, 0, 0, -36*flip);
337
338         state->new_tile(state, vs, 3, depth);
339     }
340 #endif
341
342     if (flip > 0) {
343         vector vs[4];
344
345         vs[0] = v_orig;
346         XFORM(1, 0, 0, -36);
347         XFORM(2, 0, 1, 0);
348         XFORM(3, 0, 0, 36);
349
350         state->new_tile(state, vs, 4, depth);
351     }
352     if (depth >= state->max_depth) return 0;
353
354     vv_orig = v_trans(v_orig, v_edge);
355
356     penrose_p3_large(state, depth+1, -flip,
357                      vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
358
359     penrose_p3_small(state, depth+1, flip,
360                      vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
361
362     vv_orig = v_trans(v_orig, v_growphi(v_edge));
363
364     penrose_p3_large(state, depth+1, flip,
365                      vv_orig, v_shrinkphi(v_rotate(v_edge, -144*flip)));
366
367
368     return 0;
369 }
370
371 static int penrose_p3_small(penrose_state *state, int depth, int flip,
372                             vector v_orig, vector v_edge)
373 {
374     vector vv_orig;
375
376 #ifdef DEBUG_PENROSE
377     {
378         vector vs[3];
379         vs[0] = v_orig;
380         XFORM(1, 0, 0, 0);
381         XFORM(2, 0, 0, -36*flip);
382
383         state->new_tile(state, vs, 3, depth);
384     }
385 #endif
386
387     if (flip > 0) {
388         vector vs[4];
389
390         vs[0] = v_orig;
391         XFORM(1, 0, 0, -36);
392         XFORM(3, 0, 0, 0);
393         XFORM(2, 3, 0, -36);
394
395         state->new_tile(state, vs, 4, depth);
396     }
397     if (depth >= state->max_depth) return 0;
398
399     /* NB these two are identical to the first two of p3_large. */
400     vv_orig = v_trans(v_orig, v_edge);
401
402     penrose_p3_large(state, depth+1, -flip,
403                      vv_orig, v_shrinkphi(v_rotate(v_edge, 180)));
404
405     penrose_p3_small(state, depth+1, flip,
406                      vv_orig, v_shrinkphi(v_rotate(v_edge, -108*flip)));
407
408     return 0;
409 }
410
411 /* -------------------------------------------------------
412  * Utility routines.
413  */
414
415 double penrose_side_length(double start_size, int depth)
416 {
417   return start_size / pow(PHI, depth);
418 }
419
420 void penrose_count_tiles(int depth, int *nlarge, int *nsmall)
421 {
422   /* Steal sgt's fibonacci thingummy. */
423 }
424
425 /*
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.
433  *
434  * (later mail: this is an overestimate by about 5%)
435  */
436
437 int penrose(penrose_state *state, int which, int angle)
438 {
439     vector vo = v_origin();
440     vector vb = v_origin();
441
442     vo.b = vo.c = -state->start_size;
443     vo = v_shrinkphi(v_shrinkphi(vo));
444
445     vb.b = state->start_size;
446
447     vo = v_rotate(vo, angle);
448     vb = v_rotate(vb, angle);
449
450     if (which == PENROSE_P2)
451         return penrose_p2_large(state, 0, 1, vo, vb);
452     else
453         return penrose_p3_small(state, 0, 1, vo, vb);
454 }
455
456 /*
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.
463  *
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.
471  */
472 #define INCIRCLE_RADIUS 0.22426 /* phi^-3 less 5%: see above */
473
474 void penrose_calculate_size(int which, int tilesize, int w, int h,
475                             double *required_radius, int *start_size, int *depth)
476 {
477     double rradius, size;
478     int n = 0;
479
480     /*
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.
484      */
485     if (which == PENROSE_P2)
486         tilesize = tilesize * 3 / 2;
487     else
488         tilesize = tilesize * 5 / 4;
489
490     rradius = tilesize * 3.11 * sqrt((double)(w*w + h*h));
491     size = tilesize;
492
493     while ((size * INCIRCLE_RADIUS) < rradius) {
494         n++;
495         size = size * PHI;
496     }
497
498     *start_size = (int)size;
499     *depth = n;
500     *required_radius = rradius;
501 }
502
503 /* -------------------------------------------------------
504  * Test code.
505  */
506
507 #ifdef TEST_PENROSE
508
509 #include <stdio.h>
510 #include <string.h>
511
512 int show_recursion = 0;
513 int ntiles, nfinal;
514
515 int test_cb(penrose_state *state, vector *vs, int n, int depth)
516 {
517     int i, xoff = 0, yoff = 0;
518     double l = penrose_side_length(state->start_size, depth);
519     double rball = l / 10.0;
520     const char *col;
521
522     ntiles++;
523     if (state->max_depth == depth) {
524         col = n == 4 ? "black" : "green";
525         nfinal++;
526     } else {
527         if (!show_recursion)
528             return 0;
529         col = n == 4 ? "red" : "blue";
530     }
531     if (n != 4) yoff = state->start_size;
532
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);
537     }
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);
541
542     return 0;
543 }
544
545 void usage_exit(void)
546 {
547     fprintf(stderr, "Usage: penrose-test [--recursion] P2|P3 SIZE DEPTH\n");
548     exit(1);
549 }
550
551 int main(int argc, char *argv[])
552 {
553     penrose_state ps;
554     int which = 0;
555
556     while (--argc > 0) {
557         char *p = *++argv;
558         if (!strcmp(p, "-h") || !strcmp(p, "--help")) {
559             usage_exit();
560         } else if (!strcmp(p, "--recursion")) {
561             show_recursion = 1;
562         } else if (*p == '-') {
563             fprintf(stderr, "Unrecognised option '%s'\n", p);
564             exit(1);
565         } else {
566             break;
567         }
568     }
569
570     if (argc < 3) usage_exit();
571
572     if (strcmp(argv[0], "P2") == 0) which = PENROSE_P2;
573     else if (strcmp(argv[0], "P3") == 0) which = PENROSE_P3;
574     else usage_exit();
575
576     ps.start_size = atoi(argv[1]);
577     ps.max_depth = atoi(argv[2]);
578     ps.new_tile = test_cb;
579
580     ntiles = nfinal = 0;
581
582     printf("\
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\
586 \n\
587 <svg xmlns=\"http://www.w3.org/2000/svg\"\n\
588 xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n\n");
589
590     printf("<g>\n");
591     penrose(&ps, which);
592     printf("</g>\n");
593
594     printf("<!-- %d tiles and %d leaf tiles total -->\n",
595            ntiles, nfinal);
596
597     printf("</svg>");
598
599     return 0;
600 }
601
602 #endif
603
604 #ifdef TEST_VECTORS
605
606 static void dbgv(const char *msg, vector v)
607 {
608     printf("%s: %s\n", msg, v_debug(v));
609 }
610
611 int main(int argc, const char *argv[])
612 {
613    vector v = v_unit();
614
615    dbgv("unit vector", v);
616    v = v_rotate(v, 36);
617    dbgv("rotated 36", v);
618    v = v_scale(v, 2);
619    dbgv("scaled x2", v);
620    v = v_shrinkphi(v);
621    dbgv("shrunk phi", v);
622    v = v_rotate(v, -36);
623    dbgv("rotated -36", v);
624
625    return 0;
626 }
627
628 #endif
629 /* vim: set shiftwidth=4 tabstop=8: */