#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <stdint.h>

#include <math.h>

typedef struct {

uint32_t state[8];

uint64_t bitlen;

unsigned char buf[64];

size_t buflen;

} Sha256;

static const uint32_t sha256_k[64] = {

0x428a2f98,0x71374491,0xb5c0fbcf,0xe9b5dba5,0x3956c25b,0x59f111f1,0x923f82a4,0xab1c5ed5,

0xd807aa98,0x12835b01,0x243185be,0x550c7dc3,0x72be5d74,0x80deb1fe,0x9bdc06a7,0xc19bf174,

0xe49b69c1,0xefbe4786,0x0fc19dc6,0x240ca1cc,0x2de92c6f,0x4a7484aa,0x5cb0a9dc,0x76f988da,

0x983e5152,0xa831c66d,0xb00327c8,0xbf597fc7,0xc6e00bf3,0xd5a79147,0x06ca6351,0x14292967,

0x27b70a85,0x2e1b2138,0x4d2c6dfc,0x53380d13,0x650a7354,0x766a0abb,0x81c2c92e,0x92722c85,

0xa2bfe8a1,0xa81a664b,0xc24b8b70,0xc76c51a3,0xd192e819,0xd6990624,0xf40e3585,0x106aa070,

0x19a4c116,0x1e376c08,0x2748774c,0x34b0bcb5,0x391c0cb3,0x4ed8aa4a,0x5b9cca4f,0x682e6ff3,

0x748f82ee,0x78a5636f,0x84c87814,0x8cc70208,0x90befffa,0xa4506ceb,0xbef9a3f7,0xc67178f2

};

static uint32_t sha256_rotr(uint32_t x, int n) { return (x >> n) | (x << (32 - n)); }

static void sha256_transform(Sha256 *ctx, const unsigned char *data) {

uint32_t m[64], a, b, c, d, e, f, g, h, t1, t2;

int i;

for (i = 0; i < 16; i++)

m[i] = ((uint32_t)data[i*4] << 24) | ((uint32_t)data[i*4+1] << 16) |

((uint32_t)data[i*4+2] << 8) | ((uint32_t)data[i*4+3]);

for (i = 16; i < 64; i++) {

uint32_t s0 = sha256_rotr(m[i-15], 7) ^ sha256_rotr(m[i-15], 18) ^ (m[i-15] >> 3);

uint32_t s1 = sha256_rotr(m[i-2], 17) ^ sha256_rotr(m[i-2], 19) ^ (m[i-2] >> 10);

m[i] = m[i-16] + s0 + m[i-7] + s1;

}

a = ctx->state[0]; b = ctx->state[1]; c = ctx->state[2]; d = ctx->state[3];

e = ctx->state[4]; f = ctx->state[5]; g = ctx->state[6]; h = ctx->state[7];

for (i = 0; i < 64; i++) {

uint32_t S1 = sha256_rotr(e,6) ^ sha256_rotr(e,11) ^ sha256_rotr(e,25);

uint32_t ch = (e & f) ^ ((~e) & g);

t1 = h + S1 + ch + sha256_k[i] + m[i];

uint32_t S0 = sha256_rotr(a,2) ^ sha256_rotr(a,13) ^ sha256_rotr(a,22);

uint32_t maj = (a & b) ^ (a & c) ^ (b & c);

t2 = S0 + maj;

h = g; g = f; f = e; e = d + t1;

d = c; c = b; b = a; a = t1 + t2;

}

ctx->state[0] += a; ctx->state[1] += b; ctx->state[2] += c; ctx->state[3] += d;

ctx->state[4] += e; ctx->state[5] += f; ctx->state[6] += g; ctx->state[7] += h;

}

static void sha256_init(Sha256 *ctx) {

ctx->state[0]=0x6a09e667; ctx->state[1]=0xbb67ae85; ctx->state[2]=0x3c6ef372; ctx->state[3]=0xa54ff53a;

ctx->state[4]=0x510e527f; ctx->state[5]=0x9b05688c; ctx->state[6]=0x1f83d9ab; ctx->state[7]=0x5be0cd19;

ctx->bitlen = 0; ctx->buflen = 0;

}

static void sha256_update(Sha256 *ctx, const unsigned char *data, size_t len) {

size_t i;

for (i = 0; i < len; i++) {

ctx->buf[ctx->buflen++] = data[i];

if (ctx->buflen == 64) {

sha256_transform(ctx, ctx->buf);

ctx->bitlen += 512;

ctx->buflen = 0;

}

}

}

static void sha256_final(Sha256 *ctx, unsigned char out[32]) {

size_t i = ctx->buflen;

ctx->bitlen += (uint64_t)ctx->buflen * 8;

ctx->buf[i++] = 0x80;

if (i > 56) {

while (i < 64) ctx->buf[i++] = 0;

sha256_transform(ctx, ctx->buf);

i = 0;

}

while (i < 56) ctx->buf[i++] = 0;

for (int k = 7; k >= 0; k--) ctx->buf[i++] = (unsigned char)(ctx->bitlen >> (k*8));

sha256_transform(ctx, ctx->buf);

for (i = 0; i < 8; i++) {

out[i*4] = (unsigned char)(ctx->state[i] >> 24);

out[i*4+1] = (unsigned char)(ctx->state[i] >> 16);

out[i*4+2] = (unsigned char)(ctx->state[i] >> 8);

out[i*4+3] = (unsigned char)(ctx->state[i]);

}

}

static void sha256_hex(const unsigned char *data, long len, char out_hex[65]) {

Sha256 ctx; unsigned char digest[32];

sha256_init(&ctx);

sha256_update(&ctx, data, (size_t)len);

sha256_final(&ctx, digest);

for (int i = 0; i < 32; i++) sprintf(out_hex + i*2, "%02x", digest[i]);

out_hex[64] = 0;

}

static uint64_t rng_next(uint64_t *state) {

uint64_t x = *state;

x ^= x << 13;

x ^= x >> 7;

x ^= x << 17;

*state = x;

return x;

}

static unsigned char *read_file(const char *path, long *out_len) {

FILE *f = fopen(path, "rb");

if (!f) { fprintf(stderr, "cannot open %s\n", path); exit(1); }

fseek(f, 0, SEEK_END);

long n = ftell(f);

fseek(f, 0, SEEK_SET);

unsigned char *buf = malloc((size_t)n);

if (fread(buf, 1, (size_t)n, f) != (size_t)n) { fprintf(stderr, "short read %s\n", path); exit(1); }

fclose(f);

*out_len = n;

return buf;

}

static void build_cipher_perm(unsigned char perm[256], uint64_t seed) {

uint64_t st = seed;

for (int i = 0; i < 256; i++) perm[i] = (unsigned char)i;

for (int i = 255; i >= 1; i--) {

uint64_t r = rng_next(&st);

int j = (int)(r % (uint64_t)(i + 1));

unsigned char t = perm[i]; perm[i] = perm[j]; perm[j] = t;

}

}

static unsigned char *build_w_iso(const unsigned char *A, long n, const unsigned char perm[256]) {

unsigned char *out = malloc((size_t)n);

for (long i = 0; i < n; i++) out[i] = perm[A[i]];

return out;

}

static unsigned char *build_w_ghost(const unsigned char *A_train, long train_len, long out_n, uint64_t seed) {

long counts[256] = {0};

for (long i = 0; i < train_len; i++) counts[A_train[i]]++;

long cum[257]; cum[0] = 0;

for (int b = 0; b < 256; b++) cum[b+1] = cum[b] + counts[b];

long total = train_len;

unsigned char *out = malloc((size_t)out_n);

uint64_t st = seed;

for (long i = 0; i < out_n; i++) {

uint64_t r = rng_next(&st);

long v = (long)(r % (uint64_t)total);

int lo = 0, hi = 255, b = 0;

/* binary search: cum[b] <= v < cum[b+1] */

while (lo <= hi) {

int mid = (lo + hi) / 2;

if (v < cum[mid+1]) { b = mid; hi = mid - 1; }

else lo = mid + 1;

}

out[i] = (unsigned char)b;

}

return out;

}

typedef struct { unsigned char word[32]; int len; long count; long first_pos; } WordFreq;

static int is_ws(unsigned char c) {

return c==' '||c=='\t'||c=='\n'||c=='\r'||c=='\v'||c=='\f';

}

static unsigned char *build_w_ff(const unsigned char *A, long n, long train_len,

long *out_n,

unsigned char rank_word[16][32], int rank_len[16]) {

/* tokenize train region by whitespace, count frequency + first occurrence word-index */

WordFreq *freqs = NULL; long nfreq = 0, cap = 0;

long widx = 0;

long i = 0;

while (i < train_len) {

if (is_ws(A[i])) { i++; continue; }

long j = i;

while (j < train_len && !is_ws(A[j])) j++;

int wlen = (int)(j - i);

if (wlen > 31) wlen = 31; /* clamp; top words are short in practice */

/* find in freqs */

long found = -1;

for (long k = 0; k < nfreq; k++) {

if (freqs[k].len == wlen && memcmp(freqs[k].word, A + i, (size_t)wlen) == 0) { found = k; break; }

}

if (found >= 0) {

freqs[found].count++;

} else {

if (nfreq == cap) {

cap = cap ? cap * 2 : 1024;

freqs = realloc(freqs, (size_t)cap * sizeof(WordFreq));

}

memcpy(freqs[nfreq].word, A + i, (size_t)wlen);

freqs[nfreq].len = wlen;

freqs[nfreq].count = 1;

freqs[nfreq].first_pos = widx;

nfreq++;

}

widx++;

i = j;

}

/* sort by count desc, tie-break first_pos asc (simple insertion since top16 needed, but do full sort for determinism) */

for (long a = 1; a < nfreq; a++) {

WordFreq key = freqs[a];

long b = a - 1;

while (b >= 0 && (freqs[b].count < key.count ||

(freqs[b].count == key.count && freqs[b].first_pos > key.first_pos))) {

freqs[b+1] = freqs[b];

b--;

}

freqs[b+1] = key;

}

int top = nfreq < 16 ? (int)nfreq : 16;

if (top < 16) { fprintf(stderr, "FATAL: fewer than 16 distinct words in A-train\n"); exit(1); }

for (int k = 0; k < 16; k++) {

rank_len[k] = freqs[k].len;

memcpy(rank_word[k], freqs[k].word, (size_t)freqs[k].len);

}

free(freqs);

/* build swap: pairs (0,1)(2,3)... */

/* apply to full A: tokenize by whitespace, replace matching words */

unsigned char *out = malloc((size_t)n * 4 + 16); /* generous upper bound */

long oi = 0;

i = 0;

while (i < n) {

if (is_ws(A[i])) { out[oi++] = A[i]; i++; continue; }

long j = i;

while (j < n && !is_ws(A[j])) j++;

int wlen = (int)(j - i);

int matched = -1;

if (wlen <= 31) {

for (int k = 0; k < 16; k++) {

if (rank_len[k] == wlen && memcmp(rank_word[k], A + i, (size_t)wlen) == 0) { matched = k; break; }

}

}

if (matched >= 0) {

int partner = matched ^ 1; /* 0<->1, 2<->3, ... */

memcpy(out + oi, rank_word[partner], (size_t)rank_len[partner]);

oi += rank_len[partner];

} else {

memcpy(out + oi, A + i, (size_t)wlen);

oi += wlen;

}

i = j;

}

*out_n = oi;

return out;

}

typedef struct { int32_t l, r; } Merge;

typedef struct {

int32_t *tok;

long n;

} TokStream;

static TokStream grow_bpe(const int32_t *init, long n0, int max_merges, int min_pair,

Merge **merges_out, int *n_merges_out) {

int32_t *tok = malloc((size_t)(n0 ? n0 : 1) * sizeof(int32_t));

memcpy(tok, init, (size_t)n0 * sizeof(int32_t));

int32_t *scratch = malloc((size_t)(n0 ? n0 : 1) * sizeof(int32_t));

long n = n0;

Merge *merges = malloc((size_t)max_merges * sizeof(Merge));

int nm = 0;

int32_t next_id = 256;

/* generation-stamped hashmap for pair -> (count, first_pos): avoids

long cap = 1;

while (cap < (n0 * 2 + 16)) cap <<= 1;

int64_t *keys = malloc(cap * sizeof(int64_t));

long *counts = malloc(cap * sizeof(long));

long *firstpos = malloc(cap * sizeof(long));

uint32_t *gen = calloc((size_t)cap, sizeof(uint32_t));

uint32_t cur_gen = 0;

for (int m = 0; m < max_merges; m++) {

if (n < 2) break;

cur_gen++;

long best_count = -1, best_first = -1;

int32_t best_l = 0, best_r = 0;

int have_best = 0;

for (long i = 0; i < n - 1; i++) {

int64_t key = ((int64_t)tok[i] << 32) ^ (int64_t)(uint32_t)tok[i+1];

uint64_t h = (uint64_t)key * 0x9E3779B97F4A7C15ULL;

long slot = (long)(h & (uint64_t)(cap - 1));

while (gen[slot] == cur_gen && keys[slot] != key) slot = (slot + 1) & (cap - 1);

if (gen[slot] != cur_gen) { gen[slot] = cur_gen; keys[slot] = key; counts[slot] = 0; firstpos[slot] = i; }

counts[slot]++;

long c = counts[slot], f = firstpos[slot];

if (c > best_count || (c == best_count && f < best_first)) {

best_count = c; best_first = f;

best_l = tok[i]; best_r = tok[i+1];

have_best = 1;

}

}

if (!have_best || best_count < min_pair) break;

/* apply merge left-to-right, non-overlapping, into the scratch

long ni = 0;

long i = 0;

while (i < n) {

if (i < n - 1 && tok[i] == best_l && tok[i+1] == best_r) {

scratch[ni++] = next_id;

i += 2;

} else {

scratch[ni++] = tok[i];

i += 1;

}

}

{ int32_t *tmp = tok; tok = scratch; scratch = tmp; }

n = ni;

merges[nm].l = best_l; merges[nm].r = best_r;

nm++;

next_id++;

}

free(scratch);

free(keys); free(counts); free(firstpos); free(gen);

*merges_out = merges;

*n_merges_out = nm;

TokStream result; result.tok = tok; result.n = n;

return result;

}

static TokStream segment_with_merges(const int32_t *init, long n0, const Merge *merges, int n_merges) {

int32_t *tok = malloc((size_t)n0 * sizeof(int32_t));

memcpy(tok, init, (size_t)n0 * sizeof(int32_t));

long n = n0;

int32_t next_id = 256;

for (int m = 0; m < n_merges; m++) {

int32_t l = merges[m].l, r = merges[m].r;

long ni = 0;

long i = 0;

int any = 0;

for (long k = 0; k < n - 1; k++) if (tok[k]==l && tok[k+1]==r) { any = 1; break; }

if (!any) { next_id++; continue; }

int32_t *nt = malloc((size_t)n * sizeof(int32_t));

while (i < n) {

if (i < n - 1 && tok[i] == l && tok[i+1] == r) {

nt[ni++] = next_id;

i += 2;

} else {

nt[ni++] = tok[i];

i += 1;

}

}

free(tok);

tok = nt;

n = ni;

next_id++;

}

TokStream result; result.tok = tok; result.n = n;

return result;

}

static int32_t *bytes_to_units(const unsigned char *b, long n) {

int32_t *u = malloc((size_t)n * sizeof(int32_t));

for (long i = 0; i < n; i++) u[i] = b[i];

return u;

}

typedef struct {

int32_t a, b;

long c;

} Pair2;

typedef struct {

int32_t a, b, c;

long n;

} Pair3;

static int cmp_pair2(const void *x, const void *y) {

const Pair2 *p = x, *q = y;

if (p->a != q->a) return (p->a < q->a) ? -1 : 1;

return (p->b < q->b) ? -1 : (p->b > q->b) ? 1 : 0;

}

static int cmp_pair3(const void *x, const void *y) {

const Pair3 *p = x, *q = y;

if (p->a != q->a) return (p->a < q->a) ? -1 : 1;

if (p->b != q->b) return (p->b < q->b) ? -1 : 1;

return (p->c < q->c) ? -1 : (p->c > q->c) ? 1 : 0;

}

typedef struct {

long *n1; /* size V, unigram counts */

long V;

Pair2 *fwd; /* sorted by (a,b), RLE counts, distance-1 forward pairs */

long n_fwd;

Pair3 *tri; /* sorted by (a,b,c), RLE counts */

long n_tri;

} CountTable;

static long *build_n1(const int32_t *tok, long n, long V) {

long *n1 = calloc((size_t)V, sizeof(long));

for (long i = 0; i < n; i++) n1[tok[i]]++;

return n1;

}

static Pair2 *rle_pairs2(int32_t *raw_a, int32_t *raw_b, long m, long *out_n) {

Pair2 *p = malloc((size_t)m * sizeof(Pair2));

for (long i = 0; i < m; i++) { p[i].a = raw_a[i]; p[i].b = raw_b[i]; p[i].c = 1; }

qsort(p, (size_t)m, sizeof(Pair2), cmp_pair2);

long w = 0;

for (long i = 0; i < m; i++) {

if (w > 0 && p[w-1].a == p[i].a && p[w-1].b == p[i].b) p[w-1].c++;

else p[w++] = p[i];

}

*out_n = w;

return p;

}

static Pair3 *rle_pairs3(int32_t *ra, int32_t *rb, int32_t *rc, long m, long *out_n) {

Pair3 *p = malloc((size_t)m * sizeof(Pair3));

for (long i = 0; i < m; i++) { p[i].a = ra[i]; p[i].b = rb[i]; p[i].c = rc[i]; p[i].n = 1; }

qsort(p, (size_t)m, sizeof(Pair3), cmp_pair3);

long w = 0;

for (long i = 0; i < m; i++) {

if (w > 0 && p[w-1].a == p[i].a && p[w-1].b == p[i].b && p[w-1].c == p[i].c) p[w-1].n++;

else p[w++] = p[i];

}

*out_n = w;

return p;

}

static CountTable build_counts(const int32_t *tok, long n, long V) {

CountTable ct;

ct.V = V;

ct.n1 = build_n1(tok, n, V);

long m2 = n > 1 ? n - 1 : 0;

int32_t *a2 = malloc((size_t)(m2?m2:1) * sizeof(int32_t));

int32_t *b2 = malloc((size_t)(m2?m2:1) * sizeof(int32_t));

for (long i = 0; i < m2; i++) { a2[i] = tok[i]; b2[i] = tok[i+1]; }

ct.fwd = rle_pairs2(a2, b2, m2, &ct.n_fwd);

free(a2); free(b2);

long m3 = n > 2 ? n - 2 : 0;

int32_t *a3 = malloc((size_t)(m3?m3:1) * sizeof(int32_t));

int32_t *b3 = malloc((size_t)(m3?m3:1) * sizeof(int32_t));

int32_t *c3 = malloc((size_t)(m3?m3:1) * sizeof(int32_t));

for (long i = 0; i < m3; i++) { a3[i] = tok[i]; b3[i] = tok[i+1]; c3[i] = tok[i+2]; }

ct.tri = rle_pairs3(a3, b3, c3, m3, &ct.n_tri);

free(a3); free(b3); free(c3);

return ct;

}

static void free_counts(CountTable *ct) {

free(ct->n1); free(ct->fwd); free(ct->tri);

}

static long fwd_lower_bound(const Pair2 *fwd, long n_fwd, int32_t a) {

long lo = 0, hi = n_fwd;

while (lo < hi) { long mid = (lo+hi)/2; if (fwd[mid].a < a) lo = mid+1; else hi = mid; }

return lo;

}

static long tri_lower_bound(const Pair3 *tri, long n_tri, int32_t a, int32_t b) {

long lo = 0, hi = n_tri;

while (lo < hi) {

long mid = (lo+hi)/2;

Pair3 *p = (Pair3*)&tri[mid];

int less = (p->a < a) || (p->a == a && p->b < b);

if (less) lo = mid+1; else hi = mid;

}

return lo;

}

static void ctx2_lookup(const CountTable *ct, int32_t c1, int32_t u, long *sum_out, long *w_out) {

long i = fwd_lower_bound(ct->fwd, ct->n_fwd, c1);

long sum = 0, w = 0;

while (i < ct->n_fwd && ct->fwd[i].a == c1) {

sum += ct->fwd[i].c;

if (ct->fwd[i].b == u) w = ct->fwd[i].c;

i++;

}

*sum_out = sum; *w_out = w;

}

static void ctx3_lookup(const CountTable *ct, int32_t c2, int32_t c1, int32_t u, long *sum_out, long *w_out) {

long i = tri_lower_bound(ct->tri, ct->n_tri, c2, c1);

long sum = 0, w = 0;

while (i < ct->n_tri && ct->tri[i].a == c2 && ct->tri[i].b == c1) {

sum += ct->tri[i].n;

if (ct->tri[i].c == u) w = ct->tri[i].n;

i++;

}

*sum_out = sum; *w_out = w;

}

#define EPS 0.1

static double P1_of(const CountTable *ct, long n1_total, int32_t u) {

long w = (u >= 0 && u < ct->V) ? ct->n1[u] : 0;

double top = (n1_total > 0) ? (1.0 - EPS) * (double)w / (double)n1_total : 0.0;

return top + EPS * (1.0 / (double)ct->V);

}

static double P2_of(const CountTable *ct, long n1_total, int32_t c1, int32_t u) {

long sum, w;

ctx2_lookup(ct, c1, u, &sum, &w);

if (sum == 0) return P1_of(ct, n1_total, u);

return (1.0 - EPS) * (double)w / (double)sum + EPS * P1_of(ct, n1_total, u);

}

static double P3_of(const CountTable *ct, long n1_total, int32_t c2, int32_t c1, int32_t u) {

long sum, w;

ctx3_lookup(ct, c2, c1, u, &sum, &w);

if (sum == 0) return P2_of(ct, n1_total, c1, u);

return (1.0 - EPS) * (double)w / (double)sum + EPS * P2_of(ct, n1_total, c1, u);

}

static double price_unit(const CountTable *ct, long n1_total, long pos_in_world,

int32_t c2, int32_t c1, int32_t u) {

if (pos_in_world == 0) return P1_of(ct, n1_total, u);

if (pos_in_world == 1) return P2_of(ct, n1_total, c1, u);

return P3_of(ct, n1_total, c2, c1, u);

}

static long *build_unit_len(const Merge *merges, int n_merges) {

long V = 256 + n_merges;

long *ulen = malloc((size_t)V * sizeof(long));

for (int i = 0; i < 256; i++) ulen[i] = 1;

for (int i = 0; i < n_merges; i++) ulen[256+i] = ulen[merges[i].l] + ulen[merges[i].r];

return ulen;

}

static Pair2 *build_dist_pairs(const int32_t *tok, long n, int d, long *out_n) {

long m = (n > d) ? n - d : 0;

int32_t *a = malloc((size_t)(m?m:1) * sizeof(int32_t));

int32_t *b = malloc((size_t)(m?m:1) * sizeof(int32_t));

for (long i = 0; i < m; i++) { a[i] = tok[i]; b[i] = tok[i+d]; }

Pair2 *p = rle_pairs2(a, b, m, out_n);

free(a); free(b);

return p;

}

static Pair2 *build_reverse_pairs(const Pair2 *fwd, long n_fwd) {

Pair2 *rev = malloc((size_t)(n_fwd?n_fwd:1) * sizeof(Pair2));

for (long i = 0; i < n_fwd; i++) { rev[i].a = fwd[i].b; rev[i].b = fwd[i].a; rev[i].c = fwd[i].c; }

qsort(rev, (size_t)n_fwd, sizeof(Pair2), cmp_pair2);

return rev;

}

#define NDESC 12

typedef struct {

long V;

unsigned char *alive; /* 0/1, size V */

double *raw; /* V*NDESC */

double *z; /* V*NDESC, z-scored over alive units only */

} Descriptors;

static void for_each_group(const Pair2 *p, long n, void (*fn)(int32_t u, long start, long len, void *ctx), void *ctx) {

long i = 0;

while (i < n) {

long j = i;

while (j < n && p[j].a == p[i].a) j++;

fn(p[i].a, i, j - i, ctx);

i = j;

}

}

typedef struct {

const Pair2 *fwd1;

double *raw; /* V*NDESC */

long V;

} EntropyCtx;

static void entropy_pass(int32_t u, long start, long len, void *ctx_) {

EntropyCtx *ctx = ctx_;

long total = 0;

for (long k = 0; k < len; k++) total += ctx->fwd1[start+k].c;

if (total == 0 || u < 0 || u >= ctx->V) return;

double H = 0.0, top1 = 0.0, sorted4[4] = {0,0,0,0};

for (long k = 0; k < len; k++) {

double p = (double)ctx->fwd1[start+k].c / (double)total;

H += -p * log2(p);

if (p > top1) top1 = p;

for (int s = 0; s < 4; s++) {

if (p > sorted4[s]) {

for (int t = 3; t > s; t--) sorted4[t] = sorted4[t-1];

sorted4[s] = p;

break;

}

}

}

double top4 = sorted4[0]+sorted4[1]+sorted4[2]+sorted4[3];

ctx->raw[u*NDESC + 2] = H; /* right transition entropy */

ctx->raw[u*NDESC + 4] = log2((double)len); /* log2 right degree */

ctx->raw[u*NDESC + 6] = top1;

ctx->raw[u*NDESC + 7] = top4;

}

typedef struct {

const Pair2 *rev1;

double *raw;

long V;

} LeftCtx;

static void left_entropy_pass(int32_t u, long start, long len, void *ctx_) {

LeftCtx *ctx = ctx_;

long total = 0;

for (long k = 0; k < len; k++) total += ctx->rev1[start+k].c;

if (total == 0 || u < 0 || u >= ctx->V) return;

double H = 0.0;

for (long k = 0; k < len; k++) {

double p = (double)ctx->rev1[start+k].c / (double)total;

H += -p * log2(p);

}

ctx->raw[u*NDESC + 3] = H; /* left transition entropy */

ctx->raw[u*NDESC + 5] = log2((double)len); /* log2 left degree */

}

typedef struct {

const Pair2 *dtab;

double *raw;

int dimidx; /* 8,9,10 for d=1,2,4 */

long V;

} ConcCtx;

static void concentration_pass(int32_t u, long start, long len, void *ctx_) {

ConcCtx *ctx = ctx_;

long total = 0;

for (long k = 0; k < len; k++) total += ctx->dtab[start+k].c;

if (total == 0 || u < 0 || u >= ctx->V) return;

double simpson = 0.0;

for (long k = 0; k < len; k++) {

double frac = (double)ctx->dtab[start+k].c / (double)total;

simpson += frac * frac;

}

ctx->raw[u*NDESC + ctx->dimidx] = simpson;

}

typedef struct {

const Pair2 *fwd1;

double *raw;

long V;

} RefineCtx;

static void refine_pass(int32_t u, long start, long len, void *ctx_) {

RefineCtx *ctx = ctx_;

long total = 0;

for (long k = 0; k < len; k++) total += ctx->fwd1[start+k].c;

if (total == 0 || u < 0 || u >= ctx->V) return;

double acc = 0.0;

for (long k = 0; k < len; k++) {

int32_t v = ctx->fwd1[start+k].b;

double p = (double)ctx->fwd1[start+k].c / (double)total;

double Hv = (v >= 0 && v < ctx->V) ? ctx->raw[v*NDESC + 2] : 0.0;

acc += p * Hv;

}

ctx->raw[u*NDESC + 11] = acc;

}

static Descriptors compute_descriptors(const int32_t *tok, long n, long V, const long *n1, const long *ulen) {

Descriptors d;

d.V = V;

d.alive = calloc((size_t)V, 1);

d.raw = calloc((size_t)V * NDESC, sizeof(double));

d.z = calloc((size_t)V * NDESC, sizeof(double));

long n_fwd1; Pair2 *fwd1 = build_dist_pairs(tok, n, 1, &n_fwd1);

long n_fwd2; Pair2 *fwd2 = build_dist_pairs(tok, n, 2, &n_fwd2);

long n_fwd4; Pair2 *fwd4 = build_dist_pairs(tok, n, 4, &n_fwd4);

Pair2 *rev1 = build_reverse_pairs(fwd1, n_fwd1);

for (long u = 0; u < V; u++) {

d.alive[u] = (n1[u] > 0);

if (n1[u] > 0) d.raw[u*NDESC + 0] = log2((double)n1[u]);

d.raw[u*NDESC + 1] = (double)ulen[u];

}

EntropyCtx ec = { fwd1, d.raw, V };

for_each_group(fwd1, n_fwd1, entropy_pass, &ec);

LeftCtx lc = { rev1, d.raw, V };

for_each_group(rev1, n_fwd1, left_entropy_pass, &lc);

ConcCtx cc1 = { fwd1, d.raw, 8, V };

for_each_group(fwd1, n_fwd1, concentration_pass, &cc1);

ConcCtx cc2 = { fwd2, d.raw, 9, V };

for_each_group(fwd2, n_fwd2, concentration_pass, &cc2);

ConcCtx cc4 = { fwd4, d.raw, 10, V };

for_each_group(fwd4, n_fwd4, concentration_pass, &cc4);

RefineCtx rc = { fwd1, d.raw, V };

for_each_group(fwd1, n_fwd1, refine_pass, &rc);

free(fwd1); free(fwd2); free(fwd4); free(rev1);

/* z-score per dimension over alive units */

for (int k = 0; k < NDESC; k++) {