/* Exhaustive boundary minimisation for P(n,k) (Stage 5 robustness check).
 *
 * For every p in [p_lo, p_hi] enumerates ALL p-subsets X of V(P(n,k))
 * (no symmetry reduction, no SAT) and tabulates the distribution of
 * |dX| = |N[X] \ X|.  Output per p: beta_p(n) = min |dX| and the histogram.
 * Independent of code/sat_or_ilp_solver.py (Lemma 3.4 / F(n,p)).
 *
 * usage: beta_brute n k p_lo p_hi [threads]      (2n <= 64)
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <pthread.h>

typedef uint64_t mask;
static int n, k, N;
static mask nb[64];
static int plo, phi, nthreads;
static unsigned long long hist[64][65];
static pthread_mutex_t mu = PTHREAD_MUTEX_INITIALIZER;
static volatile long next_item = 0;

static void build(void) {
    N = 2 * n;
    memset(nb, 0, sizeof nb);
    for (int i = 0; i < n; i++) {
        int u = i, v = n + i;
        nb[u] |= (mask)1 << ((i + n - 1) % n); nb[u] |= (mask)1 << ((i + 1) % n); nb[u] |= (mask)1 << v;
        nb[v] |= (mask)1 << (n + (i + n - k) % n); nb[v] |= (mask)1 << (n + (i + k) % n); nb[v] |= (mask)1 << u;
    }
    for (int w = 0; w < N; w++) if (__builtin_popcountll(nb[w]) != 3) { fprintf(stderr, "degree != 3\n"); exit(3); }
}

static void *worker(void *a) {
    (void)a;
    unsigned long long local[64][65]; memset(local, 0, sizeof local);
    int idx[64];
    for (;;) {
        pthread_mutex_lock(&mu); long item = next_item++; pthread_mutex_unlock(&mu);
        int p = plo + (int)(item / N), first = (int)(item % N);
        if (p > phi) break;
        if (p == 0) { if (first == 0) local[p][0]++; continue; }
        if (N - first < p) continue;
        for (int i = 0; i < p; i++) idx[i] = first + i;
        for (;;) {
            mask X = 0, NX = 0;
            for (int i = 0; i < p; i++) { X |= (mask)1 << idx[i]; NX |= nb[idx[i]]; }
            local[p][__builtin_popcountll(NX & ~X)]++;
            int i = p - 1;
            while (i >= 1 && idx[i] == N - p + i) i--;
            if (i < 1) break;
            idx[i]++;
            for (int j = i + 1; j < p; j++) idx[j] = idx[j - 1] + 1;
        }
    }
    pthread_mutex_lock(&mu);
    for (int p = 0; p < 64; p++) for (int b = 0; b <= 64; b++) hist[p][b] += local[p][b];
    pthread_mutex_unlock(&mu);
    return NULL;
}

int main(int argc, char **argv) {
    if (argc < 5) { fprintf(stderr, "usage: beta_brute n k p_lo p_hi [threads]\n"); return 2; }
    n = atoi(argv[1]); k = atoi(argv[2]); plo = atoi(argv[3]); phi = atoi(argv[4]);
    nthreads = argc > 5 ? atoi(argv[5]) : 1;
    if (2 * n > 64 || n < 2 * k + 1 || phi >= 64) { fprintf(stderr, "bad arguments\n"); return 2; }
    build();
    pthread_t th[64];
    for (int i = 0; i < nthreads; i++) pthread_create(&th[i], NULL, worker, NULL);
    for (int i = 0; i < nthreads; i++) pthread_join(th[i], NULL);
    for (int p = plo; p <= phi; p++) {
        unsigned long long total = 0; int beta = -1;
        for (int b = 0; b <= 64; b++) { total += hist[p][b]; if (hist[p][b] && beta < 0) beta = b; }
        printf("n=%d k=%d p=%d subsets=%llu beta=%d hist=", n, k, p, total, beta);
        for (int b = 0; b <= 64; b++) if (hist[p][b]) printf("%d:%llu ", b, hist[p][b]);
        printf("\n");
    }
    return 0;
}