grid.c

   1 /*
   2  * (c) Lambros Lambrou 2008
   3  *
   4  * Code for working with general grids, which can be any planar graph
   5  * with faces, edges and vertices (dots).  Includes generators for a few
   6  * types of grid, including square, hexagonal, triangular and others.
   7  */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <assert.h>
  13 #include <ctype.h>
  14 #include <math.h>
  15
  16 #include "puzzles.h"
  17 #include "tree234.h"
  18 #include "grid.h"
  19
  20 /* Debugging options */
  21
  22 /*
  23 #define DEBUG_GRID
  24 */
  25
  26 /* ----------------------------------------------------------------------
  27  * Deallocate or dereference a grid
  28  */
  29 void grid_free(grid *g)
  30 {
  31     assert(g->refcount);
  32
  33     g->refcount--;
  34     if (g->refcount == 0) {
  35         int i;
  36         for (i = 0; i < g->num_faces; i++) {
  37             sfree(g->faces[i].dots);
  38             sfree(g->faces[i].edges);
  39         }
  40         for (i = 0; i < g->num_dots; i++) {
  41             sfree(g->dots[i].faces);
  42             sfree(g->dots[i].edges);
  43         }
  44         sfree(g->faces);
  45         sfree(g->edges);
  46         sfree(g->dots);
  47         sfree(g);
  48     }
  49 }
  50
  51 /* Used by the other grid generators.  Create a brand new grid with nothing
  52  * initialised (all lists are NULL) */
  53 static grid *grid_new(void)
  54 {
  55     grid *g = snew(grid);
  56     g->faces = NULL;
  57     g->edges = NULL;
  58     g->dots = NULL;
  59     g->num_faces = g->num_edges = g->num_dots = 0;
  60     g->refcount = 1;
  61     g->lowest_x = g->lowest_y = g->highest_x = g->highest_y = 0;
  62     return g;
  63 }
  64
  65 /* Helper function to calculate perpendicular distance from
  66  * a point P to a line AB.  A and B mustn't be equal here.
  67  *
  68  * Well-known formula for area A of a triangle:
  69  *                             /  1   1   1 \
  70  * 2A = determinant of matrix  | px  ax  bx |
  71  *                             \ py  ay  by /
  72  *
  73  * Also well-known: 2A = base * height
  74  *                     = perpendicular distance * line-length.
  75  *
  76  * Combining gives: distance = determinant / line-length(a,b)
  77  */
  78 static double point_line_distance(long px, long py,
  79                                   long ax, long ay,
  80                                   long bx, long by)
  81 {
  82     long det = ax*by - bx*ay + bx*py - px*by + px*ay - ax*py;
  83     double len;
  84     det = max(det, -det);
  85     len = sqrt(SQ(ax - bx) + SQ(ay - by));
  86     return det / len;
  87 }
  88
  89 /* Determine nearest edge to where the user clicked.
  90  * (x, y) is the clicked location, converted to grid coordinates.
  91  * Returns the nearest edge, or NULL if no edge is reasonably
  92  * near the position.
  93  *
  94  * Just judging edges by perpendicular distance is not quite right -
  95  * the edge might be "off to one side". So we insist that the triangle
  96  * with (x,y) has acute angles at the edge's dots.
  97  *
  98  *     edge1
  99  *  *---------*------
 100  *            |
 101  *            |      *(x,y)
 102  *      edge2 |
 103  *            |   edge2 is OK, but edge1 is not, even though
 104  *            |   edge1 is perpendicularly closer to (x,y)
 105  *            *
 106  *
 107  */
 108 grid_edge *grid_nearest_edge(grid *g, int x, int y)
 109 {
 110     grid_edge *best_edge;
 111     double best_distance = 0;
 112     int i;
 113
 114     best_edge = NULL;
 115
 116     for (i = 0; i < g->num_edges; i++) {
 117         grid_edge *e = &g->edges[i];
 118         long e2; /* squared length of edge */
 119         long a2, b2; /* squared lengths of other sides */
 120         double dist;
 121
 122         /* See if edge e is eligible - the triangle must have acute angles
 123          * at the edge's dots.
 124          * Pythagoras formula h^2 = a^2 + b^2 detects right-angles,
 125          * so detect acute angles by testing for h^2 < a^2 + b^2 */
 126         e2 = SQ((long)e->dot1->x - (long)e->dot2->x) + SQ((long)e->dot1->y - (long)e->dot2->y);
 127         a2 = SQ((long)e->dot1->x - (long)x) + SQ((long)e->dot1->y - (long)y);
 128         b2 = SQ((long)e->dot2->x - (long)x) + SQ((long)e->dot2->y - (long)y);
 129         if (a2 >= e2 + b2) continue;
 130         if (b2 >= e2 + a2) continue;
 131
 132         /* e is eligible so far.  Now check the edge is reasonably close
 133          * to where the user clicked.  Don't want to toggle an edge if the
 134          * click was way off the grid.
 135          * There is room for experimentation here.  We could check the
 136          * perpendicular distance is within a certain fraction of the length
 137          * of the edge.  That amounts to testing a rectangular region around
 138          * the edge.
 139          * Alternatively, we could check that the angle at the point is obtuse.
 140          * That would amount to testing a circular region with the edge as
 141          * diameter. */
 142         dist = point_line_distance((long)x, (long)y,
 143                                    (long)e->dot1->x, (long)e->dot1->y,
 144                                    (long)e->dot2->x, (long)e->dot2->y);
 145         /* Is dist more than half edge length ? */
 146         if (4 * SQ(dist) > e2)
 147             continue;
 148
 149         if (best_edge == NULL || dist < best_distance) {
 150             best_edge = e;
 151             best_distance = dist;
 152         }
 153     }
 154     return best_edge;
 155 }
 156
 157 /* ----------------------------------------------------------------------
 158  * Grid generation
 159  */
 160
 161 #ifdef DEBUG_GRID
 162 /* Show the basic grid information, before doing grid_make_consistent */
 163 static void grid_print_basic(grid *g)
 164 {
 165     /* TODO: Maybe we should generate an SVG image of the dots and lines
 166      * of the grid here, before grid_make_consistent.
 167      * Would help with debugging grid generation. */
 168     int i;
 169     printf("--- Basic Grid Data ---\n");
 170     for (i = 0; i < g->num_faces; i++) {
 171         grid_face *f = g->faces + i;
 172         printf("Face %d: dots[", i);
 173         int j;
 174         for (j = 0; j < f->order; j++) {
 175             grid_dot *d = f->dots[j];
 176             printf("%s%d", j ? "," : "", (int)(d - g->dots));
 177         }
 178         printf("]\n");
 179     }
 180 }
 181 /* Show the derived grid information, computed by grid_make_consistent */
 182 static void grid_print_derived(grid *g)
 183 {
 184     /* edges */
 185     int i;
 186     printf("--- Derived Grid Data ---\n");
 187     for (i = 0; i < g->num_edges; i++) {
 188         grid_edge *e = g->edges + i;
 189         printf("Edge %d: dots[%d,%d] faces[%d,%d]\n",
 190             i, (int)(e->dot1 - g->dots), (int)(e->dot2 - g->dots),
 191             e->face1 ? (int)(e->face1 - g->faces) : -1,
 192             e->face2 ? (int)(e->face2 - g->faces) : -1);
 193     }
 194     /* faces */
 195     for (i = 0; i < g->num_faces; i++) {
 196         grid_face *f = g->faces + i;
 197         int j;
 198         printf("Face %d: faces[", i);
 199         for (j = 0; j < f->order; j++) {
 200             grid_edge *e = f->edges[j];
 201             grid_face *f2 = (e->face1 == f) ? e->face2 : e->face1;
 202             printf("%s%d", j ? "," : "", f2 ? (int)(f2 - g->faces) : -1);
 203         }
 204         printf("]\n");
 205     }
 206     /* dots */
 207     for (i = 0; i < g->num_dots; i++) {
 208         grid_dot *d = g->dots + i;
 209         int j;
 210         printf("Dot %d: dots[", i);
 211         for (j = 0; j < d->order; j++) {
 212             grid_edge *e = d->edges[j];
 213             grid_dot *d2 = (e->dot1 == d) ? e->dot2 : e->dot1;
 214             printf("%s%d", j ? "," : "", (int)(d2 - g->dots));
 215         }
 216         printf("] faces[");
 217         for (j = 0; j < d->order; j++) {
 218             grid_face *f = d->faces[j];
 219             printf("%s%d", j ? "," : "", f ? (int)(f - g->faces) : -1);
 220         }
 221         printf("]\n");
 222     }
 223 }
 224 #endif /* DEBUG_GRID */
 225
 226 /* Helper function for building incomplete-edges list in
 227  * grid_make_consistent() */
 228 static int grid_edge_bydots_cmpfn(void *v1, void *v2)
 229 {
 230     grid_edge *a = v1;
 231     grid_edge *b = v2;
 232     grid_dot *da, *db;
 233
 234     /* Pointer subtraction is valid here, because all dots point into the
 235      * same dot-list (g->dots).
 236      * Edges are not "normalised" - the 2 dots could be stored in any order,
 237      * so we need to take this into account when comparing edges. */
 238
 239     /* Compare first dots */
 240     da = (a->dot1 < a->dot2) ? a->dot1 : a->dot2;
 241     db = (b->dot1 < b->dot2) ? b->dot1 : b->dot2;
 242     if (da != db)
 243         return db - da;
 244     /* Compare last dots */
 245     da = (a->dot1 < a->dot2) ? a->dot2 : a->dot1;
 246     db = (b->dot1 < b->dot2) ? b->dot2 : b->dot1;
 247     if (da != db)
 248         return db - da;
 249
 250     return 0;
 251 }
 252
 253 /* Input: grid has its dots and faces initialised:
 254  * - dots have (optionally) x and y coordinates, but no edges or faces
 255  * (pointers are NULL).
 256  * - edges not initialised at all
 257  * - faces initialised and know which dots they have (but no edges yet).  The
 258  * dots around each face are assumed to be clockwise.
 259  *
 260  * Output: grid is complete and valid with all relationships defined.
 261  */
 262 static void grid_make_consistent(grid *g)
 263 {
 264     int i;
 265     tree234 *incomplete_edges;
 266     grid_edge *next_new_edge; /* Where new edge will go into g->edges */
 267
 268 #ifdef DEBUG_GRID
 269     grid_print_basic(g);
 270 #endif
 271
 272     /* ====== Stage 1 ======
 273      * Generate edges
 274      */
 275
 276     /* We know how many dots and faces there are, so we can find the exact
 277      * number of edges from Euler's polyhedral formula: F + V = E + 2 .
 278      * We use "-1", not "-2" here, because Euler's formula includes the
 279      * infinite face, which we don't count. */
 280     g->num_edges = g->num_faces + g->num_dots - 1;
 281     g->edges = snewn(g->num_edges, grid_edge);
 282     next_new_edge = g->edges;
 283
 284     /* Iterate over faces, and over each face's dots, generating edges as we
 285      * go.  As we find each new edge, we can immediately fill in the edge's
 286      * dots, but only one of the edge's faces.  Later on in the iteration, we
 287      * will find the same edge again (unless it's on the border), but we will
 288      * know the other face.
 289      * For efficiency, maintain a list of the incomplete edges, sorted by
 290      * their dots. */
 291     incomplete_edges = newtree234(grid_edge_bydots_cmpfn);
 292     for (i = 0; i < g->num_faces; i++) {
 293         grid_face *f = g->faces + i;
 294         int j;
 295         for (j = 0; j < f->order; j++) {
 296             grid_edge e; /* fake edge for searching */
 297             grid_edge *edge_found;
 298             int j2 = j + 1;
 299             if (j2 == f->order)
 300                 j2 = 0;
 301             e.dot1 = f->dots[j];
 302             e.dot2 = f->dots[j2];
 303             /* Use del234 instead of find234, because we always want to
 304              * remove the edge if found */
 305             edge_found = del234(incomplete_edges, &e);
 306             if (edge_found) {
 307                 /* This edge already added, so fill out missing face.
 308                  * Edge is already removed from incomplete_edges. */
 309                 edge_found->face2 = f;
 310             } else {
 311                 assert(next_new_edge - g->edges < g->num_edges);
 312                 next_new_edge->dot1 = e.dot1;
 313                 next_new_edge->dot2 = e.dot2;
 314                 next_new_edge->face1 = f;
 315                 next_new_edge->face2 = NULL; /* potentially infinite face */
 316                 add234(incomplete_edges, next_new_edge);
 317                 ++next_new_edge;
 318             }
 319         }
 320     }
 321     freetree234(incomplete_edges);
 322
 323     /* ====== Stage 2 ======
 324      * For each face, build its edge list.
 325      */
 326
 327     /* Allocate space for each edge list.  Can do this, because each face's
 328      * edge-list is the same size as its dot-list. */
 329     for (i = 0; i < g->num_faces; i++) {
 330         grid_face *f = g->faces + i;
 331         int j;
 332         f->edges = snewn(f->order, grid_edge*);
 333         /* Preload with NULLs, to help detect potential bugs. */
 334         for (j = 0; j < f->order; j++)
 335             f->edges[j] = NULL;
 336     }
 337
 338     /* Iterate over each edge, and over both its faces.  Add this edge to
 339      * the face's edge-list, after finding where it should go in the
 340      * sequence. */
 341     for (i = 0; i < g->num_edges; i++) {
 342         grid_edge *e = g->edges + i;
 343         int j;
 344         for (j = 0; j < 2; j++) {
 345             grid_face *f = j ? e->face2 : e->face1;
 346             int k, k2;
 347             if (f == NULL) continue;
 348             /* Find one of the dots around the face */
 349             for (k = 0; k < f->order; k++) {
 350                 if (f->dots[k] == e->dot1)
 351                     break; /* found dot1 */
 352             }
 353             assert(k != f->order); /* Must find the dot around this face */
 354
 355             /* Labelling scheme: as we walk clockwise around the face,
 356              * starting at dot0 (f->dots[0]), we hit:
 357              * (dot0), edge0, dot1, edge1, dot2,...
 358              *
 359              *     0
 360              *  0-----1
 361              *        |
 362              *        |1
 363              *        |
 364              *  3-----2
 365              *     2
 366              *
 367              * Therefore, edgeK joins dotK and dot{K+1}
 368              */
 369
 370             /* Around this face, either the next dot or the previous dot
 371              * must be e->dot2.  Otherwise the edge is wrong. */
 372             k2 = k + 1;
 373             if (k2 == f->order)
 374                 k2 = 0;
 375             if (f->dots[k2] == e->dot2) {
 376                 /* dot1(k) and dot2(k2) go clockwise around this face, so add
 377                  * this edge at position k (see diagram). */
 378                 assert(f->edges[k] == NULL);
 379                 f->edges[k] = e;
 380                 continue;
 381             }
 382             /* Try previous dot */
 383             k2 = k - 1;
 384             if (k2 == -1)
 385                 k2 = f->order - 1;
 386             if (f->dots[k2] == e->dot2) {
 387                 /* dot1(k) and dot2(k2) go anticlockwise around this face. */
 388                 assert(f->edges[k2] == NULL);
 389                 f->edges[k2] = e;
 390                 continue;
 391             }
 392             assert(!"Grid broken: bad edge-face relationship");
 393         }
 394     }
 395
 396     /* ====== Stage 3 ======
 397      * For each dot, build its edge-list and face-list.
 398      */
 399
 400     /* We don't know how many edges/faces go around each dot, so we can't
 401      * allocate the right space for these lists.  Pre-compute the sizes by
 402      * iterating over each edge and recording a tally against each dot. */
 403     for (i = 0; i < g->num_dots; i++) {
 404         g->dots[i].order = 0;
 405     }
 406     for (i = 0; i < g->num_edges; i++) {
 407         grid_edge *e = g->edges + i;
 408         ++(e->dot1->order);
 409         ++(e->dot2->order);
 410     }
 411     /* Now we have the sizes, pre-allocate the edge and face lists. */
 412     for (i = 0; i < g->num_dots; i++) {
 413         grid_dot *d = g->dots + i;
 414         int j;
 415         assert(d->order >= 2); /* sanity check */
 416         d->edges = snewn(d->order, grid_edge*);
 417         d->faces = snewn(d->order, grid_face*);
 418         for (j = 0; j < d->order; j++) {
 419             d->edges[j] = NULL;
 420             d->faces[j] = NULL;
 421         }
 422     }
 423     /* For each dot, need to find a face that touches it, so we can seed
 424      * the edge-face-edge-face process around each dot. */
 425     for (i = 0; i < g->num_faces; i++) {
 426         grid_face *f = g->faces + i;
 427         int j;
 428         for (j = 0; j < f->order; j++) {
 429             grid_dot *d = f->dots[j];
 430             d->faces[0] = f;
 431         }
 432     }
 433     /* Each dot now has a face in its first slot.  Generate the remaining
 434      * faces and edges around the dot, by searching both clockwise and
 435      * anticlockwise from the first face.  Need to do both directions,
 436      * because of the possibility of hitting the infinite face, which
 437      * blocks progress.  But there's only one such face, so we will
 438      * succeed in finding every edge and face this way. */
 439     for (i = 0; i < g->num_dots; i++) {
 440         grid_dot *d = g->dots + i;
 441         int current_face1 = 0; /* ascends clockwise */
 442         int current_face2 = 0; /* descends anticlockwise */
 443
 444         /* Labelling scheme: as we walk clockwise around the dot, starting
 445          * at face0 (d->faces[0]), we hit:
 446          * (face0), edge0, face1, edge1, face2,...
 447          *
 448          *       0
 449          *       |
 450          *    0  |  1
 451          *       |
 452          *  -----d-----1
 453          *       |
 454          *       |  2
 455          *       |
 456          *       2
 457          *
 458          * So, for example, face1 should be joined to edge0 and edge1,
 459          * and those edges should appear in an anticlockwise sense around
 460          * that face (see diagram). */
 461
 462         /* clockwise search */
 463         while (TRUE) {
 464             grid_face *f = d->faces[current_face1];
 465             grid_edge *e;
 466             int j;
 467             assert(f != NULL);
 468             /* find dot around this face */
 469             for (j = 0; j < f->order; j++) {
 470                 if (f->dots[j] == d)
 471                     break;
 472             }
 473             assert(j != f->order); /* must find dot */
 474
 475             /* Around f, required edge is anticlockwise from the dot.  See
 476              * the other labelling scheme higher up, for why we subtract 1
 477              * from j. */
 478             j--;
 479             if (j == -1)
 480                 j = f->order - 1;
 481             e = f->edges[j];
 482             d->edges[current_face1] = e; /* set edge */
 483             current_face1++;
 484             if (current_face1 == d->order)
 485                 break;
 486             else {
 487                 /* set face */
 488                 d->faces[current_face1] =
 489                     (e->face1 == f) ? e->face2 : e->face1;
 490                 if (d->faces[current_face1] == NULL)
 491                     break; /* cannot progress beyond infinite face */
 492             }
 493         }
 494         /* If the clockwise search made it all the way round, don't need to
 495          * bother with the anticlockwise search. */
 496         if (current_face1 == d->order)
 497             continue; /* this dot is complete, move on to next dot */
 498
 499         /* anticlockwise search */
 500         while (TRUE) {
 501             grid_face *f = d->faces[current_face2];
 502             grid_edge *e;
 503             int j;
 504             assert(f != NULL);
 505             /* find dot around this face */
 506             for (j = 0; j < f->order; j++) {
 507                 if (f->dots[j] == d)
 508                     break;
 509             }
 510             assert(j != f->order); /* must find dot */
 511
 512             /* Around f, required edge is clockwise from the dot. */
 513             e = f->edges[j];
 514
 515             current_face2--;
 516             if (current_face2 == -1)
 517                 current_face2 = d->order - 1;
 518             d->edges[current_face2] = e; /* set edge */
 519
 520             /* set face */
 521             if (current_face2 == current_face1)
 522                 break;
 523             d->faces[current_face2] =
 524                     (e->face1 == f) ? e->face2 : e->face1;
 525             /* There's only 1 infinite face, so we must get all the way
 526              * to current_face1 before we hit it. */
 527             assert(d->faces[current_face2]);
 528         }
 529     }
 530
 531     /* ====== Stage 4 ======
 532      * Compute other grid settings
 533      */
 534
 535     /* Bounding rectangle */
 536     for (i = 0; i < g->num_dots; i++) {
 537         grid_dot *d = g->dots + i;
 538         if (i == 0) {
 539             g->lowest_x = g->highest_x = d->x;
 540             g->lowest_y = g->highest_y = d->y;
 541         } else {
 542             g->lowest_x = min(g->lowest_x, d->x);
 543             g->highest_x = max(g->highest_x, d->x);
 544             g->lowest_y = min(g->lowest_y, d->y);
 545             g->highest_y = max(g->highest_y, d->y);
 546         }
 547     }
 548
 549 #ifdef DEBUG_GRID
 550     grid_print_derived(g);
 551 #endif
 552 }
 553
 554 /* Helpers for making grid-generation easier.  These functions are only
 555  * intended for use during grid generation. */
 556
 557 /* Comparison function for the (tree234) sorted dot list */
 558 static int grid_point_cmp_fn(void *v1, void *v2)
 559 {
 560     grid_dot *p1 = v1;
 561     grid_dot *p2 = v2;
 562     if (p1->y != p2->y)
 563         return p2->y - p1->y;
 564     else
 565         return p2->x - p1->x;
 566 }
 567 /* Add a new face to the grid, with its dot list allocated.
 568  * Assumes there's enough space allocated for the new face in grid->faces */
 569 static void grid_face_add_new(grid *g, int face_size)
 570 {
 571     int i;
 572     grid_face *new_face = g->faces + g->num_faces;
 573     new_face->order = face_size;
 574     new_face->dots = snewn(face_size, grid_dot*);
 575     for (i = 0; i < face_size; i++)
 576         new_face->dots[i] = NULL;
 577     new_face->edges = NULL;
 578     g->num_faces++;
 579 }
 580 /* Assumes dot list has enough space */
 581 static grid_dot *grid_dot_add_new(grid *g, int x, int y)
 582 {
 583     grid_dot *new_dot = g->dots + g->num_dots;
 584     new_dot->order = 0;
 585     new_dot->edges = NULL;
 586     new_dot->faces = NULL;
 587     new_dot->x = x;
 588     new_dot->y = y;
 589     g->num_dots++;
 590     return new_dot;
 591 }
 592 /* Retrieve a dot with these (x,y) coordinates.  Either return an existing dot
 593  * in the dot_list, or add a new dot to the grid (and the dot_list) and
 594  * return that.
 595  * Assumes g->dots has enough capacity allocated */
 596 static grid_dot *grid_get_dot(grid *g, tree234 *dot_list, int x, int y)
 597 {
 598     grid_dot test, *ret;
 599
 600     test.order = 0;
 601     test.edges = NULL;
 602     test.faces = NULL;
 603     test.x = x;
 604     test.y = y;
 605     ret = find234(dot_list, &test, NULL);
 606     if (ret)
 607         return ret;
 608
 609     ret = grid_dot_add_new(g, x, y);
 610     add234(dot_list, ret);
 611     return ret;
 612 }
 613
 614 /* Sets the last face of the grid to include this dot, at this position
 615  * around the face. Assumes num_faces is at least 1 (a new face has
 616  * previously been added, with the required number of dots allocated) */
 617 static void grid_face_set_dot(grid *g, grid_dot *d, int position)
 618 {
 619     grid_face *last_face = g->faces + g->num_faces - 1;
 620     last_face->dots[position] = d;
 621 }
 622
 623 /*
 624  * Helper routines for grid_find_incentre.
 625  */
 626 static int solve_2x2_matrix(double mx[4], double vin[2], double vout[2])
 627 {
 628     double inv[4];
 629     double det;
 630     det = (mx[0]*mx[3] - mx[1]*mx[2]);
 631     if (det == 0)
 632         return FALSE;
 633
 634     inv[0] = mx[3] / det;
 635     inv[1] = -mx[1] / det;
 636     inv[2] = -mx[2] / det;
 637     inv[3] = mx[0] / det;
 638
 639     vout[0] = inv[0]*vin[0] + inv[1]*vin[1];
 640     vout[1] = inv[2]*vin[0] + inv[3]*vin[1];
 641
 642     return TRUE;
 643 }
 644 static int solve_3x3_matrix(double mx[9], double vin[3], double vout[3])
 645 {
 646     double inv[9];
 647     double det;
 648
 649     det = (mx[0]*mx[4]*mx[8] + mx[1]*mx[5]*mx[6] + mx[2]*mx[3]*mx[7] -
 650            mx[0]*mx[5]*mx[7] - mx[1]*mx[3]*mx[8] - mx[2]*mx[4]*mx[6]);
 651     if (det == 0)
 652         return FALSE;
 653
 654     inv[0] = (mx[4]*mx[8] - mx[5]*mx[7]) / det;
 655     inv[1] = (mx[2]*mx[7] - mx[1]*mx[8]) / det;
 656     inv[2] = (mx[1]*mx[5] - mx[2]*mx[4]) / det;
 657     inv[3] = (mx[5]*mx[6] - mx[3]*mx[8]) / det;
 658     inv[4] = (mx[0]*mx[8] - mx[2]*mx[6]) / det;
 659     inv[5] = (mx[2]*mx[3] - mx[0]*mx[5]) / det;
 660     inv[6] = (mx[3]*mx[7] - mx[4]*mx[6]) / det;
 661     inv[7] = (mx[1]*mx[6] - mx[0]*mx[7]) / det;
 662     inv[8] = (mx[0]*mx[4] - mx[1]*mx[3]) / det;
 663
 664     vout[0] = inv[0]*vin[0] + inv[1]*vin[1] + inv[2]*vin[2];
 665     vout[1] = inv[3]*vin[0] + inv[4]*vin[1] + inv[5]*vin[2];
 666     vout[2] = inv[6]*vin[0] + inv[7]*vin[1] + inv[8]*vin[2];
 667
 668     return TRUE;
 669 }
 670
 671 void grid_find_incentre(grid_face *f)
 672 {
 673     double xbest, ybest, bestdist;
 674     int i, j, k, m;
 675     grid_dot *edgedot1[3], *edgedot2[3];
 676     grid_dot *dots[3];
 677     int nedges, ndots;
 678
 679     if (f->has_incentre)
 680         return;
 681
 682     /*
 683      * Find the point in the polygon with the maximum distance to any
 684      * edge or corner.
 685      *
 686      * Such a point must exist which is in contact with at least three
 687      * edges and/or vertices. (Proof: if it's only in contact with two
 688      * edges and/or vertices, it can't even be at a _local_ maximum -
 689      * any such circle can always be expanded in some direction.) So
 690      * we iterate through all 3-subsets of the combined set of edges
 691      * and vertices; for each subset we generate one or two candidate
 692      * points that might be the incentre, and then we vet each one to
 693      * see if it's inside the polygon and what its maximum radius is.
 694      *
 695      * (There's one case which this algorithm will get noticeably
 696      * wrong, and that's when a continuum of equally good answers
 697      * exists due to parallel edges. Consider a long thin rectangle,
 698      * for instance, or a parallelogram. This algorithm will pick a
 699      * point near one end, and choose the end arbitrarily; obviously a
 700      * nicer point to choose would be in the centre. To fix this I
 701      * would have to introduce a special-case system which detected
 702      * parallel edges in advance, set aside all candidate points
 703      * generated using both edges in a parallel pair, and generated
 704      * some additional candidate points half way between them. Also,
 705      * of course, I'd have to cope with rounding error making such a
 706      * point look worse than one of its endpoints. So I haven't done
 707      * this for the moment, and will cross it if necessary when I come
 708      * to it.)
 709      *
 710      * We don't actually iterate literally over _edges_, in the sense
 711      * of grid_edge structures. Instead, we fill in edgedot1[] and
 712      * edgedot2[] with a pair of dots adjacent in the face's list of
 713      * vertices. This ensures that we get the edges in consistent
 714      * orientation, which we could not do from the grid structure
 715      * alone. (A moment's consideration of an order-3 vertex should
 716      * make it clear that if a notional arrow was written on each
 717      * edge, _at least one_ of the three faces bordering that vertex
 718      * would have to have the two arrows tip-to-tip or tail-to-tail
 719      * rather than tip-to-tail.)
 720      */
 721     nedges = ndots = 0;
 722     bestdist = 0;
 723     xbest = ybest = 0;
 724
 725     for (i = 0; i+2 < 2*f->order; i++) {
 726         if (i < f->order) {
 727             edgedot1[nedges] = f->dots[i];
 728             edgedot2[nedges++] = f->dots[(i+1)%f->order];
 729         } else
 730             dots[ndots++] = f->dots[i - f->order];
 731
 732         for (j = i+1; j+1 < 2*f->order; j++) {
 733             if (j < f->order) {
 734                 edgedot1[nedges] = f->dots[j];
 735                 edgedot2[nedges++] = f->dots[(j+1)%f->order];
 736             } else
 737                 dots[ndots++] = f->dots[j - f->order];
 738
 739             for (k = j+1; k < 2*f->order; k++) {
 740                 double cx[2], cy[2];   /* candidate positions */
 741                 int cn = 0;            /* number of candidates */
 742
 743                 if (k < f->order) {
 744                     edgedot1[nedges] = f->dots[k];
 745                     edgedot2[nedges++] = f->dots[(k+1)%f->order];
 746                 } else
 747                     dots[ndots++] = f->dots[k - f->order];
 748
 749                 /*
 750                  * Find a point, or pair of points, equidistant from
 751                  * all the specified edges and/or vertices.
 752                  */
 753                 if (nedges == 3) {
 754                     /*
 755                      * Three edges. This is a linear matrix equation:
 756                      * each row of the matrix represents the fact that
 757                      * the point (x,y) we seek is at distance r from
 758                      * that edge, and we solve three of those
 759                      * simultaneously to obtain x,y,r. (We ignore r.)
 760                      */
 761                     double matrix[9], vector[3], vector2[3];
 762                     int m;
 763
 764                     for (m = 0; m < 3; m++) {
 765                         int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
 766                         int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
 767                         int dx = x2-x1, dy = y2-y1;
 768
 769                         /*
 770                          * ((x,y) - (x1,y1)) . (dy,-dx) = r |(dx,dy)|
 771                          *
 772                          * => x dy - y dx - r |(dx,dy)| = (x1 dy - y1 dx)
 773                          */
 774                         matrix[3*m+0] = dy;
 775                         matrix[3*m+1] = -dx;
 776                         matrix[3*m+2] = -sqrt((double)dx*dx+(double)dy*dy);
 777                         vector[m] = (double)x1*dy - (double)y1*dx;
 778                     }
 779
 780                     if (solve_3x3_matrix(matrix, vector, vector2)) {
 781                         cx[cn] = vector2[0];
 782                         cy[cn] = vector2[1];
 783                         cn++;
 784                     }
 785                 } else if (nedges == 2) {
 786                     /*
 787                      * Two edges and a dot. This will end up in a
 788                      * quadratic equation.
 789                      *
 790                      * First, look at the two edges. Having our point
 791                      * be some distance r from both of them gives rise
 792                      * to a pair of linear equations in x,y,r of the
 793                      * form
 794                      *
 795                      *   (x-x1) dy - (y-y1) dx = r sqrt(dx^2+dy^2)
 796                      *
 797                      * We eliminate r between those equations to give
 798                      * us a single linear equation in x,y describing
 799                      * the locus of points equidistant from both lines
 800                      * - i.e. the angle bisector.
 801                      *
 802                      * We then choose one of x,y to be a parameter t,
 803                      * and derive linear formulae for x,y,r in terms
 804                      * of t. This enables us to write down the
 805                      * circular equation (x-xd)^2+(y-yd)^2=r^2 as a
 806                      * quadratic in t; solving that and substituting
 807                      * in for x,y gives us two candidate points.
 808                      */
 809                     double eqs[2][4];  /* a,b,c,d : ax+by+cr=d */
 810                     double eq[3];      /* a,b,c: ax+by=c */
 811                     double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
 812                     double q[3];                /* a,b,c: at^2+bt+c=0 */
 813                     double disc;
 814
 815                     /* Find equations of the two input lines. */
 816                     for (m = 0; m < 2; m++) {
 817                         int x1 = edgedot1[m]->x, x2 = edgedot2[m]->x;
 818                         int y1 = edgedot1[m]->y, y2 = edgedot2[m]->y;
 819                         int dx = x2-x1, dy = y2-y1;
 820
 821                         eqs[m][0] = dy;
 822                         eqs[m][1] = -dx;
 823                         eqs[m][2] = -sqrt(dx*dx+dy*dy);
 824                         eqs[m][3] = x1*dy - y1*dx;
 825                     }
 826
 827                     /* Derive the angle bisector by eliminating r. */
 828                     eq[0] = eqs[0][0]*eqs[1][2] - eqs[1][0]*eqs[0][2];
 829                     eq[1] = eqs[0][1]*eqs[1][2] - eqs[1][1]*eqs[0][2];
 830                     eq[2] = eqs[0][3]*eqs[1][2] - eqs[1][3]*eqs[0][2];
 831
 832                     /* Parametrise x and y in terms of some t. */
 833                     if (abs(eq[0]) < abs(eq[1])) {
 834                         /* Parameter is x. */
 835                         xt[0] = 1; xt[1] = 0;
 836                         yt[0] = -eq[0]/eq[1]; yt[1] = eq[2]/eq[1];
 837                     } else {
 838                         /* Parameter is y. */
 839                         yt[0] = 1; yt[1] = 0;
 840                         xt[0] = -eq[1]/eq[0]; xt[1] = eq[2]/eq[0];
 841                     }
 842
 843                     /* Find a linear representation of r using eqs[0]. */
 844                     rt[0] = -(eqs[0][0]*xt[0] + eqs[0][1]*yt[0])/eqs[0][2];
 845                     rt[1] = (eqs[0][3] - eqs[0][0]*xt[1] -
 846                              eqs[0][1]*yt[1])/eqs[0][2];
 847
 848                     /* Construct the quadratic equation. */
 849                     q[0] = -rt[0]*rt[0];
 850                     q[1] = -2*rt[0]*rt[1];
 851                     q[2] = -rt[1]*rt[1];
 852                     q[0] += xt[0]*xt[0];
 853                     q[1] += 2*xt[0]*(xt[1]-dots[0]->x);
 854                     q[2] += (xt[1]-dots[0]->x)*(xt[1]-dots[0]->x);
 855                     q[0] += yt[0]*yt[0];
 856                     q[1] += 2*yt[0]*(yt[1]-dots[0]->y);
 857                     q[2] += (yt[1]-dots[0]->y)*(yt[1]-dots[0]->y);
 858
 859                     /* And solve it. */
 860                     disc = q[1]*q[1] - 4*q[0]*q[2];
 861                     if (disc >= 0) {
 862                         double t;
 863
 864                         disc = sqrt(disc);
 865
 866                         t = (-q[1] + disc) / (2*q[0]);
 867                         cx[cn] = xt[0]*t + xt[1];
 868                         cy[cn] = yt[0]*t + yt[1];
 869                         cn++;
 870
 871                         t = (-q[1] - disc) / (2*q[0]);
 872                         cx[cn] = xt[0]*t + xt[1];
 873                         cy[cn] = yt[0]*t + yt[1];
 874                         cn++;
 875                     }
 876                 } else if (nedges == 1) {
 877                     /*
 878                      * Two dots and an edge. This one's another
 879                      * quadratic equation.
 880                      *
 881                      * The point we want must lie on the perpendicular
 882                      * bisector of the two dots; that much is obvious.
 883                      * So we can construct a parametrisation of that
 884                      * bisecting line, giving linear formulae for x,y
 885                      * in terms of t. We can also express the distance
 886                      * from the edge as such a linear formula.
 887                      *
 888                      * Then we set that equal to the radius of the
 889                      * circle passing through the two points, which is
 890                      * a Pythagoras exercise; that gives rise to a
 891                      * quadratic in t, which we solve.
 892                      */
 893                     double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
 894                     double q[3];                /* a,b,c: at^2+bt+c=0 */
 895                     double disc;
 896                     double halfsep;
 897
 898                     /* Find parametric formulae for x,y. */
 899                     {
 900                         int x1 = dots[0]->x, x2 = dots[1]->x;
 901                         int y1 = dots[0]->y, y2 = dots[1]->y;
 902                         int dx = x2-x1, dy = y2-y1;
 903                         double d = sqrt((double)dx*dx + (double)dy*dy);
 904
 905                         xt[1] = (x1+x2)/2.0;
 906                         yt[1] = (y1+y2)/2.0;
 907                         /* It's convenient if we have t at standard scale. */
 908                         xt[0] = -dy/d;
 909                         yt[0] = dx/d;
 910
 911                         /* Also note down half the separation between
 912                          * the dots, for use in computing the circle radius. */
 913                         halfsep = 0.5*d;
 914                     }
 915
 916                     /* Find a parametric formula for r. */
 917                     {
 918                         int x1 = edgedot1[0]->x, x2 = edgedot2[0]->x;
 919                         int y1 = edgedot1[0]->y, y2 = edgedot2[0]->y;
 920                         int dx = x2-x1, dy = y2-y1;
 921                         double d = sqrt((double)dx*dx + (double)dy*dy);
 922                         rt[0] = (xt[0]*dy - yt[0]*dx) / d;
 923                         rt[1] = ((xt[1]-x1)*dy - (yt[1]-y1)*dx) / d;
 924                     }
 925
 926                     /* Construct the quadratic equation. */
 927                     q[0] = rt[0]*rt[0];
 928                     q[1] = 2*rt[0]*rt[1];
 929                     q[2] = rt[1]*rt[1];
 930                     q[0] -= 1;
 931                     q[2] -= halfsep*halfsep;
 932
 933                     /* And solve it. */
 934                     disc = q[1]*q[1] - 4*q[0]*q[2];
 935                     if (disc >= 0) {
 936                         double t;
 937
 938                         disc = sqrt(disc);
 939
 940                         t = (-q[1] + disc) / (2*q[0]);
 941                         cx[cn] = xt[0]*t + xt[1];
 942                         cy[cn] = yt[0]*t + yt[1];
 943                         cn++;
 944
 945                         t = (-q[1] - disc) / (2*q[0]);
 946                         cx[cn] = xt[0]*t + xt[1];
 947                         cy[cn] = yt[0]*t + yt[1];
 948                         cn++;
 949                     }
 950                 } else if (nedges == 0) {
 951                     /*
 952                      * Three dots. This is another linear matrix
 953                      * equation, this time with each row of the matrix
 954                      * representing the perpendicular bisector between
 955                      * two of the points. Of course we only need two
 956                      * such lines to find their intersection, so we
 957                      * need only solve a 2x2 matrix equation.
 958                      */
 959
 960                     double matrix[4], vector[2], vector2[2];
 961                     int m;
 962
 963                     for (m = 0; m < 2; m++) {
 964                         int x1 = dots[m]->x, x2 = dots[m+1]->x;
 965                         int y1 = dots[m]->y, y2 = dots[m+1]->y;
 966                         int dx = x2-x1, dy = y2-y1;
 967
 968                         /*
 969                          * ((x,y) - (x1,y1)) . (dx,dy) = 1/2 |(dx,dy)|^2
 970                          *
 971                          * => 2x dx + 2y dy = dx^2+dy^2 + (2 x1 dx + 2 y1 dy)
 972                          */
 973                         matrix[2*m+0] = 2*dx;
 974                         matrix[2*m+1] = 2*dy;
 975                         vector[m] = ((double)dx*dx + (double)dy*dy +
 976                                      2.0*x1*dx + 2.0*y1*dy);
 977                     }
 978
 979                     if (solve_2x2_matrix(matrix, vector, vector2)) {
 980                         cx[cn] = vector2[0];
 981                         cy[cn] = vector2[1];
 982                         cn++;
 983                     }
 984                 }
 985
 986                 /*
 987                  * Now go through our candidate points and see if any
 988                  * of them are better than what we've got so far.
 989                  */
 990                 for (m = 0; m < cn; m++) {
 991                     double x = cx[m], y = cy[m];
 992
 993                     /*
 994                      * First, disqualify the point if it's not inside
 995                      * the polygon, which we work out by counting the
 996                      * edges to the right of the point. (For
 997                      * tiebreaking purposes when edges start or end on
 998                      * our y-coordinate or go right through it, we
 999                      * consider our point to be offset by a small
1000                      * _positive_ epsilon in both the x- and
1001                      * y-direction.)
1002                      */
1003                     int e, in = 0;
1004                     for (e = 0; e < f->order; e++) {
1005                         int xs = f->edges[e]->dot1->x;
1006                         int xe = f->edges[e]->dot2->x;
1007                         int ys = f->edges[e]->dot1->y;
1008                         int ye = f->edges[e]->dot2->y;
1009                         if ((y >= ys && y < ye) || (y >= ye && y < ys)) {
1010                             /*
1011                              * The line goes past our y-position. Now we need
1012                              * to know if its x-coordinate when it does so is
1013                              * to our right.
1014                              *
1015                              * The x-coordinate in question is mathematically
1016                              * (y - ys) * (xe - xs) / (ye - ys), and we want
1017                              * to know whether (x - xs) >= that. Of course we
1018                              * avoid the division, so we can work in integers;
1019                              * to do this we must multiply both sides of the
1020                              * inequality by ye - ys, which means we must
1021                              * first check that's not negative.
1022                              */
1023                             int num = xe - xs, denom = ye - ys;
1024                             if (denom < 0) {
1025                                 num = -num;
1026                                 denom = -denom;
1027                             }
1028                             if ((x - xs) * denom >= (y - ys) * num)
1029                                 in ^= 1;
1030                         }
1031                     }
1032
1033                     if (in) {
1034                         double mindist = HUGE_VAL;
1035                         int e, d;
1036
1037                         /*
1038                          * This point is inside the polygon, so now we check
1039                          * its minimum distance to every edge and corner.
1040                          * First the corners ...
1041                          */
1042                         for (d = 0; d < f->order; d++) {
1043                             int xp = f->dots[d]->x;
1044                             int yp = f->dots[d]->y;
1045                             double dx = x - xp, dy = y - yp;
1046                             double dist = dx*dx + dy*dy;
1047                             if (mindist > dist)
1048                                 mindist = dist;
1049                         }
1050
1051                         /*
1052                          * ... and now also check the perpendicular distance
1053                          * to every edge, if the perpendicular lies between
1054                          * the edge's endpoints.
1055                          */
1056                         for (e = 0; e < f->order; e++) {
1057                             int xs = f->edges[e]->dot1->x;
1058                             int xe = f->edges[e]->dot2->x;
1059                             int ys = f->edges[e]->dot1->y;
1060                             int ye = f->edges[e]->dot2->y;
1061
1062                             /*
1063                              * If s and e are our endpoints, and p our
1064                              * candidate circle centre, the foot of a
1065                              * perpendicular from p to the line se lies
1066                              * between s and e if and only if (p-s).(e-s) lies
1067                              * strictly between 0 and (e-s).(e-s).
1068                              */
1069                             int edx = xe - xs, edy = ye - ys;
1070                             double pdx = x - xs, pdy = y - ys;
1071                             double pde = pdx * edx + pdy * edy;
1072                             long ede = (long)edx * edx + (long)edy * edy;
1073                             if (0 < pde && pde < ede) {
1074                                 /*
1075                                  * Yes, the nearest point on this edge is
1076                                  * closer than either endpoint, so we must
1077                                  * take it into account by measuring the
1078                                  * perpendicular distance to the edge and
1079                                  * checking its square against mindist.
1080                                  */
1081
1082                                 double pdre = pdx * edy - pdy * edx;
1083                                 double sqlen = pdre * pdre / ede;
1084
1085                                 if (mindist > sqlen)
1086                                     mindist = sqlen;
1087                             }
1088                         }
1089
1090                         /*
1091                          * Right. Now we know the biggest circle around this
1092                          * point, so we can check it against bestdist.
1093                          */
1094                         if (bestdist < mindist) {
1095                             bestdist = mindist;
1096                             xbest = x;
1097                             ybest = y;
1098                         }
1099                     }
1100                 }
1101
1102                 if (k < f->order)
1103                     nedges--;
1104                 else
1105                     ndots--;
1106             }
1107             if (j < f->order)
1108                 nedges--;
1109             else
1110                 ndots--;
1111         }
1112         if (i < f->order)
1113             nedges--;
1114         else
1115             ndots--;
1116     }
1117
1118     assert(bestdist > 0);
1119
1120     f->has_incentre = TRUE;
1121     f->ix = xbest + 0.5;               /* round to nearest */
1122     f->iy = ybest + 0.5;
1123 }
1124
1125 /* ------ Generate various types of grid ------ */
1126
1127 /* General method is to generate faces, by calculating their dot coordinates.
1128  * As new faces are added, we keep track of all the dots so we can tell when
1129  * a new face reuses an existing dot.  For example, two squares touching at an
1130  * edge would generate six unique dots: four dots from the first face, then
1131  * two additional dots for the second face, because we detect the other two
1132  * dots have already been taken up.  This list is stored in a tree234
1133  * called "points".  No extra memory-allocation needed here - we store the
1134  * actual grid_dot* pointers, which all point into the g->dots list.
1135  * For this reason, we have to calculate coordinates in such a way as to
1136  * eliminate any rounding errors, so we can detect when a dot on one
1137  * face precisely lands on a dot of a different face.  No floating-point
1138  * arithmetic here!
1139  */
1140
1141 grid *grid_new_square(int width, int height)
1142 {
1143     int x, y;
1144     /* Side length */
1145     int a = 20;
1146
1147     /* Upper bounds - don't have to be exact */
1148     int max_faces = width * height;
1149     int max_dots = (width + 1) * (height + 1);
1150
1151     tree234 *points;
1152
1153     grid *g = grid_new();
1154     g->tilesize = a;
1155     g->faces = snewn(max_faces, grid_face);
1156     g->dots = snewn(max_dots, grid_dot);
1157
1158     points = newtree234(grid_point_cmp_fn);
1159
1160     /* generate square faces */
1161     for (y = 0; y < height; y++) {
1162         for (x = 0; x < width; x++) {
1163             grid_dot *d;
1164             /* face position */
1165             int px = a * x;
1166             int py = a * y;
1167
1168             grid_face_add_new(g, 4);
1169             d = grid_get_dot(g, points, px, py);
1170             grid_face_set_dot(g, d, 0);
1171             d = grid_get_dot(g, points, px + a, py);
1172             grid_face_set_dot(g, d, 1);
1173             d = grid_get_dot(g, points, px + a, py + a);
1174             grid_face_set_dot(g, d, 2);
1175             d = grid_get_dot(g, points, px, py + a);
1176             grid_face_set_dot(g, d, 3);
1177         }
1178     }
1179
1180     freetree234(points);
1181     assert(g->num_faces <= max_faces);
1182     assert(g->num_dots <= max_dots);
1183
1184     grid_make_consistent(g);
1185     return g;
1186 }
1187
1188 grid *grid_new_honeycomb(int width, int height)
1189 {
1190     int x, y;
1191     /* Vector for side of hexagon - ratio is close to sqrt(3) */
1192     int a = 15;
1193     int b = 26;
1194
1195     /* Upper bounds - don't have to be exact */
1196     int max_faces = width * height;
1197     int max_dots = 2 * (width + 1) * (height + 1);
1198
1199     tree234 *points;
1200
1201     grid *g = grid_new();
1202     g->tilesize = 3 * a;
1203     g->faces = snewn(max_faces, grid_face);
1204     g->dots = snewn(max_dots, grid_dot);
1205
1206     points = newtree234(grid_point_cmp_fn);
1207
1208     /* generate hexagonal faces */
1209     for (y = 0; y < height; y++) {
1210         for (x = 0; x < width; x++) {
1211             grid_dot *d;
1212             /* face centre */
1213             int cx = 3 * a * x;
1214             int cy = 2 * b * y;
1215             if (x % 2)
1216                 cy += b;
1217             grid_face_add_new(g, 6);
1218
1219             d = grid_get_dot(g, points, cx - a, cy - b);
1220             grid_face_set_dot(g, d, 0);
1221             d = grid_get_dot(g, points, cx + a, cy - b);
1222             grid_face_set_dot(g, d, 1);
1223             d = grid_get_dot(g, points, cx + 2*a, cy);
1224             grid_face_set_dot(g, d, 2);
1225             d = grid_get_dot(g, points, cx + a, cy + b);
1226             grid_face_set_dot(g, d, 3);
1227             d = grid_get_dot(g, points, cx - a, cy + b);
1228             grid_face_set_dot(g, d, 4);
1229             d = grid_get_dot(g, points, cx - 2*a, cy);
1230             grid_face_set_dot(g, d, 5);
1231         }
1232     }
1233
1234     freetree234(points);
1235     assert(g->num_faces <= max_faces);
1236     assert(g->num_dots <= max_dots);
1237
1238     grid_make_consistent(g);
1239     return g;
1240 }
1241
1242 /* Doesn't use the previous method of generation, it pre-dates it!
1243  * A triangular grid is just about simple enough to do by "brute force" */
1244 grid *grid_new_triangular(int width, int height)
1245 {
1246     int x,y;
1247
1248     /* Vector for side of triangle - ratio is close to sqrt(3) */
1249     int vec_x = 15;
1250     int vec_y = 26;
1251
1252     int index;
1253
1254     /* convenient alias */
1255     int w = width + 1;
1256
1257     grid *g = grid_new();
1258     g->tilesize = 18; /* adjust to your taste */
1259
1260     g->num_faces = width * height * 2;
1261     g->num_dots = (width + 1) * (height + 1);
1262     g->faces = snewn(g->num_faces, grid_face);
1263     g->dots = snewn(g->num_dots, grid_dot);
1264
1265     /* generate dots */
1266     index = 0;
1267     for (y = 0; y <= height; y++) {
1268         for (x = 0; x <= width; x++) {
1269             grid_dot *d = g->dots + index;
1270             /* odd rows are offset to the right */
1271             d->order = 0;
1272             d->edges = NULL;
1273             d->faces = NULL;
1274             d->x = x * 2 * vec_x + ((y % 2) ? vec_x : 0);
1275             d->y = y * vec_y;
1276             index++;
1277         }
1278     }
1279
1280     /* generate faces */
1281     index = 0;
1282     for (y = 0; y < height; y++) {
1283         for (x = 0; x < width; x++) {
1284             /* initialise two faces for this (x,y) */
1285             grid_face *f1 = g->faces + index;
1286             grid_face *f2 = f1 + 1;
1287             f1->edges = NULL;
1288             f1->order = 3;
1289             f1->dots = snewn(f1->order, grid_dot*);
1290             f2->edges = NULL;
1291             f2->order = 3;
1292             f2->dots = snewn(f2->order, grid_dot*);
1293
1294             /* face descriptions depend on whether the row-number is
1295              * odd or even */
1296             if (y % 2) {
1297                 f1->dots[0] = g->dots + y       * w + x;
1298                 f1->dots[1] = g->dots + (y + 1) * w + x + 1;
1299                 f1->dots[2] = g->dots + (y + 1) * w + x;
1300                 f2->dots[0] = g->dots + y       * w + x;
1301                 f2->dots[1] = g->dots + y       * w + x + 1;
1302                 f2->dots[2] = g->dots + (y + 1) * w + x + 1;
1303             } else {
1304                 f1->dots[0] = g->dots + y       * w + x;
1305                 f1->dots[1] = g->dots + y       * w + x + 1;
1306                 f1->dots[2] = g->dots + (y + 1) * w + x;
1307                 f2->dots[0] = g->dots + y       * w + x + 1;
1308                 f2->dots[1] = g->dots + (y + 1) * w + x + 1;
1309                 f2->dots[2] = g->dots + (y + 1) * w + x;
1310             }
1311             index += 2;
1312         }
1313     }
1314
1315     grid_make_consistent(g);
1316     return g;
1317 }
1318
1319 grid *grid_new_snubsquare(int width, int height)
1320 {
1321     int x, y;
1322     /* Vector for side of triangle - ratio is close to sqrt(3) */
1323     int a = 15;
1324     int b = 26;
1325
1326     /* Upper bounds - don't have to be exact */
1327     int max_faces = 3 * width * height;
1328     int max_dots = 2 * (width + 1) * (height + 1);
1329
1330     tree234 *points;
1331
1332     grid *g = grid_new();
1333     g->tilesize = 18;
1334     g->faces = snewn(max_faces, grid_face);
1335     g->dots = snewn(max_dots, grid_dot);
1336
1337     points = newtree234(grid_point_cmp_fn);
1338
1339     for (y = 0; y < height; y++) {
1340         for (x = 0; x < width; x++) {
1341             grid_dot *d;
1342             /* face position */
1343             int px = (a + b) * x;
1344             int py = (a + b) * y;
1345
1346             /* generate square faces */
1347             grid_face_add_new(g, 4);
1348             if ((x + y) % 2) {
1349                 d = grid_get_dot(g, points, px + a, py);
1350                 grid_face_set_dot(g, d, 0);
1351                 d = grid_get_dot(g, points, px + a + b, py + a);
1352                 grid_face_set_dot(g, d, 1);
1353                 d = grid_get_dot(g, points, px + b, py + a + b);
1354                 grid_face_set_dot(g, d, 2);
1355                 d = grid_get_dot(g, points, px, py + b);
1356                 grid_face_set_dot(g, d, 3);
1357             } else {
1358                 d = grid_get_dot(g, points, px + b, py);
1359                 grid_face_set_dot(g, d, 0);
1360                 d = grid_get_dot(g, points, px + a + b, py + b);
1361                 grid_face_set_dot(g, d, 1);
1362                 d = grid_get_dot(g, points, px + a, py + a + b);
1363                 grid_face_set_dot(g, d, 2);
1364                 d = grid_get_dot(g, points, px, py + a);
1365                 grid_face_set_dot(g, d, 3);
1366             }
1367
1368             /* generate up/down triangles */
1369             if (x > 0) {
1370                 grid_face_add_new(g, 3);
1371                 if ((x + y) % 2) {
1372                     d = grid_get_dot(g, points, px + a, py);
1373                     grid_face_set_dot(g, d, 0);
1374                     d = grid_get_dot(g, points, px, py + b);
1375                     grid_face_set_dot(g, d, 1);
1376                     d = grid_get_dot(g, points, px - a, py);
1377                     grid_face_set_dot(g, d, 2);
1378                 } else {
1379                     d = grid_get_dot(g, points, px, py + a);
1380                     grid_face_set_dot(g, d, 0);
1381                     d = grid_get_dot(g, points, px + a, py + a + b);
1382                     grid_face_set_dot(g, d, 1);
1383                     d = grid_get_dot(g, points, px - a, py + a + b);
1384                     grid_face_set_dot(g, d, 2);
1385                 }
1386             }
1387
1388             /* generate left/right triangles */
1389             if (y > 0) {
1390                 grid_face_add_new(g, 3);
1391                 if ((x + y) % 2) {
1392                     d = grid_get_dot(g, points, px + a, py);
1393                     grid_face_set_dot(g, d, 0);
1394                     d = grid_get_dot(g, points, px + a + b, py - a);
1395                     grid_face_set_dot(g, d, 1);
1396                     d = grid_get_dot(g, points, px + a + b, py + a);
1397                     grid_face_set_dot(g, d, 2);
1398                 } else {
1399                     d = grid_get_dot(g, points, px, py - a);
1400                     grid_face_set_dot(g, d, 0);
1401                     d = grid_get_dot(g, points, px + b, py);
1402                     grid_face_set_dot(g, d, 1);
1403                     d = grid_get_dot(g, points, px, py + a);
1404                     grid_face_set_dot(g, d, 2);
1405                 }
1406             }
1407         }
1408     }
1409
1410     freetree234(points);
1411     assert(g->num_faces <= max_faces);
1412     assert(g->num_dots <= max_dots);
1413
1414     grid_make_consistent(g);
1415     return g;
1416 }
1417
1418 grid *grid_new_cairo(int width, int height)
1419 {
1420     int x, y;
1421     /* Vector for side of pentagon - ratio is close to (4+sqrt(7))/3 */
1422     int a = 14;
1423     int b = 31;
1424
1425     /* Upper bounds - don't have to be exact */
1426     int max_faces = 2 * width * height;
1427     int max_dots = 3 * (width + 1) * (height + 1);
1428
1429     tree234 *points;
1430
1431     grid *g = grid_new();
1432     g->tilesize = 40;
1433     g->faces = snewn(max_faces, grid_face);
1434     g->dots = snewn(max_dots, grid_dot);
1435
1436     points = newtree234(grid_point_cmp_fn);
1437
1438     for (y = 0; y < height; y++) {
1439         for (x = 0; x < width; x++) {
1440             grid_dot *d;
1441             /* cell position */
1442             int px = 2 * b * x;
1443             int py = 2 * b * y;
1444
1445             /* horizontal pentagons */
1446             if (y > 0) {
1447                 grid_face_add_new(g, 5);
1448                 if ((x + y) % 2) {
1449                     d = grid_get_dot(g, points, px + a, py - b);
1450                     grid_face_set_dot(g, d, 0);
1451                     d = grid_get_dot(g, points, px + 2*b - a, py - b);
1452                     grid_face_set_dot(g, d, 1);
1453                     d = grid_get_dot(g, points, px + 2*b, py);
1454                     grid_face_set_dot(g, d, 2);
1455                     d = grid_get_dot(g, points, px + b, py + a);
1456                     grid_face_set_dot(g, d, 3);
1457                     d = grid_get_dot(g, points, px, py);
1458                     grid_face_set_dot(g, d, 4);
1459                 } else {
1460                     d = grid_get_dot(g, points, px, py);
1461                     grid_face_set_dot(g, d, 0);
1462                     d = grid_get_dot(g, points, px + b, py - a);
1463                     grid_face_set_dot(g, d, 1);
1464                     d = grid_get_dot(g, points, px + 2*b, py);
1465                     grid_face_set_dot(g, d, 2);
1466                     d = grid_get_dot(g, points, px + 2*b - a, py + b);
1467                     grid_face_set_dot(g, d, 3);
1468                     d = grid_get_dot(g, points, px + a, py + b);
1469                     grid_face_set_dot(g, d, 4);
1470                 }
1471             }
1472             /* vertical pentagons */
1473             if (x > 0) {
1474                 grid_face_add_new(g, 5);
1475                 if ((x + y) % 2) {
1476                     d = grid_get_dot(g, points, px, py);
1477                     grid_face_set_dot(g, d, 0);
1478                     d = grid_get_dot(g, points, px + b, py + a);
1479                     grid_face_set_dot(g, d, 1);
1480                     d = grid_get_dot(g, points, px + b, py + 2*b - a);
1481                     grid_face_set_dot(g, d, 2);
1482                     d = grid_get_dot(g, points, px, py + 2*b);
1483                     grid_face_set_dot(g, d, 3);
1484                     d = grid_get_dot(g, points, px - a, py + b);
1485                     grid_face_set_dot(g, d, 4);
1486                 } else {
1487                     d = grid_get_dot(g, points, px, py);
1488                     grid_face_set_dot(g, d, 0);
1489                     d = grid_get_dot(g, points, px + a, py + b);
1490                     grid_face_set_dot(g, d, 1);
1491                     d = grid_get_dot(g, points, px, py + 2*b);
1492                     grid_face_set_dot(g, d, 2);
1493                     d = grid_get_dot(g, points, px - b, py + 2*b - a);
1494                     grid_face_set_dot(g, d, 3);
1495                     d = grid_get_dot(g, points, px - b, py + a);
1496                     grid_face_set_dot(g, d, 4);
1497                 }
1498             }
1499         }
1500     }
1501
1502     freetree234(points);
1503     assert(g->num_faces <= max_faces);
1504     assert(g->num_dots <= max_dots);
1505
1506     grid_make_consistent(g);
1507     return g;
1508 }
1509
1510 grid *grid_new_greathexagonal(int width, int height)
1511 {
1512     int x, y;
1513     /* Vector for side of triangle - ratio is close to sqrt(3) */
1514     int a = 15;
1515     int b = 26;
1516
1517     /* Upper bounds - don't have to be exact */
1518     int max_faces = 6 * (width + 1) * (height + 1);
1519     int max_dots = 6 * width * height;
1520
1521     tree234 *points;
1522
1523     grid *g = grid_new();
1524     g->tilesize = 18;
1525     g->faces = snewn(max_faces, grid_face);
1526     g->dots = snewn(max_dots, grid_dot);
1527
1528     points = newtree234(grid_point_cmp_fn);
1529
1530     for (y = 0; y < height; y++) {
1531         for (x = 0; x < width; x++) {
1532             grid_dot *d;
1533             /* centre of hexagon */
1534             int px = (3*a + b) * x;
1535             int py = (2*a + 2*b) * y;
1536             if (x % 2)
1537                 py += a + b;
1538
1539             /* hexagon */
1540             grid_face_add_new(g, 6);
1541             d = grid_get_dot(g, points, px - a, py - b);
1542             grid_face_set_dot(g, d, 0);
1543             d = grid_get_dot(g, points, px + a, py - b);
1544             grid_face_set_dot(g, d, 1);
1545             d = grid_get_dot(g, points, px + 2*a, py);
1546             grid_face_set_dot(g, d, 2);
1547             d = grid_get_dot(g, points, px + a, py + b);
1548             grid_face_set_dot(g, d, 3);
1549             d = grid_get_dot(g, points, px - a, py + b);
1550             grid_face_set_dot(g, d, 4);
1551             d = grid_get_dot(g, points, px - 2*a, py);
1552             grid_face_set_dot(g, d, 5);
1553
1554             /* square below hexagon */
1555             if (y < height - 1) {
1556                 grid_face_add_new(g, 4);
1557                 d = grid_get_dot(g, points, px - a, py + b);
1558                 grid_face_set_dot(g, d, 0);
1559                 d = grid_get_dot(g, points, px + a, py + b);
1560                 grid_face_set_dot(g, d, 1);
1561                 d = grid_get_dot(g, points, px + a, py + 2*a + b);
1562                 grid_face_set_dot(g, d, 2);
1563                 d = grid_get_dot(g, points, px - a, py + 2*a + b);
1564                 grid_face_set_dot(g, d, 3);
1565             }
1566
1567             /* square below right */
1568             if ((x < width - 1) && (((x % 2) == 0) || (y < height - 1))) {
1569                 grid_face_add_new(g, 4);
1570                 d = grid_get_dot(g, points, px + 2*a, py);
1571                 grid_face_set_dot(g, d, 0);
1572                 d = grid_get_dot(g, points, px + 2*a + b, py + a);
1573                 grid_face_set_dot(g, d, 1);
1574                 d = grid_get_dot(g, points, px + a + b, py + a + b);
1575                 grid_face_set_dot(g, d, 2);
1576                 d = grid_get_dot(g, points, px + a, py + b);
1577                 grid_face_set_dot(g, d, 3);
1578             }
1579
1580             /* square below left */
1581             if ((x > 0) && (((x % 2) == 0) || (y < height - 1))) {
1582                 grid_face_add_new(g, 4);
1583                 d = grid_get_dot(g, points, px - 2*a, py);
1584                 grid_face_set_dot(g, d, 0);
1585                 d = grid_get_dot(g, points, px - a, py + b);
1586                 grid_face_set_dot(g, d, 1);
1587                 d = grid_get_dot(g, points, px - a - b, py + a + b);
1588                 grid_face_set_dot(g, d, 2);
1589                 d = grid_get_dot(g, points, px - 2*a - b, py + a);
1590                 grid_face_set_dot(g, d, 3);
1591             }
1592
1593             /* Triangle below right */
1594             if ((x < width - 1) && (y < height - 1)) {
1595                 grid_face_add_new(g, 3);
1596                 d = grid_get_dot(g, points, px + a, py + b);
1597                 grid_face_set_dot(g, d, 0);
1598                 d = grid_get_dot(g, points, px + a + b, py + a + b);
1599                 grid_face_set_dot(g, d, 1);
1600                 d = grid_get_dot(g, points, px + a, py + 2*a + b);
1601                 grid_face_set_dot(g, d, 2);
1602             }
1603
1604             /* Triangle below left */
1605             if ((x > 0) && (y < height - 1)) {
1606                 grid_face_add_new(g, 3);
1607                 d = grid_get_dot(g, points, px - a, py + b);
1608                 grid_face_set_dot(g, d, 0);
1609                 d = grid_get_dot(g, points, px - a, py + 2*a + b);
1610                 grid_face_set_dot(g, d, 1);
1611                 d = grid_get_dot(g, points, px - a - b, py + a + b);
1612                 grid_face_set_dot(g, d, 2);
1613             }
1614         }
1615     }
1616
1617     freetree234(points);
1618     assert(g->num_faces <= max_faces);
1619     assert(g->num_dots <= max_dots);
1620
1621     grid_make_consistent(g);
1622     return g;
1623 }
1624
1625 grid *grid_new_octagonal(int width, int height)
1626 {
1627     int x, y;
1628     /* b/a approx sqrt(2) */
1629     int a = 29;
1630     int b = 41;
1631
1632     /* Upper bounds - don't have to be exact */
1633     int max_faces = 2 * width * height;
1634     int max_dots = 4 * (width + 1) * (height + 1);
1635
1636     tree234 *points;
1637
1638     grid *g = grid_new();
1639     g->tilesize = 40;
1640     g->faces = snewn(max_faces, grid_face);
1641     g->dots = snewn(max_dots, grid_dot);
1642
1643     points = newtree234(grid_point_cmp_fn);
1644
1645     for (y = 0; y < height; y++) {
1646         for (x = 0; x < width; x++) {
1647             grid_dot *d;
1648             /* cell position */
1649             int px = (2*a + b) * x;
1650             int py = (2*a + b) * y;
1651             /* octagon */
1652             grid_face_add_new(g, 8);
1653             d = grid_get_dot(g, points, px + a, py);
1654             grid_face_set_dot(g, d, 0);
1655             d = grid_get_dot(g, points, px + a + b, py);
1656             grid_face_set_dot(g, d, 1);
1657             d = grid_get_dot(g, points, px + 2*a + b, py + a);
1658             grid_face_set_dot(g, d, 2);
1659             d = grid_get_dot(g, points, px + 2*a + b, py + a + b);
1660             grid_face_set_dot(g, d, 3);
1661             d = grid_get_dot(g, points, px + a + b, py + 2*a + b);
1662             grid_face_set_dot(g, d, 4);
1663             d = grid_get_dot(g, points, px + a, py + 2*a + b);
1664             grid_face_set_dot(g, d, 5);
1665             d = grid_get_dot(g, points, px, py + a + b);
1666             grid_face_set_dot(g, d, 6);
1667             d = grid_get_dot(g, points, px, py + a);
1668             grid_face_set_dot(g, d, 7);
1669
1670             /* diamond */
1671             if ((x > 0) && (y > 0)) {
1672                 grid_face_add_new(g, 4);
1673                 d = grid_get_dot(g, points, px, py - a);
1674                 grid_face_set_dot(g, d, 0);
1675                 d = grid_get_dot(g, points, px + a, py);
1676                 grid_face_set_dot(g, d, 1);
1677                 d = grid_get_dot(g, points, px, py + a);
1678                 grid_face_set_dot(g, d, 2);
1679                 d = grid_get_dot(g, points, px - a, py);
1680                 grid_face_set_dot(g, d, 3);
1681             }
1682         }
1683     }
1684
1685     freetree234(points);
1686     assert(g->num_faces <= max_faces);
1687     assert(g->num_dots <= max_dots);
1688
1689     grid_make_consistent(g);
1690     return g;
1691 }
1692
1693 grid *grid_new_kites(int width, int height)
1694 {
1695     int x, y;
1696     /* b/a approx sqrt(3) */
1697     int a = 15;
1698     int b = 26;
1699
1700     /* Upper bounds - don't have to be exact */
1701     int max_faces = 6 * width * height;
1702     int max_dots = 6 * (width + 1) * (height + 1);
1703
1704     tree234 *points;
1705
1706     grid *g = grid_new();
1707     g->tilesize = 40;
1708     g->faces = snewn(max_faces, grid_face);
1709     g->dots = snewn(max_dots, grid_dot);
1710
1711     points = newtree234(grid_point_cmp_fn);
1712
1713     for (y = 0; y < height; y++) {
1714         for (x = 0; x < width; x++) {
1715             grid_dot *d;
1716             /* position of order-6 dot */
1717             int px = 4*b * x;
1718             int py = 6*a * y;
1719             if (y % 2)
1720                 px += 2*b;
1721
1722             /* kite pointing up-left */
1723             grid_face_add_new(g, 4);
1724             d = grid_get_dot(g, points, px, py);
1725             grid_face_set_dot(g, d, 0);
1726             d = grid_get_dot(g, points, px + 2*b, py);
1727             grid_face_set_dot(g, d, 1);
1728             d = grid_get_dot(g, points, px + 2*b, py + 2*a);
1729             grid_face_set_dot(g, d, 2);
1730             d = grid_get_dot(g, points, px + b, py + 3*a);
1731             grid_face_set_dot(g, d, 3);
1732
1733             /* kite pointing up */
1734             grid_face_add_new(g, 4);
1735             d = grid_get_dot(g, points, px, py);
1736             grid_face_set_dot(g, d, 0);
1737             d = grid_get_dot(g, points, px + b, py + 3*a);
1738             grid_face_set_dot(g, d, 1);
1739             d = grid_get_dot(g, points, px, py + 4*a);
1740             grid_face_set_dot(g, d, 2);
1741             d = grid_get_dot(g, points, px - b, py + 3*a);
1742             grid_face_set_dot(g, d, 3);
1743
1744             /* kite pointing up-right */
1745             grid_face_add_new(g, 4);
1746             d = grid_get_dot(g, points, px, py);
1747             grid_face_set_dot(g, d, 0);
1748             d = grid_get_dot(g, points, px - b, py + 3*a);
1749             grid_face_set_dot(g, d, 1);
1750             d = grid_get_dot(g, points, px - 2*b, py + 2*a);
1751             grid_face_set_dot(g, d, 2);
1752             d = grid_get_dot(g, points, px - 2*b, py);
1753             grid_face_set_dot(g, d, 3);
1754
1755             /* kite pointing down-right */
1756             grid_face_add_new(g, 4);
1757             d = grid_get_dot(g, points, px, py);
1758             grid_face_set_dot(g, d, 0);
1759             d = grid_get_dot(g, points, px - 2*b, py);
1760             grid_face_set_dot(g, d, 1);
1761             d = grid_get_dot(g, points, px - 2*b, py - 2*a);
1762             grid_face_set_dot(g, d, 2);
1763             d = grid_get_dot(g, points, px - b, py - 3*a);
1764             grid_face_set_dot(g, d, 3);
1765
1766             /* kite pointing down */
1767             grid_face_add_new(g, 4);
1768             d = grid_get_dot(g, points, px, py);
1769             grid_face_set_dot(g, d, 0);
1770             d = grid_get_dot(g, points, px - b, py - 3*a);
1771             grid_face_set_dot(g, d, 1);
1772             d = grid_get_dot(g, points, px, py - 4*a);
1773             grid_face_set_dot(g, d, 2);
1774             d = grid_get_dot(g, points, px + b, py - 3*a);
1775             grid_face_set_dot(g, d, 3);
1776
1777             /* kite pointing down-left */
1778             grid_face_add_new(g, 4);
1779             d = grid_get_dot(g, points, px, py);
1780             grid_face_set_dot(g, d, 0);
1781             d = grid_get_dot(g, points, px + b, py - 3*a);
1782             grid_face_set_dot(g, d, 1);
1783             d = grid_get_dot(g, points, px + 2*b, py - 2*a);
1784             grid_face_set_dot(g, d, 2);
1785             d = grid_get_dot(g, points, px + 2*b, py);
1786             grid_face_set_dot(g, d, 3);
1787         }
1788     }
1789
1790     freetree234(points);
1791     assert(g->num_faces <= max_faces);
1792     assert(g->num_dots <= max_dots);
1793
1794     grid_make_consistent(g);
1795     return g;
1796 }
1797
1798 grid *grid_new_floret(int width, int height)
1799 {
1800     int x, y;
1801     /* Vectors for sides; weird numbers needed to keep puzzle aligned with window
1802      * -py/px is close to tan(30 - atan(sqrt(3)/9))
1803      * using py=26 makes everything lean to the left, rather than right
1804      */
1805     int px = 75, py = -26;       /* |( 75, -26)| = 79.43 */
1806     int qx = 4*px/5, qy = -py*2; /* |( 60,  52)| = 79.40 */
1807     int rx = qx-px, ry = qy-py;  /* |(-15,  78)| = 79.38 */
1808
1809     /* Upper bounds - don't have to be exact */
1810     int max_faces = 6 * width * height;
1811     int max_dots = 9 * (width + 1) * (height + 1);
1812
1813     tree234 *points;
1814
1815     grid *g = grid_new();
1816     g->tilesize = 2 * px;
1817     g->faces = snewn(max_faces, grid_face);
1818     g->dots = snewn(max_dots, grid_dot);
1819
1820     points = newtree234(grid_point_cmp_fn);
1821
1822     /* generate pentagonal faces */
1823     for (y = 0; y < height; y++) {
1824         for (x = 0; x < width; x++) {
1825             grid_dot *d;
1826             /* face centre */
1827             int cx = (6*px+3*qx)/2 * x;
1828             int cy = (4*py-5*qy) * y;
1829             if (x % 2)
1830                 cy -= (4*py-5*qy)/2;
1831             else if (y && y == height-1)
1832                 continue; /* make better looking grids?  try 3x3 for instance */
1833
1834             grid_face_add_new(g, 5);
1835             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
1836             d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 1);
1837             d = grid_get_dot(g, points, cx+2*rx+qx, cy+2*ry+qy); grid_face_set_dot(g, d, 2);
1838             d = grid_get_dot(g, points, cx+2*qx+rx, cy+2*qy+ry); grid_face_set_dot(g, d, 3);
1839             d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 4);
1840
1841             grid_face_add_new(g, 5);
1842             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
1843             d = grid_get_dot(g, points, cx+2*qx   , cy+2*qy   ); grid_face_set_dot(g, d, 1);
1844             d = grid_get_dot(g, points, cx+2*qx+px, cy+2*qy+py); grid_face_set_dot(g, d, 2);
1845             d = grid_get_dot(g, points, cx+2*px+qx, cy+2*py+qy); grid_face_set_dot(g, d, 3);
1846             d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 4);
1847
1848             grid_face_add_new(g, 5);
1849             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
1850             d = grid_get_dot(g, points, cx+2*px   , cy+2*py   ); grid_face_set_dot(g, d, 1);
1851             d = grid_get_dot(g, points, cx+2*px-rx, cy+2*py-ry); grid_face_set_dot(g, d, 2);
1852             d = grid_get_dot(g, points, cx-2*rx+px, cy-2*ry+py); grid_face_set_dot(g, d, 3);
1853             d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 4);
1854
1855             grid_face_add_new(g, 5);
1856             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
1857             d = grid_get_dot(g, points, cx-2*rx   , cy-2*ry   ); grid_face_set_dot(g, d, 1);
1858             d = grid_get_dot(g, points, cx-2*rx-qx, cy-2*ry-qy); grid_face_set_dot(g, d, 2);
1859             d = grid_get_dot(g, points, cx-2*qx-rx, cy-2*qy-ry); grid_face_set_dot(g, d, 3);
1860             d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 4);
1861
1862             grid_face_add_new(g, 5);
1863             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
1864             d = grid_get_dot(g, points, cx-2*qx   , cy-2*qy   ); grid_face_set_dot(g, d, 1);
1865             d = grid_get_dot(g, points, cx-2*qx-px, cy-2*qy-py); grid_face_set_dot(g, d, 2);
1866             d = grid_get_dot(g, points, cx-2*px-qx, cy-2*py-qy); grid_face_set_dot(g, d, 3);
1867             d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 4);
1868
1869             grid_face_add_new(g, 5);
1870             d = grid_get_dot(g, points, cx        , cy        ); grid_face_set_dot(g, d, 0);
1871             d = grid_get_dot(g, points, cx-2*px   , cy-2*py   ); grid_face_set_dot(g, d, 1);
1872             d = grid_get_dot(g, points, cx-2*px+rx, cy-2*py+ry); grid_face_set_dot(g, d, 2);
1873             d = grid_get_dot(g, points, cx+2*rx-px, cy+2*ry-py); grid_face_set_dot(g, d, 3);
1874             d = grid_get_dot(g, points, cx+2*rx   , cy+2*ry   ); grid_face_set_dot(g, d, 4);
1875         }
1876     }
1877
1878     freetree234(points);
1879     assert(g->num_faces <= max_faces);
1880     assert(g->num_dots <= max_dots);
1881
1882     grid_make_consistent(g);
1883     return g;
1884 }
1885
1886 grid *grid_new_dodecagonal(int width, int height)
1887 {
1888     int x, y;
1889     /* Vector for side of triangle - ratio is close to sqrt(3) */
1890     int a = 15;
1891     int b = 26;
1892
1893     /* Upper bounds - don't have to be exact */
1894     int max_faces = 3 * width * height;
1895     int max_dots = 14 * width * height;
1896
1897     tree234 *points;
1898
1899     grid *g = grid_new();
1900     g->tilesize = b;
1901     g->faces = snewn(max_faces, grid_face);
1902     g->dots = snewn(max_dots, grid_dot);
1903
1904     points = newtree234(grid_point_cmp_fn);
1905
1906     for (y = 0; y < height; y++) {
1907         for (x = 0; x < width; x++) {
1908             grid_dot *d;
1909             /* centre of dodecagon */
1910             int px = (4*a + 2*b) * x;
1911             int py = (3*a + 2*b) * y;
1912             if (y % 2)
1913                 px += 2*a + b;
1914
1915             /* dodecagon */
1916             grid_face_add_new(g, 12);
1917             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
1918             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
1919             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
1920             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
1921             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
1922             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
1923             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
1924             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
1925             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
1926             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
1927             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
1928             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
1929
1930             /* triangle below dodecagon */
1931             if ((y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
1932                 grid_face_add_new(g, 3);
1933                 d = grid_get_dot(g, points, px + a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
1934                 d = grid_get_dot(g, points, px    , py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
1935                 d = grid_get_dot(g, points, px - a, py + (2*a +   b)); grid_face_set_dot(g, d, 2);
1936             }
1937
1938             /* triangle above dodecagon */
1939             if ((y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2)))) {
1940                 grid_face_add_new(g, 3);
1941                 d = grid_get_dot(g, points, px - a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
1942                 d = grid_get_dot(g, points, px    , py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
1943                 d = grid_get_dot(g, points, px + a, py - (2*a +   b)); grid_face_set_dot(g, d, 2);
1944             }
1945         }
1946     }
1947
1948     freetree234(points);
1949     assert(g->num_faces <= max_faces);
1950     assert(g->num_dots <= max_dots);
1951
1952     grid_make_consistent(g);
1953     return g;
1954 }
1955
1956 grid *grid_new_greatdodecagonal(int width, int height)
1957 {
1958     int x, y;
1959     /* Vector for side of triangle - ratio is close to sqrt(3) */
1960     int a = 15;
1961     int b = 26;
1962
1963     /* Upper bounds - don't have to be exact */
1964     int max_faces = 30 * width * height;
1965     int max_dots = 200 * width * height;
1966
1967     tree234 *points;
1968
1969     grid *g = grid_new();
1970     g->tilesize = b;
1971     g->faces = snewn(max_faces, grid_face);
1972     g->dots = snewn(max_dots, grid_dot);
1973
1974     points = newtree234(grid_point_cmp_fn);
1975
1976     for (y = 0; y < height; y++) {
1977         for (x = 0; x < width; x++) {
1978             grid_dot *d;
1979             /* centre of dodecagon */
1980             int px = (6*a + 2*b) * x;
1981             int py = (3*a + 3*b) * y;
1982             if (y % 2)
1983                 px += 3*a + b;
1984
1985             /* dodecagon */
1986             grid_face_add_new(g, 12);
1987             d = grid_get_dot(g, points, px + (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 0);
1988             d = grid_get_dot(g, points, px + (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 1);
1989             d = grid_get_dot(g, points, px + (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 2);
1990             d = grid_get_dot(g, points, px + (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 3);
1991             d = grid_get_dot(g, points, px + (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 4);
1992             d = grid_get_dot(g, points, px + (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 5);
1993             d = grid_get_dot(g, points, px - (  a    ), py + (2*a + b)); grid_face_set_dot(g, d, 6);
1994             d = grid_get_dot(g, points, px - (  a + b), py + (  a + b)); grid_face_set_dot(g, d, 7);
1995             d = grid_get_dot(g, points, px - (2*a + b), py + (  a    )); grid_face_set_dot(g, d, 8);
1996             d = grid_get_dot(g, points, px - (2*a + b), py - (  a    )); grid_face_set_dot(g, d, 9);
1997             d = grid_get_dot(g, points, px - (  a + b), py - (  a + b)); grid_face_set_dot(g, d, 10);
1998             d = grid_get_dot(g, points, px - (  a    ), py - (2*a + b)); grid_face_set_dot(g, d, 11);
1999
2000             /* hexagon below dodecagon */
2001             if (y < height - 1 && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2002                 grid_face_add_new(g, 6);
2003                 d = grid_get_dot(g, points, px +   a, py + (2*a +   b)); grid_face_set_dot(g, d, 0);
2004                 d = grid_get_dot(g, points, px + 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2005                 d = grid_get_dot(g, points, px +   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 2);
2006                 d = grid_get_dot(g, points, px -   a, py + (2*a + 3*b)); grid_face_set_dot(g, d, 3);
2007                 d = grid_get_dot(g, points, px - 2*a, py + (2*a + 2*b)); grid_face_set_dot(g, d, 4);
2008                 d = grid_get_dot(g, points, px -   a, py + (2*a +   b)); grid_face_set_dot(g, d, 5);
2009             }
2010
2011             /* hexagon above dodecagon */
2012             if (y && (x < width - 1 || !(y % 2)) && (x > 0 || (y % 2))) {
2013                 grid_face_add_new(g, 6);
2014                 d = grid_get_dot(g, points, px -   a, py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2015                 d = grid_get_dot(g, points, px - 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2016                 d = grid_get_dot(g, points, px -   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 2);
2017                 d = grid_get_dot(g, points, px +   a, py - (2*a + 3*b)); grid_face_set_dot(g, d, 3);
2018                 d = grid_get_dot(g, points, px + 2*a, py - (2*a + 2*b)); grid_face_set_dot(g, d, 4);
2019                 d = grid_get_dot(g, points, px +   a, py - (2*a +   b)); grid_face_set_dot(g, d, 5);
2020             }
2021
2022             /* square on right of dodecagon */
2023             if (x < width - 1) {
2024                 grid_face_add_new(g, 4);
2025                 d = grid_get_dot(g, points, px + 2*a + b, py - a); grid_face_set_dot(g, d, 0);
2026                 d = grid_get_dot(g, points, px + 4*a + b, py - a); grid_face_set_dot(g, d, 1);
2027                 d = grid_get_dot(g, points, px + 4*a + b, py + a); grid_face_set_dot(g, d, 2);
2028                 d = grid_get_dot(g, points, px + 2*a + b, py + a); grid_face_set_dot(g, d, 3);
2029             }
2030
2031             /* square on top right of dodecagon */
2032             if (y && (x < width - 1 || !(y % 2))) {
2033                 grid_face_add_new(g, 4);
2034                 d = grid_get_dot(g, points, px + (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 0);
2035                 d = grid_get_dot(g, points, px + (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 1);
2036                 d = grid_get_dot(g, points, px + (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 2);
2037                 d = grid_get_dot(g, points, px + (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 3);
2038             }
2039
2040             /* square on top left of dodecagon */
2041             if (y && (x || (y % 2))) {
2042                 grid_face_add_new(g, 4);
2043                 d = grid_get_dot(g, points, px - (  a + b), py - (  a +   b)); grid_face_set_dot(g, d, 0);
2044                 d = grid_get_dot(g, points, px - (2*a + b), py - (  a + 2*b)); grid_face_set_dot(g, d, 1);
2045                 d = grid_get_dot(g, points, px - (2*a    ), py - (2*a + 2*b)); grid_face_set_dot(g, d, 2);
2046                 d = grid_get_dot(g, points, px - (  a    ), py - (2*a +   b)); grid_face_set_dot(g, d, 3);
2047             }
2048         }
2049     }
2050
2051     freetree234(points);
2052     assert(g->num_faces <= max_faces);
2053     assert(g->num_dots <= max_dots);
2054
2055     grid_make_consistent(g);
2056     return g;
2057 }
2058
2059 /* ----------- End of grid generators ------------- */