Fix false positive in bfzIII word detection

[clav.git] / mutation-finder.c

/*
 * Copyright (c) 2016-2019, S. Gilles <sgilles@math.umd.edu>
 *
 * Permission to use, copy, modify, and/or distribute this software
 * for any purpose with or without fee is hereby granted, provided
 * that the above copyright notice and this permission notice appear
 * in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
 * WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE
 * AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
 * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
 * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <errno.h>
#include <limits.h>
#include <locale.h>
#include <signal.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include "macros.h"
#include "quiver.h"

/* Flag for whether to print info */
static volatile sig_atomic_t status_requested = 0;

/* Signal handler */
static void handle(int sig);

/* Signal handler */
static void
handle(int sig)
{
        /* XXX: I'd like to check for SIGINFO here */
        if (sig == SIGUSR1) {
                status_requested = 1;
                signal(sig, handle);
        }
}

/* Print usage */
static void
usage(char *argv0)
{
        printf("Usage: %s -s <startfile> -e <endfile>\n", argv0);
        printf("       [ -v ] # Be verbose\n");
        printf("       [ -R ] # Require all non-frozen in sequence\n");
        printf("       [ -H ] # Test via graph heuristics, not exactness\n");
        printf("       [ -c <edgecap> ]\n");
        printf("       [ -f <frozen> [ -f <frozen> [...] ] ]\n");
        printf("       [ -l <startlen> ] [ -m <maxlen> ]\n");
        printf("\n");
        printf("       startfile: a quiver file representing the start "
               "state\n");
        printf("       endfile:   a quiver file representing the desired end "
               "state\n");
        printf("       edgecap:   ignore mutations if the cause edges with"
               "weights above this\n");
        printf("       startlen:  start with sequences of this length\n");
        printf("       maxlen:    halt after exhausting sequences of this "
               "length\n");
        printf("       frozen:    a vertex with this name will not be "
               "mutated\n");
}

static int
cmp(const void *a, const void *b)
{
        struct vertex **va = (struct vertex **) a;
        struct vertex **vb = (struct vertex **) b;

        return strcmp((*va)->name, (*vb)->name);
}

/* This one just counts the number of edges */
static uint_fast8_t
dumb_checksum1(struct quiver const * const q)
{
        uint_fast8_t c = 0;
        struct rational *r = 0;
        struct rational *s = 0;

        for (size_t j = 0; j < q->v_num; ++j) {
                for (size_t k = 0; k < j; ++k) {
                        r = &q->e[j * q->v_len + k];
                        s = &q->e[k * q->v_len + j];

                        if (!!(r->p ||
                               s->p)) {
                                c++;
                        }
                }
        }

        return c;
}

/* This one counts the number of each p and q */
static uint64_t
dumb_checksum2(struct quiver const * const q)
{
        uint64_t c = 0;
        struct rational *r = 0;
        struct rational *s = 0;
        size_t by_p[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
        size_t by_q[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };

        for (size_t j = 0; j < q->v_num; ++j) {
                for (size_t k = 0; k < q->v_num; ++k) {
                        r = &q->e[j * q->v_len + k];
                        s = &q->e[k * q->v_len + j];
                        by_p[(((4 + r->p) % 8) + 8) % 8]++;
                        by_q[r->q % 8]++;
                        by_p[(((4 + s->p) % 8) + 8) % 8]++;
                        by_q[s->q % 8]++;
                }
        }

        c = 8 * 0 + by_p[0] + 4 * by_q[0];
        c = 8 * c + by_p[1] + 4 * by_q[1];
        c = 8 * c + by_p[2] + 4 * by_q[2];
        c = 8 * c + by_p[3] + 4 * by_q[3];
        c = 8 * c + by_p[4] + 4 * by_q[4];
        c = 8 * c + by_p[5] + 4 * by_q[5];
        c = 8 * c + by_p[6] + 4 * by_q[6];
        c = 8 * c + by_p[7] + 4 * by_q[7];

        return c;
}

/* This one counts the number of vertices of each valence */
static uint64_t
dumb_checksum3(struct quiver const * const q)
{
        uint64_t c = 0;
        struct rational *r = 0;
        struct rational *s = 0;
        size_t by_v[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };

        for (size_t j = 0; j < q->v_num; ++j) {
                size_t v = 0;

                for (size_t k = 0; k < q->v_num; ++k) {
                        r = &q->e[j * q->v_len + k];
                        s = &q->e[k * q->v_len + j];

                        if (r->p ||
                            s->p) {
                                v++;
                        }
                }

                by_v[v % 8] += (1 << q->v[j].fatness);
        }

        c = by_v[0];
        c = 8 * c + by_v[1];
        c = 8 * c + by_v[2];
        c = 8 * c + by_v[3];
        c = 8 * c + by_v[4];
        c = 8 * c + by_v[5];
        c = 8 * c + by_v[6];
        c = 8 * c + by_v[7];

        return c;
}

/* This one counts the valence types (in/out/frozen/non-frozen) mod 4 */
static uint64_t
dumb_checksum4(struct quiver const * const q, size_t num_nonfrozen)
{
        uint64_t c = 0;
        struct rational *r = 0;
        size_t by_i[4] = { 0, 0, 0, 0 };
        size_t by_o[4] = { 0, 0, 0, 0 };
        size_t by_n[4] = { 0, 0, 0, 0 };
        size_t by_f[4] = { 0, 0, 0, 0 };

        for (size_t j = 0; j < q->v_num; ++j) {
                size_t v_i = 0;
                size_t v_o = 0;
                size_t v_n = 0;
                size_t v_f = 0;

                for (size_t k = 0; k < q->v_num; ++k) {
                        r = &q->e[j * q->v_len + k];

                        if (r->p > 0) {
                                v_o++;
                        } else if (r->p < 0) {
                                v_i++;
                        } else {
                                continue;
                        }

                        if (j >= num_nonfrozen &&
                            k >= num_nonfrozen) {
                                v_f++;
                        } else if (j < num_nonfrozen &&
                                   k < num_nonfrozen) {
                                v_n++;
                        }
                }

                by_i[v_i % 4] += (1 << q->v[j].fatness);
                by_o[v_o % 4] += (1 << q->v[j].fatness);
                by_n[v_n % 4] += (1 << q->v[j].fatness);
                by_f[v_f % 4] += (1 << q->v[j].fatness);
        }

        c = 0;

        for (size_t j = 0; j < 4; ++j) {
                c = (c << 8);
                c = c + by_i[j] + 4 * (by_o[j] + 4 * (by_n[j] + 4 * by_f[j]));
        }

        return c;
}

/* This one counts number of floor(p*6/q) */
static uint64_t
dumb_checksum5(struct quiver const * const q)
{
        uint64_t c = 0;
        struct rational *r = 0;
        struct rational *s = 0;
        size_t by_e[16] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };

        for (size_t j = 0; j < q->v_num; ++j) {
                for (size_t k = 0; k < q->v_num; ++k) {
                        r = &q->e[j * q->v_len + k];
                        s = &q->e[k * q->v_len + j];
                        by_e[(((((r->p * 6) / r->q) + 3) % 16) + 16) % 16]++;
                        by_e[(((((s->p * 6) / s->q) + 3) % 16) + 16) % 16]++;
                }
        }

        c = 8 * 0 + by_e[0];
        c = 8 * c + by_e[1];
        c = 8 * c + by_e[2];
        c = 8 * c + by_e[3];
        c = 8 * c + by_e[4];
        c = 8 * c + by_e[5];
        c = 8 * c + by_e[6];
        c = 8 * c + by_e[7];

        return c;
}

/* This one computes rk(H_1(G)), so actually takes some effort */
static uint_fast8_t
dumb_checksum6(struct quiver const *q, int oriented)
{
        uint_fast8_t r = 0;
        uint_fast8_t *touched = 0;
        size_t *frontier1 = 0;
        size_t *frontier2 = 0;
        size_t *temp = 0;

        /*
         * This DFS implementation is absolute garbage, for the
         * sake of simplicity over speed. It's not in a hot path,
         * anyway.
         */
        if (!(touched = calloc(q->v_num, sizeof *touched))) {
                goto done;
        }

        if (!(frontier1 = calloc(q->v_num + 1, sizeof *frontier1))) {
                goto done;
        }

        if (!(frontier2 = calloc(q->v_num + 1, sizeof *frontier2))) {
                goto done;
        }

        for (size_t k = 0; k <= q->v_num; ++k) {
                frontier1[k] = (size_t) -1;
                frontier2[k] = (size_t) -1;
        }

        for (size_t j = 0; j < q->v_num; ++j) {
                int need_another_pass = 1;

                if (touched[j]) {
                        continue;
                }

                touched[j] = 1;
                frontier1[0] = j;
                frontier1[1] = (size_t) -1;

                while (need_another_pass) {
                        size_t fidx = 0;

                        need_another_pass = 0;

                        for (size_t k = 0; k < q->v_num; ++k) {
                                size_t kidx = frontier1[k];

                                if (kidx == (size_t) -1) {
                                        break;
                                }

                                touched[kidx] = 1;

                                for (size_t l = 0; l < q->v_num; ++l) {
                                        if (l == kidx) {
                                                continue;
                                        }

                                        struct rational *e = &q->e[kidx *
                                                                   q->v_len +
                                                                   l];
                                        struct rational *f = &q->e[l *
                                                                   q->v_len +
                                                                   kidx];

                                        if (e->p ||
                                            (oriented &&
                                             f->p)) {
                                                if (touched[l]) {
                                                        r++;
                                                } else {
                                                        frontier2[fidx] = l;
                                                        fidx++;
                                                        frontier2[fidx] =
                                                                (size_t) -1;
                                                        need_another_pass = 1;
                                                        touched[l] = 1;
                                                }
                                        }
                                }
                        }

                        temp = frontier1;
                        frontier1 = frontier2;
                        frontier2 = temp;
                }
        }

        /*
         * Note that r is not actually the rank in the non-oriented
         * case. We've double-counted every edge exactly once. No
         * matter, we only care about using this as an invariant.
         */
done:
        free(touched);
        free(frontier1);
        free(frontier2);

        return r;
}

/* Load start/end quivers, figure out frozen vertices, order things, etc. */
static int
setup(char *start_file, char *end_file, char **frozen, size_t num_frozen, struct
      quiver *out_qstart, struct quiver *out_qend, size_t *out_num_nonfrozen)
{
        struct quiver a = { 0 };
        struct quiver b = { 0 };
        FILE *fa = 0;
        FILE *fb = 0;
        const char *errstr = 0;
        int ret = 0;
        struct vertex **averts = 0;
        struct vertex **bverts = 0;
        size_t *a_to_start = 0;
        size_t *b_to_end = 0;

        if (!out_qstart ||
            !out_qend) {
                fprintf(stderr, L("invalid quivers"));
                ret = EINVAL;
                goto done;
        }

        if (out_qstart->v_num ||
            out_qend->v_num) {
                fprintf(stderr, L("output quivers are not empty\n"));
                ret = EINVAL;
                goto done;
        }

        if (!(fa = fopen(start_file, "r"))) {
                ret = errno;
                perror(L("fopen"));
                fprintf(stderr, "cannot open file %s\n", start_file);
                goto done;
        }

        if (!(fb = fopen(end_file, "r"))) {
                ret = errno;
                perror(L("fopen"));
                fprintf(stderr, "cannot open file %s\n", end_file);
                goto done;
        }

        if ((ret = quiver_load_from_file(&a, fa, &errstr))) {
                fprintf(stderr, "%s\n", errstr);
                goto done;
        }

        if ((ret = quiver_load_from_file(&b, fb, &errstr))) {
                fprintf(stderr, "%s\n", errstr);
                goto done;
        }

        if (a.v_num != b.v_num) {
                fprintf(stderr, L(
                                "Quivers do not have same number of vertices\n"));
                ret = EINVAL;
                goto done;
        }

        /*
         * Make sure there are no duplicate vertices, but don't sort in
         * place - that would change ordering, and ordering of vertices
         * must match user input - tuning input so that known prefixes
         * arise early is sometimes crucial.
         */
        if (!(averts = calloc(a.v_num, sizeof *averts))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        if (!(bverts = calloc(b.v_num, sizeof *bverts))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        for (size_t j = 0; j < a.v_num; ++j) {
                averts[j] = &(a.v[j]);
                bverts[j] = &(b.v[j]);
        }

        qsort(averts, a.v_num, sizeof(averts[0]), cmp);
        qsort(bverts, b.v_num, sizeof(bverts[0]), cmp);

        for (size_t j = 0; j < a.v_num; ++j) {
                if (j &&
                    !strcmp(averts[j]->name, averts[j - 1]->name)) {
                        fprintf(stderr, L("Duplicate vertex \"%s\"\n"),
                                averts[j]->name);
                        ret = EINVAL;
                        goto done;
                }

                int d = strcmp(averts[j]->name, bverts[j]->name);

                if (d > 0) {
                        fprintf(stderr, L("Vertex \"%s\" not in both files\n"),
                                bverts[j]->name);
                        ret = EINVAL;
                        goto done;
                }

                if (d < 0) {
                        fprintf(stderr, L("Vertex \"%s\" not in both files\n"),
                                averts[j]->name);
                        ret = EINVAL;
                        goto done;
                }

                if (averts[j]->fatness != bverts[j]->fatness) {
                        fprintf(stderr, L(
                                        "Vertex \"%s\" has inconsistent fatness\n"),
                                averts[j]->name);
                        ret = EINVAL;
                        goto done;
                }
        }

        /* a_to_s[i] = j <=> "vertex i in a is vertex j in out_qstart" */
        if (!(a_to_start = calloc(a.v_num, sizeof(*a_to_start)))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        if (!(b_to_end = calloc(b.v_num, sizeof(*b_to_end)))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        /*
         * Order so that frozen vertices are all at the end - required
         * for incrementation algorithm
         */
        for (uint_fast8_t want_frozen = 0; want_frozen <= 1; want_frozen++) {
                for (size_t j = 0; j < a.v_num; ++j) {
                        uint_fast8_t is_frozen = 0;

                        for (size_t k = 0; k < num_frozen; ++k) {
                                is_frozen = is_frozen ||
                                            !strcmp(a.v[j].name, frozen[k]);

                                if (is_frozen) {
                                        k = num_frozen;
                                }
                        }

                        if (is_frozen != want_frozen) {
                                continue;
                        }

                        if ((ret = quiver_add_vertex(out_qstart,
                                                     &(a_to_start[j]),
                                                     a.v[j].name,
                                                     a.v[j].fatness, 0, 0, 0,
                                                     &errstr))) {
                                fprintf(stderr, "%s\n", errstr);
                                goto done;
                        }

                        if (!want_frozen &&
                            out_num_nonfrozen) {
                                *out_num_nonfrozen = a_to_start[j] + 1;
                        }

                        size_t tj = 0;

                        if ((ret = quiver_add_vertex(out_qend, &tj, a.v[j].name,
                                                     a.v[j].fatness, 0, 0, 0,
                                                     &errstr))) {
                                fprintf(stderr, "%s\n", errstr);
                                goto done;
                        }

                        for (size_t k = 0; k < b.v_num; ++k) {
                                if (!(strcmp(b.v[k].name,
                                             out_qend->v[tj].name))) {
                                        b_to_end[k] = tj;
                                }
                        }
                }
        }

        /* Add the edges in correct order */
        for (size_t j = 0; j < a.v_num; ++j) {
                size_t tj = a_to_start[j];

                for (size_t k = 0; k < a.v_num; ++k) {
                        size_t tk = a_to_start[k];
                        size_t idx_o = tj * out_qstart->v_len + tk;
                        size_t idx_a = j * a.v_len + k;

                        out_qstart->e[idx_o] = a.e[idx_a];
                }
        }

        for (size_t j = 0; j < b.v_num; ++j) {
                size_t tj = b_to_end[j];

                for (size_t k = 0; k < b.v_num; ++k) {
                        size_t tk = b_to_end[k];
                        size_t idx_o = tj * out_qend->v_len + tk;
                        size_t idx_b = j * b.v_len + k;

                        out_qend->e[idx_o] = b.e[idx_b];
                }
        }

done:
        quiver_clean(&a);
        quiver_clean(&b);

        if (fa) {
                fclose(fa);
        }

        if (fb) {
                fclose(fb);
        }

        free(averts);
        free(bverts);
        free(a_to_start);
        free(b_to_end);

        return ret;
}

/* Main loop */
int
main(int argc, char **argv)
{
        int ret = 0;
        struct quiver q = { 0 };
        struct quiver qend = { 0 };
        char *frozen[argc];
        size_t fidx = 0;
        char *start_file = 0;
        char *end_file = 0;
        long start_len_l = 0;
        long max_len_l = 0;
        size_t start_len = 0;
        size_t max_len = (size_t) -1;
        size_t warn_level = (size_t) -1;
        struct rational *target = 0;
        char *p = 0;
        const char *errstr = 0;
        uint_fast8_t verbose = 0;
        uint_fast8_t require_all_nonfrozen = 0;
        uint_fast8_t check_heuristics_only = 0;
        int *seen_in_seq = 0;
        uint_fast8_t desired_dumb_checksum1 = 0;
        uint64_t desired_dumb_checksum2 = 0;
        uint64_t desired_dumb_checksum3 = 0;
        uint64_t desired_dumb_checksum4 = 0;
        uint64_t desired_dumb_checksum5 = 0;
        uint8_t desired_dumb_checksum6u = 0;
        uint8_t desired_dumb_checksum6o = 0;
        uint_fast8_t this_dumb_checksum1 = 0;
        uint64_t this_dumb_checksum2 = 0;
        uint64_t this_dumb_checksum3 = 0;
        uint64_t this_dumb_checksum4 = 0;
        uint64_t this_dumb_checksum5 = 0;
        uint8_t this_dumb_checksum6u = 0;
        uint8_t this_dumb_checksum6o = 0;

        setlocale(LC_ALL, "");
        int opt = 0;

        for (int j = 0; j < argc; ++j) {
                frozen[j] = 0;
        }

        while ((opt = getopt(argc, argv, "vhRHs:e:c:m:l:f:")) != -1) {
                switch (opt) {
                case 'v':
                        verbose++;
                        break;
                case 'h':
                        usage(argv[0]);
                        ret = 0;
                        goto done;
                case 'R':
                        require_all_nonfrozen = 1;
                        break;
                case 'H':
                        check_heuristics_only = 1;
                        break;
                case 's':
                        start_file = optarg;
                        break;
                case 'e':
                        end_file = optarg;
                        break;
                case 'c':
                        warn_level = strtoll(optarg, &p, 0);

                        if (p &&
                            *p) {
                                fprintf(stderr, "Invalid edge cap: %s\n",
                                        optarg);
                                ret = EINVAL;
                                goto done;
                        }

                        break;
                case 'l':
                        start_len_l = strtoll(optarg, &p, 0);

                        if ((start_len_l <= 0) ||
                            (p &&
                             *p)) {
                                fprintf(stderr, "Invalid start length: %s\n",
                                        optarg);
                                ret = EINVAL;
                                goto done;
                        }

                        start_len = start_len_l;
                        break;
                case 'm':
                        max_len_l = strtoll(optarg, &p, 0);

                        if ((max_len_l <= 0) ||
                            (p &&
                             *p)) {
                                fprintf(stderr, "Invalid max length: %s\n",
                                        optarg);
                                ret = EINVAL;
                                goto done;
                        }

                        max_len = max_len_l;
                        break;
                case 'f':
                        frozen[fidx] = optarg;
                        fidx++;
                        break;
                default:
                        usage(argv[0]);
                        ret = EINVAL;
                        goto done;
                }
        }

        if (!start_file ||
            !end_file) {
                usage(argv[0]);
                ret = EINVAL;
                goto done;
        }

        size_t num_nonfrozen = 0;

        if (setup(start_file, end_file, frozen, fidx, &q, &qend,
                  &num_nonfrozen)) {
                goto done;
        }

        if (warn_level == (size_t) -1) {
                warn_level = 0;

                for (size_t k = 0; k < q.v_num; ++k) {
                        if (warn_level < q.v[k].fatness) {
                                warn_level = q.v[k].fatness;
                        }
                }

                warn_level = warn_level + 1;
        }

        if (warn_level >= UINT_MAX) {
                warn_level = UINT_MAX - 1;
        }

        if ((ret = quiver_set_warning_threshold(&q, warn_level, &errstr))) {
                fprintf(stderr, "%s\n", errstr);
                goto done;
        }

        /* Save off qend's vertices.  XXX: Why not just USE qend? */
        if ((qend.v_num * qend.v_num) / qend.v_num != qend.v_num) {
                /* Should never happen - isn't qend.v_num capped? */
                ret = errno = EOVERFLOW;
                perror(L(""));
                goto done;
        }

        if (!(target = calloc(qend.v_num * qend.v_num, sizeof *target))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        if (!(seen_in_seq = calloc(qend.v_num + 1, sizeof *seen_in_seq))) {
                ret = errno;
                perror(L("calloc"));
                goto done;
        }

        for (size_t i = 0; i < qend.v_num * qend.v_num; ++i) {
                target[i] = (struct rational) { .p = 0, .q = 1 };
        }

        for (size_t i = 0; i < qend.v_num; ++i) {
                for (size_t j = 0; j < qend.v_num; ++j) {
                        target[i * qend.v_num + j] = qend.e[i * qend.v_len + j];
                }
        }

        if (check_heuristics_only) {
                desired_dumb_checksum1 = dumb_checksum1(&qend);
                desired_dumb_checksum2 = dumb_checksum2(&qend);
                desired_dumb_checksum3 = dumb_checksum3(&qend);
                desired_dumb_checksum4 = dumb_checksum4(&qend, num_nonfrozen);
                desired_dumb_checksum5 = dumb_checksum5(&qend);
                desired_dumb_checksum6o = dumb_checksum6(&qend, 1);
                desired_dumb_checksum6u = dumb_checksum6(&qend, 0);
        }

        /* Now, the long haul*/
        size_t len = 0;

        if (start_len) {
                len = start_len - 1;
        }

        signal(SIGUSR1, handle);
        size_t s = 0;
        size_t c_idx;
        int warn = 0;
        uint_fast8_t correct = 0;
        size_t *sequence = 0;

        while (len < max_len) {
                len++;
                free(sequence);
                sequence = calloc(len + 1, sizeof *sequence);

                for (size_t k = 0; k < len; ++k) {
                        seen_in_seq[k] = 0;
                }

                for (size_t k = 0; k < len; ++k) {
                        /* This will get pruned, it's okay */
                        sequence[k] = 0;
                        seen_in_seq[sequence[k]]++;
                }

                c_idx = 0;
have_full_sequence:

                if (status_requested) {
                        status_requested = 0;
                        printf("Current sequence %s", q.v[sequence[0]].name);

                        for (size_t k = 1; k < len; ++k) {
                                printf(", %s", q.v[sequence[k]].name);
                        }

                        printf("\n");
                        fflush(stdout);
                }

                /*
                 * Here, we have computed a full sequence that we want
                 * to check.  c_idx is in [0, len - 1], and the quiver
                 * has been mutated following the sequence up to (but
                 * not including) c_idx.
                 *
                 * This code is duplicated from increment_sequence
                 * because it is reasonable to start the whole process
                 * with a sequence that may, in fact, produce warnings.
                 */
                warn = 0;

                while (c_idx < len) {
                        /* Prune duplicates (mutations are involutions) */
                        if (c_idx &&
                            sequence[c_idx] == sequence[c_idx - 1]) {
                                goto increment_sequence;
                        }

                        /* Proceed one more mutation */
                        quiver_mutate(&q, sequence[c_idx], &warn, 0);

                        if (warn) {
                                /* That mutation caused a warning: skip it */
                                quiver_mutate(&q, sequence[c_idx], 0, 0);

                                /* And try the sequence, starting at c_idx */
                                goto increment_sequence;
                        }

                        /*
                         * Skip ... X, Y ... when X > Y and the edge
                         * from X to Y is trivial
                         */
                        if (c_idx &&
                            sequence[c_idx] < sequence[c_idx - 1]) {
                                size_t x = sequence[c_idx - 1];
                                size_t y = sequence[c_idx];
                                struct rational *xy = &(q.e[x * q.v_len + y]);
                                struct rational *yx = &(q.e[y * q.v_len + x]);

                                if (!xy->p &&
                                    !yx->p) {
                                        /* Undo mutation */
                                        quiver_mutate(&q, sequence[c_idx], 0,
                                                      0);

                                        /* And try the sequence, starting at c_idx */
                                        goto increment_sequence;
                                }
                        }

                        c_idx++;
                }

                /*
                 * If we are looking for a flip, we're probably
                 * going to require that all non-frozen vertices
                 * participate in the mutation.
                 */
                if (require_all_nonfrozen) {
                        for (size_t k = 0; k < num_nonfrozen; ++k) {
                                if (!seen_in_seq[k]) {
                                        goto finished_check;
                                }
                        }
                }

                /*
                 * Here, there's a full sequence, and the quiver has
                 * been mutated according to it (all of it).  We want
                 * to check if the result matches the target.
                 */
                if (check_heuristics_only) {
                        /* Take a short-cut */
                        this_dumb_checksum1 = dumb_checksum1(&q);

                        if (this_dumb_checksum1 != desired_dumb_checksum1) {
                                goto finished_check;
                        }

                        this_dumb_checksum2 = dumb_checksum2(&q);

                        if (this_dumb_checksum2 != desired_dumb_checksum2) {
                                goto finished_check;
                        }

                        this_dumb_checksum3 = dumb_checksum3(&q);

                        if (this_dumb_checksum3 != desired_dumb_checksum3) {
                                goto finished_check;
                        }

                        this_dumb_checksum4 = dumb_checksum4(&q, num_nonfrozen);

                        if (this_dumb_checksum4 != desired_dumb_checksum4) {
                                goto finished_check;
                        }

                        this_dumb_checksum5 = dumb_checksum5(&q);

                        if (this_dumb_checksum5 != desired_dumb_checksum5) {
                                goto finished_check;
                        }

                        this_dumb_checksum6o = dumb_checksum6(&q, 1);

                        if (this_dumb_checksum6o != desired_dumb_checksum6o) {
                                goto finished_check;
                        }

                        this_dumb_checksum6u = dumb_checksum6(&q, 0);

                        if (this_dumb_checksum6u != desired_dumb_checksum6u) {
                                goto finished_check;
                        }

                        goto print_this;
                }

                correct = 1;

                for (size_t i = 0; i < q.v_num; ++i) {
                        for (size_t j = 0; j < q.v_num; ++j) {
                                struct rational *e = &(q.e[i * q.v_len + j]);
                                struct rational *f = &(target[i * q.v_num + j]);

                                if (e->p * f->q != f->p * e->q) {
                                        correct = 0;
                                        i = q.v_num + 1;
                                        j = i;
                                }
                        }
                }

print_this:

                if (correct ||
                    check_heuristics_only) {
                        printf("SEQUENCE:     %s", q.v[sequence[0]].name);

                        for (size_t k = 1; k < len; ++k) {
                                printf(", %s", q.v[sequence[k]].name);
                        }

                        printf("\n");
                        fflush(stdout);
                }

finished_check:
                c_idx = len - 1;
                quiver_mutate(&q, sequence[c_idx], 0, 0);
increment_sequence:

                /*
                 * Here, we have some kind of sequence. The quiver agrees
                 * with it up to (but not including) c_idx, and we want
                 * to increment it at c_idx.
                 */
                s++;
                seen_in_seq[sequence[c_idx]]--;
                sequence[c_idx]++;
                seen_in_seq[sequence[c_idx]]++;

                if (verbose &&
                    c_idx == 2) {
                        double completed = sequence[0] * num_nonfrozen *
                                           num_nonfrozen + sequence[1] *
                                           num_nonfrozen +
                                           sequence[2] - 1;
                        double total = num_nonfrozen * num_nonfrozen *
                                       num_nonfrozen;

                        printf("Length %zu: completed %.3g%%\n", len, (100.0 *
                                                                       completed)
                               / total);
                        fflush(stdout);
                }

                /* Skip over repeats: each mutation is an involution */
                if (c_idx > 0 &&
                    sequence[c_idx] == sequence[c_idx - 1]) {
                        seen_in_seq[sequence[c_idx]]--;
                        sequence[c_idx]++;
                        seen_in_seq[sequence[c_idx]]++;
                }

                /* Wrap back to 0 */
                if (sequence[c_idx] >= num_nonfrozen) {
                        if (!c_idx) {
                                goto exhausted_len;
                        }

                        c_idx--;
                        quiver_mutate(&q, sequence[c_idx], 0, 0);
                        goto increment_sequence;
                }

                if (c_idx != len - 1) {
                        /*
                         * Here, we now have a start of a sequence.
                         * The quiver agrees with our start up to (but
                         * not including) c_idx, and we want to possibly
                         * complete this to a full sequence. First,
                         * however, we need to make sure that we don't
                         * introduce any warnings at c_idx.
                         *
                         * This single-mutation testing block could be
                         * ignored, but I have a vague suspicion it will
                         * be a hot path.
                         */
                        warn = 0;
                        quiver_mutate(&q, sequence[c_idx], &warn, 0);

                        if (warn) {
                                /* That mutation caused a warning: skip it */
                                quiver_mutate(&q, sequence[c_idx], 0, 0);

                                /* And try the sequence, starting at c_idx */
                                goto increment_sequence;
                        }

                        /*
                         * Quivers almost commute: in particular if
                         * the edge v_i -> v_j is zero, then mutation
                         * at v_i and v_j commute. So if we're at X, Y
                         * in the list, and X > Y, then we've actually
                         * already computed that branch of the tree,
                         * back when it was Y, X. Therefore, increment X.
                         */
                        if (c_idx &&
                            sequence[c_idx] < sequence[c_idx - 1]) {
                                size_t x = sequence[c_idx - 1];
                                size_t y = sequence[c_idx];
                                struct rational *xy = &(q.e[x * q.v_len + y]);
                                struct rational *yx = &(q.e[y * q.v_len + x]);

                                if (!xy->p &&
                                    !yx->p) {
                                        /* Undo mutation */
                                        quiver_mutate(&q, sequence[c_idx], 0,
                                                      0);

                                        /* And try the sequence, starting at c_idx */
                                        goto increment_sequence;
                                }
                        }

                        /*
                         * Here the quiver agrees with our start up
                         * to and including c_idx, and we need to fill
                         * it out into a full sequence so that it can
                         * be tested.  We don't need to perform all the
                         * mutations of that full sequence, but we do
                         * need to generate the minimal one.
                         */
                        for (size_t i = c_idx + 1; i < len; ++i) {
                                seen_in_seq[sequence[i]]--;
                                sequence[i] = 1 - !!((i - c_idx) % 2);
                                seen_in_seq[sequence[i]]++;
                        }

                        c_idx++;
                }

                goto have_full_sequence;
exhausted_len:

                /*
                 * Here, we have c_idx = 0, and the quiver has not been
                 * mutated by anything in sequence
                 */
                printf("Exhausted length %zu\n", len);
        }

        free(sequence);
        fflush(stdout);
done:
        free(target);
        free(seen_in_seq);

        return ret;
}