/* * misc.c: Miscellaneous helpful functions. */ #include #include #ifdef NO_TGMATH_H # include #else # include #endif #include #include #include #include "puzzles.h" char MOVE_UI_UPDATE[] = ""; char MOVE_NO_EFFECT[] = ""; char MOVE_UNUSED[] = ""; void free_cfg(config_item *cfg) { config_item *i; for (i = cfg; i->type != C_END; i++) if (i->type == C_STRING) sfree(i->u.string.sval); sfree(cfg); } void free_keys(key_label *keys, int nkeys) { int i; for(i = 0; i < nkeys; i++) sfree(keys[i].label); sfree(keys); } /* * The Mines (among others) game descriptions contain the location of every * mine, and can therefore be used to cheat. * * It would be pointless to attempt to _prevent_ this form of * cheating by encrypting the description, since Mines is * open-source so anyone can find out the encryption key. However, * I think it is worth doing a bit of gentle obfuscation to prevent * _accidental_ spoilers: if you happened to note that the game ID * starts with an F, for example, you might be unable to put the * knowledge of those mines out of your mind while playing. So, * just as discussions of film endings are rot13ed to avoid * spoiling it for people who don't want to be told, we apply a * keyless, reversible, but visually completely obfuscatory masking * function to the mine bitmap. */ void obfuscate_bitmap(unsigned char *bmp, int bits, bool decode) { int bytes, firsthalf, secondhalf; struct step { unsigned char *seedstart; int seedlen; unsigned char *targetstart; int targetlen; } steps[2]; int i, j; /* * My obfuscation algorithm is similar in concept to the OAEP * encoding used in some forms of RSA. Here's a specification * of it: * * + We have a `masking function' which constructs a stream of * pseudorandom bytes from a seed of some number of input * bytes. * * + We pad out our input bit stream to a whole number of * bytes by adding up to 7 zero bits on the end. (In fact * the bitmap passed as input to this function will already * have had this done in practice.) * * + We divide the _byte_ stream exactly in half, rounding the * half-way position _down_. So an 81-bit input string, for * example, rounds up to 88 bits or 11 bytes, and then * dividing by two gives 5 bytes in the first half and 6 in * the second half. * * + We generate a mask from the second half of the bytes, and * XOR it over the first half. * * + We generate a mask from the (encoded) first half of the * bytes, and XOR it over the second half. Any null bits at * the end which were added as padding are cleared back to * zero even if this operation would have made them nonzero. * * To de-obfuscate, the steps are precisely the same except * that the final two are reversed. * * Finally, our masking function. Given an input seed string of * bytes, the output mask consists of concatenating the SHA-1 * hashes of the seed string and successive decimal integers, * starting from 0. */ bytes = (bits + 7) / 8; firsthalf = bytes / 2; secondhalf = bytes - firsthalf; steps[decode ? 1 : 0].seedstart = bmp + firsthalf; steps[decode ? 1 : 0].seedlen = secondhalf; steps[decode ? 1 : 0].targetstart = bmp; steps[decode ? 1 : 0].targetlen = firsthalf; steps[decode ? 0 : 1].seedstart = bmp; steps[decode ? 0 : 1].seedlen = firsthalf; steps[decode ? 0 : 1].targetstart = bmp + firsthalf; steps[decode ? 0 : 1].targetlen = secondhalf; for (i = 0; i < 2; i++) { SHA_State base, final; unsigned char digest[20]; char numberbuf[80]; int digestpos = 20, counter = 0; SHA_Init(&base); SHA_Bytes(&base, steps[i].seedstart, steps[i].seedlen); for (j = 0; j < steps[i].targetlen; j++) { if (digestpos >= 20) { sprintf(numberbuf, "%d", counter++); final = base; SHA_Bytes(&final, numberbuf, strlen(numberbuf)); SHA_Final(&final, digest); digestpos = 0; } steps[i].targetstart[j] ^= digest[digestpos++]; } /* * Mask off the pad bits in the final byte after both steps. */ if (bits % 8) bmp[bits / 8] &= 0xFF & (0xFF00 >> (bits % 8)); } } /* err, yeah, these two pretty much rely on unsigned char being 8 bits. * Platforms where this is not the case probably have bigger problems * than just making these two work, though... */ char *bin2hex(const unsigned char *in, int inlen) { char *ret = snewn(inlen*2 + 1, char), *p = ret; int i; for (i = 0; i < inlen*2; i++) { int v = in[i/2]; if (i % 2 == 0) v >>= 4; *p++ = "0123456789abcdef"[v & 0xF]; } *p = '\0'; return ret; } unsigned char *hex2bin(const char *in, int outlen) { unsigned char *ret = snewn(outlen, unsigned char); int i; memset(ret, 0, outlen*sizeof(unsigned char)); for (i = 0; i < outlen*2; i++) { int c = in[i]; int v; assert(c != 0); if (c >= '0' && c <= '9') v = c - '0'; else if (c >= 'a' && c <= 'f') v = c - 'a' + 10; else if (c >= 'A' && c <= 'F') v = c - 'A' + 10; else v = 0; ret[i / 2] |= v << (4 * (1 - (i % 2))); } return ret; } char *fgetline(FILE *fp) { char *ret = snewn(512, char); int size = 512, len = 0; while (fgets(ret + len, size - len, fp)) { len += strlen(ret + len); if (ret[len-1] == '\n') break; /* got a newline, we're done */ size = len + 512; ret = sresize(ret, size, char); } if (len == 0) { /* first fgets returned NULL */ sfree(ret); return NULL; } ret[len] = '\0'; return ret; } int getenv_bool(const char *name, int dflt) { char *env = getenv(name); if (env == NULL) return dflt; if (strchr("yYtT", env[0])) return true; return false; } /* Utility functions for colour manipulation. */ static float colour_distance(const float a[3], const float b[3]) { return (float)sqrt((a[0]-b[0]) * (a[0]-b[0]) + (a[1]-b[1]) * (a[1]-b[1]) + (a[2]-b[2]) * (a[2]-b[2])); } void colour_mix(const float src1[3], const float src2[3], float p, float dst[3]) { int i; for (i = 0; i < 3; i++) dst[i] = src1[i] * (1.0F - p) + src2[i] * p; } void game_mkhighlight_specific(frontend *fe, float *ret, int background, int highlight, int lowlight) { static const float black[3] = { 0.0F, 0.0F, 0.0F }; static const float white[3] = { 1.0F, 1.0F, 1.0F }; float db, dw; int i; /* * New geometric highlight-generation algorithm: Draw a line from * the base colour to white. The point K distance along this line * from the base colour is the highlight colour. Similarly, draw * a line from the base colour to black. The point on this line * at a distance K from the base colour is the shadow. If either * of these colours is imaginary (for reasonable K at most one * will be), _extrapolate_ the base colour along the same line * until it's a distance K from white (or black) and start again * with that as the base colour. * * This preserves the hue of the base colour, ensures that of the * three the base colour is the most saturated, and only ever * flattens the highlight and shadow to pure white or pure black. * * K must be at most sqrt(3)/2, or mid grey would be too close to * both white and black. Here K is set to sqrt(3)/6 so that this * code produces the same results as the former code in the common * case where the background is grey and the highlight saturates * to white. */ const float k = sqrt(3)/6.0F; if (lowlight >= 0) { db = colour_distance(&ret[background*3], black); if (db < k) { for (i = 0; i < 3; i++) ret[lowlight*3+i] = black[i]; if (db == 0.0F) colour_mix(black, white, k/sqrt(3), &ret[background*3]); else colour_mix(black, &ret[background*3], k/db, &ret[background*3]); } else { colour_mix(&ret[background*3], black, k/db, &ret[lowlight*3]); } } if (highlight >= 0) { dw = colour_distance(&ret[background*3], white); if (dw < k) { for (i = 0; i < 3; i++) ret[highlight*3+i] = white[i]; if (dw == 0.0F) colour_mix(white, black, k/sqrt(3), &ret[background*3]); else colour_mix(white, &ret[background*3], k/dw, &ret[background*3]); /* Background has changed; recalculate lowlight. */ if (lowlight >= 0) colour_mix(&ret[background*3], black, k/db, &ret[lowlight*3]); } else { colour_mix(&ret[background*3], white, k/dw, &ret[highlight*3]); } } } void game_mkhighlight(frontend *fe, float *ret, int background, int highlight, int lowlight) { frontend_default_colour(fe, &ret[background * 3]); game_mkhighlight_specific(fe, ret, background, highlight, lowlight); } static void memswap(void *av, void *bv, int size) { char tmpbuf[512]; char *a = av, *b = bv; while (size > 0) { int thislen = min(size, sizeof(tmpbuf)); memcpy(tmpbuf, a, thislen); memcpy(a, b, thislen); memcpy(b, tmpbuf, thislen); a += thislen; b += thislen; size -= thislen; } } void shuffle(void *array, int nelts, int eltsize, random_state *rs) { char *carray = (char *)array; int i; for (i = nelts; i-- > 1 ;) { int j = random_upto(rs, i+1); if (j != i) memswap(carray + eltsize * i, carray + eltsize * j, eltsize); } } void draw_rect_outline(drawing *dr, int x, int y, int w, int h, int colour) { int x0 = x, x1 = x+w-1, y0 = y, y1 = y+h-1; int coords[8]; coords[0] = x0; coords[1] = y0; coords[2] = x0; coords[3] = y1; coords[4] = x1; coords[5] = y1; coords[6] = x1; coords[7] = y0; draw_polygon(dr, coords, 4, -1, colour); } void draw_rect_corners(drawing *dr, int cx, int cy, int r, int col) { draw_line(dr, cx - r, cy - r, cx - r, cy - r/2, col); draw_line(dr, cx - r, cy - r, cx - r/2, cy - r, col); draw_line(dr, cx - r, cy + r, cx - r, cy + r/2, col); draw_line(dr, cx - r, cy + r, cx - r/2, cy + r, col); draw_line(dr, cx + r, cy - r, cx + r, cy - r/2, col); draw_line(dr, cx + r, cy - r, cx + r/2, cy - r, col); draw_line(dr, cx + r, cy + r, cx + r, cy + r/2, col); draw_line(dr, cx + r, cy + r, cx + r/2, cy + r, col); } void move_cursor(int button, int *x, int *y, int maxw, int maxh, bool wrap) { int dx = 0, dy = 0; switch (button) { case CURSOR_UP: dy = -1; break; case CURSOR_DOWN: dy = 1; break; case CURSOR_RIGHT: dx = 1; break; case CURSOR_LEFT: dx = -1; break; default: return; } if (wrap) { *x = (*x + dx + maxw) % maxw; *y = (*y + dy + maxh) % maxh; } else { *x = min(max(*x+dx, 0), maxw - 1); *y = min(max(*y+dy, 0), maxh - 1); } } /* Used in netslide.c and sixteen.c for cursor movement around edge. */ int c2pos(int w, int h, int cx, int cy) { if (cy == -1) return cx; /* top row, 0 .. w-1 (->) */ else if (cx == w) return w + cy; /* R col, w .. w+h -1 (v) */ else if (cy == h) return w + h + (w-cx-1); /* bottom row, w+h .. w+h+w-1 (<-) */ else if (cx == -1) return w + h + w + (h-cy-1); /* L col, w+h+w .. w+h+w+h-1 (^) */ assert(!"invalid cursor pos!"); return -1; /* not reached */ } int c2diff(int w, int h, int cx, int cy, int button) { int diff = 0; assert(IS_CURSOR_MOVE(button)); /* Obvious moves around edge. */ if (cy == -1) diff = (button == CURSOR_RIGHT) ? +1 : (button == CURSOR_LEFT) ? -1 : diff; if (cy == h) diff = (button == CURSOR_RIGHT) ? -1 : (button == CURSOR_LEFT) ? +1 : diff; if (cx == -1) diff = (button == CURSOR_UP) ? +1 : (button == CURSOR_DOWN) ? -1 : diff; if (cx == w) diff = (button == CURSOR_UP) ? -1 : (button == CURSOR_DOWN) ? +1 : diff; if (button == CURSOR_LEFT && cx == w && (cy == 0 || cy == h-1)) diff = (cy == 0) ? -1 : +1; if (button == CURSOR_RIGHT && cx == -1 && (cy == 0 || cy == h-1)) diff = (cy == 0) ? +1 : -1; if (button == CURSOR_DOWN && cy == -1 && (cx == 0 || cx == w-1)) diff = (cx == 0) ? -1 : +1; if (button == CURSOR_UP && cy == h && (cx == 0 || cx == w-1)) diff = (cx == 0) ? +1 : -1; debug(("cx,cy = %d,%d; w%d h%d, diff = %d", cx, cy, w, h, diff)); return diff; } void pos2c(int w, int h, int pos, int *cx, int *cy) { int max = w+h+w+h; pos = (pos + max) % max; if (pos < w) { *cx = pos; *cy = -1; return; } pos -= w; if (pos < h) { *cx = w; *cy = pos; return; } pos -= h; if (pos < w) { *cx = w-pos-1; *cy = h; return; } pos -= w; if (pos < h) { *cx = -1; *cy = h-pos-1; return; } assert(!"invalid pos, huh?"); /* limited by % above! */ } void draw_text_outline(drawing *dr, int x, int y, int fonttype, int fontsize, int align, int text_colour, int outline_colour, const char *text) { if (outline_colour > -1) { draw_text(dr, x-1, y, fonttype, fontsize, align, outline_colour, text); draw_text(dr, x+1, y, fonttype, fontsize, align, outline_colour, text); draw_text(dr, x, y-1, fonttype, fontsize, align, outline_colour, text); draw_text(dr, x, y+1, fonttype, fontsize, align, outline_colour, text); } draw_text(dr, x, y, fonttype, fontsize, align, text_colour, text); } /* kludge for sprintf() in Rockbox not supporting "%-8.8s" */ void copy_left_justified(char *buf, size_t sz, const char *str) { size_t len = strlen(str); assert(sz > 0); memset(buf, ' ', sz - 1); assert(len <= sz - 1); memcpy(buf, str, len); buf[sz - 1] = 0; } /* Returns a dynamically allocated label for a generic button. * Game-specific buttons should go into the `label' field of key_label * instead. */ char *button2label(int button) { /* check if it's a keyboard button */ if(('A' <= button && button <= 'Z') || ('a' <= button && button <= 'z') || ('0' <= button && button <= '9') ) { char str[2]; str[0] = button; str[1] = '\0'; return dupstr(str); } switch(button) { case CURSOR_UP: return dupstr("Up"); case CURSOR_DOWN: return dupstr("Down"); case CURSOR_LEFT: return dupstr("Left"); case CURSOR_RIGHT: return dupstr("Right"); case CURSOR_SELECT: return dupstr("Select"); case '\b': return dupstr("Clear"); default: fatal("unknown generic key"); } /* should never get here */ return NULL; } char *make_prefs_path(const char *dir, const char *sep, const game *game, const char *suffix) { size_t dirlen = strlen(dir); size_t seplen = strlen(sep); size_t gamelen = strlen(game->name); size_t suffixlen = strlen(suffix); char *path, *p; const char *q; if (!dir || !sep || !game || !suffix) return NULL; path = snewn(dirlen + seplen + gamelen + suffixlen + 1, char); p = path; memcpy(p, dir, dirlen); p += dirlen; memcpy(p, sep, seplen); p += seplen; for (q = game->name; *q; q++) if (*q != ' ') *p++ = tolower((unsigned char)*q); memcpy(p, suffix, suffixlen); p += suffixlen; *p = '\0'; return path; } /* vim: set shiftwidth=4 tabstop=8: */ 6 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 
/*
 * Code to generate patches of the aperiodic 'hat' tiling discovered
 * in 2023.
 *
 * This uses the 'combinatorial coordinates' system documented in my
 * public article
 * https://www.chiark.greenend.org.uk/~sgtatham/quasiblog/aperiodic-tilings/
 *
 * The internal document auxiliary/doc/hats.html also contains an
 * explanation of the basic ideas of this algorithm (less polished but
 * containing more detail).
 *
 * Neither of those documents can really be put in a source file,
 * because they just have too many complicated diagrams. So read at
 * least one of those first; the comments in here will refer to it.
 *
 * Discoverers' website: https://cs.uwaterloo.ca/~csk/hat/
 * Preprint of paper:    https://arxiv.org/abs/2303.10798
 */

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

#include "puzzles.h"
#include "hat.h"
#include "hat-internal.h"

void hat_kiteenum_first(KiteEnum *s, int w, int h)
{
    Kite start = { {0,0}, {0, 3}, {3, 0}, {2, 2} };
    size_t i;

    for (i = 0; i < KE_NKEEP; i++)
        s->recent[i] = start;          /* initialise to *something* */
    s->curr_index = 0;
    s->curr = &s->recent[s->curr_index];
    s->state = 1;
    s->w = w;
    s->h = h;
    s->x = 0;
    s->y = 0;
}

bool hat_kiteenum_next(KiteEnum *s)
{
    unsigned lastbut1 = s->last_index;
    s->last_index = s->curr_index;
    s->curr_index = (s->curr_index + 1) % KE_NKEEP;
    s->curr = &s->recent[s->curr_index];

    switch (s->state) {
        /* States 1,2,3 walk rightwards along the upper side of a
         * horizontal grid line with a pointy kite end at the start
         * point */
      case 1:
        s->last_step = KS_F_RIGHT;
        s->state = 2;
        break;

      case 2:
        if (s->x+1 >= s->w) {
            s->last_step = KS_F_RIGHT;
            s->state = 4;
            break;
        }
        s->last_step = KS_RIGHT;
        s->state = 3;
        s->x++;
        break;

      case 3:
        s->last_step = KS_RIGHT;
        s->state = 1;
        break;

        /* State 4 is special: we've just moved up into a row below a
         * grid line, but we can't produce the rightmost tile of that
         * row because it's not adjacent any tile so far emitted. So
         * instead, emit the second-rightmost tile, and next time,
         * we'll emit the rightmost. */
      case 4:
        s->last_step = KS_LEFT;
        s->state = 5;
        break;

        /* And now we have to emit the third-rightmost tile relative
         * to the last but one tile we emitted (the one from state 2,
         * not state 4). */
      case 5:
        s->last_step = KS_RIGHT;
        s->last_index = lastbut1;
        s->state = 6;
        break;

        /* Now states 6-8 handle the general case of walking leftwards
         * along the lower side of a line, starting from a
         * right-angled kite end. */
      case 6:
        if (s->x <= 0) {
            if (s->y+1 >= s->h) {
                s->state = 0;
                return false;
            }
            s->last_step = KS_RIGHT;
            s->state = 9;
            s->y++;
            break;
        }
        s->last_step = KS_F_RIGHT;
        s->state = 7;
        s->x--;
        break;

      case 7:
        s->last_step = KS_RIGHT;
        s->state = 8;
        break;

      case 8:
        s->last_step = KS_RIGHT;
        s->state = 6;
        break;

        /* States 9,10,11 walk rightwards along the upper side of a
         * horizontal grid line with a right-angled kite end at the
         * start point. This time there's no awkward transition from
         * the previous row. */
      case 9:
        s->last_step = KS_RIGHT;
        s->state = 10;
        break;

      case 10:
        s->last_step = KS_RIGHT;
        s->state = 11;
        break;

      case 11:
        if (s->x+1 >= s->w) {
            /* Another awkward transition to the next row, where we
             * have to generate it based on the previous state-9 tile.
             * But this time at least we generate the rightmost tile
             * of the new row, so the next states will be simple. */
            s->last_step = KS_F_RIGHT;
            s->last_index = lastbut1;
            s->state = 12;
            break;
        }
        s->last_step = KS_F_RIGHT;
        s->state = 9;
        s->x++;
        break;

        /* States 12,13,14 walk leftwards along the upper edge of a
         * horizontal grid line with a pointy kite end at the start
         * point */
      case 12:
        s->last_step = KS_F_RIGHT;
        s->state = 13;
        break;

      case 13:
        if (s->x <= 0) {
            if (s->y+1 >= s->h) {
                s->state = 0;
                return false;
            }
            s->last_step = KS_LEFT;
            s->state = 1;
            s->y++;
            break;
        }
        s->last_step = KS_RIGHT;
        s->state = 14;
        s->x--;
        break;

      case 14:
        s->last_step = KS_RIGHT;
        s->state = 12;
        break;

      default:
        return false;
    }

    *s->curr = kite_step(s->recent[s->last_index], s->last_step);
    return true;
}

/*
 * The actual tables.
 */
#include "hat-tables.h"

/*
 * One set of tables that we write by hand: the permitted ways to
 * extend the coordinate system outwards from a given metatile.
 *
 * One obvious approach would be to make a table of all the places
 * each metatile can appear in the expansion of another (e.g. H can be
 * subtile 0, 1 or 2 of another H, subtile 0 of a T, or 0 or 1 of a P
 * or an F), and when we need to decide what our current topmost tile
 * turns out to be a subtile of, choose equiprobably at random from
 * those options.
 *
 * That's what I did originally, but a better approach is to skew the
 * probabilities. We'd like to generate our patch of actual tiling
 * uniformly at random, in the sense that if you selected uniformly
 * from a very large region of the plane, the distribution of possible
 * finite patches of tiling would converge to some limit as that
 * region tended to infinity, and we'd be picking from that limiting
 * distribution on finite patches.
 *
 * For this we have to refer back to the original paper, which
 * indicates the subset of each metatile's expansion that can be
 * considered to 'belong' to that metatile, such that every subtile
 * belongs to exactly one parent metatile, and the overlaps are
 * eliminated. Reading out the diagrams from their Figure 2.8:
 *
 * - H: we discard three of the outer F subtiles, in the symmetric
 *    positions index by our coordinates as 7, 10, 11. So we keep the
 *    remaining subtiles {0,1,2,3,4,5,6,8,9,12}, which consist of
 *    three H, one T, three P and three F.
 *
 * - T: only the central H expanded from a T is considered to belong
 *   to it, so we just keep {0}, a single H.
 *
 * - P: we discard everything intersected by a long edge of the
 *   parallelogram, leaving the central three tiles and the endmost
 *   pair of F. That is, we keep {0,1,4,5,10}, consisting of two H,
 *   one P and two F.
 *
 * - F: looks like P at one end, and we retain the corresponding set
 *   of tiles there, but at the other end we keep the two F on either
 *   side of the endmost one. So we keep {0,1,3,6,8,10}, consisting of
 *   two H, one P and _three_ F.
 *
 * Adding up the tile numbers gives us this matrix system:
 *
 * (H_1)   (3 1 2 2)(H_0)
 * (T_1) = (1 0 0 0)(T_0)
 * (P_1)   (3 0 1 1)(P_0)
 * (F_1)   (3 0 2 3)(F_0)
 *
 * which says that if you have a patch of metatiling consisting of H_0
 * H tiles, T_0 T tiles etc, then this matrix shows the number H_1 of
 * smaller H tiles, etc, expanded from it.
 *
 * If you expand _many_ times, that's equivalent to raising the matrix
 * to a power:
 *
 *                  n
 * (H_n)   (3 1 2 2) (H_0)
 * (T_n) = (1 0 0 0) (T_0)
 * (P_n)   (3 0 1 1) (P_0)
 * (F_n)   (3 0 2 3) (F_0)
 *
 * The limiting distribution of metatiles is obtained by looking at
 * the four-way ratio between H_n, T_n, P_n and F_n as n tends to
 * infinity. To calculate this, we find the eigenvalues and
 * eigenvectors of the matrix, and extract the eigenvector
 * corresponding to the eigenvalue of largest magnitude. (Things get
 * more complicated in cases where there isn't a _unique_ eigenvalue
 * of largest magnitude, but here, there is.)
 *
 * That eigenvector is
 *
 *   [          1          ]      [ 1                      ]
 *   [ (7 - 3 sqrt(5)) / 2 ]  ~=  [ 0.14589803375031545538 ]
 *   [    3 sqrt(5) - 6    ]      [ 0.70820393249936908922 ]
 *   [ (9 - 3 sqrt(5)) / 2 ]      [ 1.14589803375031545538 ]
 *
 * So those are the limiting relative proportions of metatiles.
 *
 * So if we have a particular metatile, how likely is it for its
 * parent to be one of those? We have to adjust by the number of
 * metatiles of each type that each tile has as its children. For
 * example, the P and F tiles have one P child each, but the H has
 * three P children. So if we have a P, the proportion of H in its
 * potential ancestry is three times what's shown here. (And T can't
 * occur at all as a parent.)
 *
 * In other words, we should choose _each coordinate_ with probability
 * corresponding to one of those numbers (scaled down so they all sum
 * to 1). Continuing to use P as an example, it will be:
 *
 *  - child 4 of H with relative probability 1
 *  - child 5 of H with relative probability 1
 *  - child 6 of H with relative probability 1
 *  - child 4 of P with relative probability 0.70820393249936908922
 *  - child 3 of F with relative probability 1.14589803375031545538
 *
 * and then we obtain the true probabilities by scaling those values
 * down so that they sum to 1.
 *
 * The tables below give a reasonable approximation in 32-bit
 * integers to these proportions.
 */

typedef struct MetatilePossibleParent {
    TileType type;
    unsigned index;
    unsigned long probability;
} MetatilePossibleParent;

/* The above probabilities scaled up by 10000000 */
#define PROB_H 10000000
#define PROB_T  1458980
#define PROB_P  7082039
#define PROB_F 11458980

static const MetatilePossibleParent parents_H[] = {
    { TT_H,  0, PROB_H },
    { TT_H,  1, PROB_H },
    { TT_H,  2, PROB_H },
    { TT_T,  0, PROB_T },
    { TT_P,  0, PROB_P },
    { TT_P,  1, PROB_P },
    { TT_F,  0, PROB_F },
    { TT_F,  1, PROB_F },
};
static const MetatilePossibleParent parents_T[] = {
    { TT_H,  3, PROB_H },
};
static const MetatilePossibleParent parents_P[] = {
    { TT_H,  4, PROB_H },
    { TT_H,  5, PROB_H },
    { TT_H,  6, PROB_H },
    { TT_P,  4, PROB_P },
    { TT_F,  3, PROB_F },
};
static const MetatilePossibleParent parents_F[] = {
    { TT_H,  8, PROB_H },
    { TT_H,  9, PROB_H },
    { TT_H, 12, PROB_H },
    { TT_P,  5, PROB_P },
    { TT_P, 10, PROB_P },
    { TT_F,  6, PROB_F },
    { TT_F,  8, PROB_F },
    { TT_F, 10, PROB_F },
};

static const MetatilePossibleParent *const possible_parents[] = {
    parents_H, parents_T, parents_P, parents_F,
};
static const size_t n_possible_parents[] = {
    lenof(parents_H), lenof(parents_T), lenof(parents_P), lenof(parents_F),
};

/*
 * Similarly, we also want to choose our absolute starting hat with
 * close to uniform probability, which again we do by looking at the
 * limiting ratio of the metatile types, and this time, scaling by the
 * number of hats in each metatile.
 *
 * We cheatingly use the same MetatilePossibleParent struct, because
 * it's got all the right fields, even if it has an inappropriate
 * name.
 */
static const MetatilePossibleParent starting_hats[] = {
    { TT_H,  0, PROB_H },
    { TT_H,  1, PROB_H },
    { TT_H,  2, PROB_H },
    { TT_H,  3, PROB_H },
    { TT_T,  0, PROB_P },
    { TT_P,  0, PROB_P },
    { TT_P,  1, PROB_P },
    { TT_F,  0, PROB_F },
    { TT_F,  1, PROB_F },
};

#undef PROB_H
#undef PROB_T
#undef PROB_P
#undef PROB_F

HatCoords *hat_coords_new(void)
{
    HatCoords *hc = snew(HatCoords);
    hc->nc = hc->csize = 0;
    hc->c = NULL;
    return hc;
}

void hat_coords_free(HatCoords *hc)
{
    if (hc) {
        sfree(hc->c);
        sfree(hc);
    }
}

void hat_coords_make_space(HatCoords *hc, size_t size)
{
    if (hc->csize < size) {
        hc->csize = hc->csize * 5 / 4 + 16;
        if (hc->csize < size)
            hc->csize = size;
        hc->c = sresize(hc->c, hc->csize, HatCoord);
    }
}

HatCoords *hat_coords_copy(HatCoords *hc_in)
{
    HatCoords *hc_out = hat_coords_new();
    hat_coords_make_space(hc_out, hc_in->nc);
    memcpy(hc_out->c, hc_in->c, hc_in->nc * sizeof(*hc_out->c));
    hc_out->nc = hc_in->nc;
    return hc_out;
}

static const MetatilePossibleParent *choose_mpp(
    random_state *rs, const MetatilePossibleParent *parents, size_t nparents)
{
    /*
     * If we needed to do this _efficiently_, we'd rewrite all those
     * tables above as cumulative frequency tables and use binary
     * search. But this happens about log n times in a grid of area n,
     * so it hardly matters, and it's easier to keep the tables
     * legible.
     */
    unsigned long limit = 0, value;
    size_t i;

    for (i = 0; i < nparents; i++)
        limit += parents[i].probability;

    value = random_upto(rs, limit);

    for (i = 0; i+1 < nparents; i++) {
        if (value < parents[i].probability)
            return &parents[i];
        value -= parents[i].probability;
    }

    assert(i == nparents - 1);
    assert(value < parents[i].probability);
    return &parents[i];
}
void hatctx_init_random(HatContext *ctx, random_state *rs)
{
    const MetatilePossibleParent *starting_hat = choose_mpp(
        rs, starting_hats, lenof(starting_hats));

    ctx->rs = rs;
    ctx->prototype = hat_coords_new();
    hat_coords_make_space(ctx->prototype, 3);
    ctx->prototype->c[2].type = starting_hat->type;
    ctx->prototype->c[2].index = -1;
    ctx->prototype->c[1].type = TT_HAT;
    ctx->prototype->c[1].index = starting_hat->index;
    ctx->prototype->c[0].type = TT_KITE;
    ctx->prototype->c[0].index = random_upto(rs, HAT_KITES);
    ctx->prototype->nc = 3;
}

static inline int metatile_char_to_enum(char metatile)
{
    return (metatile == 'H' ? TT_H :
            metatile == 'T' ? TT_T :
            metatile == 'P' ? TT_P :
            metatile == 'F' ? TT_F : -1);
}

static void init_coords_params(HatContext *ctx,
                               const struct HatPatchParams *hp)
{
    size_t i;

    ctx->rs = NULL;
    ctx->prototype = hat_coords_new();

    assert(hp->ncoords >= 3);

    hat_coords_make_space(ctx->prototype, hp->ncoords + 1);
    ctx->prototype->nc = hp->ncoords + 1;

    for (i = 0; i < hp->ncoords; i++)
        ctx->prototype->c[i].index = hp->coords[i];

    ctx->prototype->c[hp->ncoords].type =
        metatile_char_to_enum(hp->final_metatile);
    ctx->prototype->c[hp->ncoords].index = -1;

    ctx->prototype->c[0].type = TT_KITE;
    ctx->prototype->c[1].type = TT_HAT;

    for (i = hp->ncoords - 1; i > 1; i--) {
        TileType metatile = ctx->prototype->c[i+1].type;
        assert(hp->coords[i] < nchildren[metatile]);
        ctx->prototype->c[i].type = children[metatile][hp->coords[i]];
    }

    assert(hp->coords[0] < 8);
}

HatCoords *hatctx_initial_coords(HatContext *ctx)
{
    return hat_coords_copy(ctx->prototype);
}

/*
 * Extend hc until it has at least n coordinates in, by copying from
 * ctx->prototype if needed, and extending ctx->prototype if needed in
 * order to do that.
 */
void hatctx_extend_coords(HatContext *ctx, HatCoords *hc, size_t n)
{
    if (ctx->prototype->nc < n) {
        hat_coords_make_space(ctx->prototype, n);
        while (ctx->prototype->nc < n) {
            TileType type = ctx->prototype->c[ctx->prototype->nc - 1].type;
            assert(ctx->prototype->c[ctx->prototype->nc - 1].index == -1);
            const MetatilePossibleParent *parent;

            if (ctx->rs)
                parent = choose_mpp(ctx->rs, possible_parents[type],
                                    n_possible_parents[type]);
            else
                parent = possible_parents[type];

            ctx->prototype->c[ctx->prototype->nc - 1].index = parent->index;
            ctx->prototype->c[ctx->prototype->nc].index = -1;
            ctx->prototype->c[ctx->prototype->nc].type = parent->type;
            ctx->prototype->nc++;
        }
    }

    hat_coords_make_space(hc, n);
    while (hc->nc < n) {
        assert(hc->c[hc->nc - 1].index == -1);
        assert(hc->c[hc->nc - 1].type == ctx->prototype->c[hc->nc - 1].type);
        hc->c[hc->nc - 1].index = ctx->prototype->c[hc->nc - 1].index;
        hc->c[hc->nc].index = -1;
        hc->c[hc->nc].type = ctx->prototype->c[hc->nc].type;
        hc->nc++;
    }
}

void hatctx_cleanup(HatContext *ctx)
{
    hat_coords_free(ctx->prototype);
}

/*
 * The actual system for finding the coordinates of an adjacent kite.
 */

/*
 * Kitemap step: ensure we have enough coordinates to know two levels
 * of meta-tiling, and use the kite map for the outer layer to move
 * around the individual kites. If this fails, return NULL.
 */
static HatCoords *try_step_coords_kitemap(
    HatContext *ctx, HatCoords *hc_in, KiteStep step)
{
    hatctx_extend_coords(ctx, hc_in, 4);
    hat_coords_debug("  try kitemap  ", hc_in, "\n");
    unsigned kite = hc_in->c[0].index;
    unsigned hat = hc_in->c[1].index;
    unsigned meta = hc_in->c[2].index;
    TileType meta2type = hc_in->c[3].type;
    const KitemapEntry *ke = &kitemap[meta2type][
        kitemap_index(step, kite, hat, meta)];
    if (ke->kite >= 0) {
        /*
         * Success! We've got coordinates for the next kite in this
         * direction.
         */
        HatCoords *hc_out = hat_coords_copy(hc_in);

        hc_out->c[2].index = ke->meta;
        hc_out->c[2].type = children[meta2type][ke->meta];
        hc_out->c[1].index = ke->hat;
        hc_out->c[1].type = TT_HAT;
        hc_out->c[0].index = ke->kite;
        hc_out->c[0].type = TT_KITE;

        hat_coords_debug("  success!     ", hc_out, "\n");
        return hc_out;
    }

    return NULL;
}

/*
 * Recursive metamap step. Try using the metamap to rewrite the
 * coordinates at hc->c[depth] and hc->c[depth+1] (using the metamap
 * for the tile type described in hc->c[depth+2]). If successful,
 * recurse back down to see if this led to a successful step via the
 * kitemap. If even that fails (so that we need to try a higher-order
 * metamap rewrite), return NULL.
 */
static HatCoords *try_step_coords_metamap(
    HatContext *ctx, HatCoords *hc_in, KiteStep step, size_t depth)
{
    HatCoords *hc_tmp = NULL, *hc_out;

    hatctx_extend_coords(ctx, hc_in, depth+3);
#ifdef HAT_COORDS_DEBUG
    fprintf(stderr, "  try meta %-4d", (int)depth);
    hat_coords_debug("", hc_in, "\n");
#endif
    unsigned meta_orig = hc_in->c[depth].index;
    unsigned meta2_orig = hc_in->c[depth+1].index;
    TileType meta3type = hc_in->c[depth+2].type;

    unsigned meta = meta_orig, meta2 = meta2_orig;

    while (true) {
        const MetamapEntry *me;
        HatCoords *hc_curr = hc_tmp ? hc_tmp : hc_in;

        if (depth > 2)
            hc_out = try_step_coords_metamap(ctx, hc_curr, step, depth - 1);
        else
            hc_out = try_step_coords_kitemap(ctx, hc_curr, step);
        if (hc_out) {
            hat_coords_free(hc_tmp);
            return hc_out;
        }

        me = &metamap[meta3type][metamap_index(meta, meta2)];
        assert(me->meta != -1);
        if (me->meta == meta_orig && me->meta2 == meta2_orig) {
            hat_coords_free(hc_tmp);
            return NULL;
        }

        meta = me->meta;
        meta2 = me->meta2;

        /*
         * We must do the rewrite in a copy of hc_in. It's not
         * _necessarily_ obvious that that's the case (any successful
         * rewrite leaves the coordinates still valid and still
         * referring to the same kite, right?). But the problem is
         * that we might do a rewrite at this level more than once,