puzzles/penrose.c

/* penrose.c
 *
 * Penrose tile generator.
 *
 * Uses half-tile technique outlined on:
 *
 * http://tartarus.org/simon/20110412-penrose/penrose.xhtml
 */

#include <assert.h>
#include <string.h>
#ifdef NO_TGMATH_H
#  include <math.h>
#else
#  include <tgmath.h>
#endif
#include <stdio.h>

#include "puzzles.h" /* for malloc routines, and PI */
#include "penrose.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 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);
    }
#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.
 *
 * (later mail: this is an overestimate by about 5%)
 */

int penrose(penrose_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_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;
}