/* Independent exhaustive zero-forcing search for P(n,k) by FIRST-FORCE
 * NORMALISATION (Stage 5 robustness check; written independently of
 * code/brute_force.c: different reduction principle, different closure
 * algorithm, different enumeration scheme).
 *
 * Principle (Lemma 3.5 argument, results/derivations.md): if S is a
 * zero-forcing set with |S| = s < 2n, some first force x -> y exists, so
 * x in S, y not in S, and the two other neighbours of x are in S.  The
 * directed edges of P(n,k) fall into four orbits under the dihedral group
 * D_n (which acts by automorphisms), with representatives
 *   u0 -> v0,  u0 -> u1,  v0 -> u0,  v0 -> v_k.
 * Hence a zero-forcing s-set exists iff for one of the four representatives
 * (x,y) some s-set S with {x} u (N(x) \ {y}) subset S, y not in S, forces.
 * We therefore enumerate all (s-3)-subsets T of V \ (N[x]) and test
 * S = {x} u (N(x)\{y}) u T.  Total work 4 * C(2n-4, s-3).
 *
 * Closure: counter/worklist algorithm (each black vertex keeps the number
 * of its white neighbours; a black vertex with count 1 is pushed and
 * performs its force), NOT the pass-based bitmask loop of brute_force.c.
 *
 * usage: ff_brute n k s [threads]
 * output: one line per representative and a total line; a witness is
 * printed for the first forcing set found (if any).
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <pthread.h>
#include <stdint.h>

static int n, k, N, s;
static int nbr[128][3];             /* neighbour lists */
static int rep_x[4], rep_y[4];

static void build_graph(void) {
    N = 2 * n;
    for (int i = 0; i < n; i++) {
        int u = i, v = n + i;
        nbr[u][0] = (i + n - 1) % n;  nbr[u][1] = (i + 1) % n;  nbr[u][2] = v;
        nbr[v][0] = n + (i + n - k) % n; nbr[v][1] = n + (i + k) % n; nbr[v][2] = u;
    }
    /* sanity: simple cubic graph, symmetric adjacency */
    for (int w = 0; w < N; w++) {
        for (int a = 0; a < 3; a++) {
            int t = nbr[w][a];
            if (t == w) { fprintf(stderr, "loop at %d\n", w); exit(3); }
            for (int b = 0; b < a; b++) if (nbr[w][b] == t) { fprintf(stderr, "multi-edge at %d\n", w); exit(3); }
            int back = 0; for (int b = 0; b < 3; b++) if (nbr[t][b] == w) back = 1;
            if (!back) { fprintf(stderr, "asymmetric adjacency %d %d\n", w, t); exit(3); }
        }
    }
    rep_x[0] = 0;  rep_y[0] = n;            /* u0 -> v0 */
    rep_x[1] = 0;  rep_y[1] = 1;            /* u0 -> u1 */
    rep_x[2] = n;  rep_y[2] = 0;            /* v0 -> u0 */
    rep_x[3] = n;  rep_y[3] = n + k;        /* v0 -> v_k */
}

/* returns 1 iff the closure of the set given by black[] is V */
static int forces_all(unsigned char *black) {
    int cnt[128], queue[128], qh = 0, qt = 0, nblack = 0;
    for (int w = 0; w < N; w++) if (black[w]) nblack++;
    for (int w = 0; w < N; w++) {
        if (!black[w]) continue;
        int c = 0;
        for (int a = 0; a < 3; a++) if (!black[nbr[w][a]]) c++;
        cnt[w] = c;
        if (c == 1) queue[qt++] = w;
    }
    while (qh < qt) {
        int b = queue[qh++];
        if (cnt[b] != 1) continue;           /* stale entry */
        int y = -1;
        for (int a = 0; a < 3; a++) if (!black[nbr[b][a]]) { y = nbr[b][a]; break; }
        if (y < 0) continue;
        black[y] = 1; nblack++;
        if (nblack == N) return 1;
        /* y became black: update counters of y's black neighbours, and y's own counter */
        int cy = 0;
        for (int a = 0; a < 3; a++) {
            int t = nbr[y][a];
            if (black[t]) { cnt[t]--; if (cnt[t] == 1) queue[qt++] = t; }
            else cy++;
        }
        cnt[y] = cy;
        if (cy == 1) queue[qt++] = y;
    }
    return nblack == N;
}

/* work items: (representative r, first element index) */
static int nthreads = 1;
static volatile long next_item = 0;
static pthread_mutex_t mu = PTHREAD_MUTEX_INITIALIZER;
static unsigned long long tot_enum[4], tot_forcing[4];
static int witness_found = 0; static unsigned char witness[128];

typedef struct { int id; } targ;

static void *worker(void *arg) {
    (void)arg;
    int free_list[128], nfree;
    int idx[128];
    unsigned char base[128], black[128];
    unsigned long long enumerated = 0, forcing = 0;
    for (;;) {
        pthread_mutex_lock(&mu); long item = next_item++; pthread_mutex_unlock(&mu);
        int r = (int)(item / N), first = (int)(item % N);
        if (r >= 4) break;
        int x = rep_x[r], y = rep_y[r];
        memset(base, 0, N);
        base[x] = 1; for (int a = 0; a < 3; a++) if (nbr[x][a] != y) base[nbr[x][a]] = 1;
        nfree = 0;
        for (int w = 0; w < N; w++) if (!base[w] && w != y) free_list[nfree++] = w;
        int t = s - 3;
        if (t < 0) continue;
        /* enumerate t-subsets of free_list whose smallest index is `first` */
        if (first >= nfree) continue;
        if (t == 0) { if (first != 0) continue; memcpy(black, base, N); enumerated++; if (forces_all(black)) { forcing++; pthread_mutex_lock(&mu); if (!witness_found) { witness_found = 1; memcpy(witness, base, N); } pthread_mutex_unlock(&mu); } continue; }
        if (nfree - first < t) continue;
        idx[0] = first; for (int i = 1; i < t; i++) idx[i] = first + i;
        for (;;) {
            memcpy(black, base, N);
            for (int i = 0; i < t; i++) black[free_list[idx[i]]] = 1;
            enumerated++;
            if (forces_all(black)) {
                forcing++;
                pthread_mutex_lock(&mu);
                if (!witness_found) { witness_found = 1; memcpy(witness, base, N); for (int i = 0; i < t; i++) witness[free_list[idx[i]]] = 1; }
                pthread_mutex_unlock(&mu);
            }
            /* next combination with idx[0] fixed */
            int i = t - 1;
            while (i >= 1 && idx[i] == nfree - t + i) i--;
            if (i < 1) break;
            idx[i]++;
            for (int j = i + 1; j < t; j++) idx[j] = idx[j - 1] + 1;
        }
        pthread_mutex_lock(&mu); tot_enum[r] += enumerated; tot_forcing[r] += forcing; pthread_mutex_unlock(&mu);
        enumerated = forcing = 0;
    }
    return NULL;
}

int main(int argc, char **argv) {
    if (argc < 4) { fprintf(stderr, "usage: ff_brute n k s [threads]\n"); return 2; }
    n = atoi(argv[1]); k = atoi(argv[2]); s = atoi(argv[3]);
    nthreads = argc > 4 ? atoi(argv[4]) : 1;
    if (n > 64 || n < 2 * k + 1) { fprintf(stderr, "need 2k+1 <= n <= 64\n"); return 2; }
    build_graph();
    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);
    const char *names[4] = {"u0->v0", "u0->u1", "v0->u0", "v0->vk"};
    unsigned long long E = 0, F = 0;
    for (int r = 0; r < 4; r++) { printf("n=%d k=%d s=%d rep=%s enumerated=%llu forcing=%llu\n", n, k, s, names[r], tot_enum[r], tot_forcing[r]); E += tot_enum[r]; F += tot_forcing[r]; }
    printf("n=%d k=%d s=%d TOTAL enumerated=%llu forcing=%llu", n, k, s, E, F);
    if (witness_found) { printf(" witness="); int first = 1; for (int w = 0; w < N; w++) if (witness[w]) { printf("%s%c%d", first ? "" : ",", w < n ? 'u' : 'v', w < n ? w : w - n); first = 0; } }
    printf("\n");
    return 0;
}