/*
 * lcg.c -- Lazy Clause Generation Phase 1: order-encoded integer
 *          variables, linear_le propagator, theory-hook plumbing to
 *          the SAT engine.
 *
 * Key data flow:
 *
 *     User: lcg_new_int_var(x, lo=0, hi=9)
 *         -> allocate 9 SAT atoms a[0..8], one per threshold:
 *                a[k] == "x <= lo+k"  for k in 0..8 (note: not lo+9 == hi,
 *                because x <= hi is trivially true and needs no atom).
 *         -> post monotonicity clauses  a[k] -> a[k+1]   (as `~a[k] v a[k+1]`).
 *
 *     User: lcg_post_linear_le(coeff, vars, n, K)
 *         -> remember the constraint in our constraint list.
 *         -> on every theory check, walk it and tighten bounds.
 *
 *     User: lcg_solve()
 *         -> register TheoryHooks with the SAT engine.
 *         -> call sat_solve, which drives the search.
 *
 *     SAT engine, during search:
 *
 *         every time SAT assigns a literal: theory_assign(lit)
 *             -> if lit corresponds to an order-encoding atom of some
 *                variable, update that variable's lo or hi via the
 *                lit->(var, kind, threshold) reverse map and push a
 *                bound-trail entry.
 *
 *         every time SAT unassigns a literal: theory_unassign(lit)
 *             -> pop bound-trail entries until the trail matches.
 *
 *         after each round of BCP: theory_check()
 *             -> propagators run, tighten bounds, generate explanations.
 *             -> if some variable's lo > hi, build the conflict clause
 *                from the bound atoms causing the wipeout and return it.
 *             -> otherwise: enqueue any literals derived but not yet
 *                on the SAT trail (these surface via theory_propagate
 *                on the next callback).
 *
 *     SAT engine: theory_propagate() returns enqueued literals one at
 *     a time, each with its explanation clause.
 *
 * Implementation notes:
 *
 *   - The reverse map (lit -> (var, threshold)) is stored as two parallel
 *     arrays indexed by SAT atom id. Negative literals never appear as
 *     order-encoding atoms (the atoms are always positive in their lit
 *     form); the negation of "x <= k" is "x > k" which the theory
 *     handles symmetrically via the same atom with sign flipped.
 *
 *   - The bound trail mirrors the EUF/LRA convention from the SMT code:
 *     a list of (var, side, old_value) records that grows on every
 *     bound tightening and shrinks on theory_unassign.
 *
 *   - Phase 1 has no propagation-queue between propagators; we run
 *     them in a single pass per check, looping until fixpoint. That
 *     is fine for the example workload.
 */
#define _POSIX_C_SOURCE 200809L
#include "lcg.h"
#include "sat.h"

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

/* ====== variable record ====== */
typedef struct {
    char *name;
    int   dom_lo, dom_hi;     /* original (immutable) domain bounds */
    int   cur_lo, cur_hi;     /* current bounds, updated by assign/unassign */
    /* atoms[k] is the SAT atom (positive int) for the literal "x <= dom_lo + k",
     * for k in 0 .. (dom_hi - dom_lo - 1). */
    int  *atoms;
    int   n_atoms;            /* == dom_hi - dom_lo */
} VarRec;

/* ====== bound trail entry ====== */
typedef struct {
    int var;
    int side;                 /* 0 = lo was raised; 1 = hi was lowered */
    int old_value;
} BoundUndo;

/* ====== reverse map entry ====== */
typedef struct {
    int var;                  /* -1 if this atom is not an order-encoding atom */
    int threshold_idx;        /* k in 0 .. (n_atoms - 1) */
} AtomBack;

/* ====== solver state ====== */
struct LcgCtx {
    Sat *sat;

    /* variables */
    VarRec *vars;
    int     n_vars, cap_vars;

    /* atom reverse map: indexed by SAT atom id (1..n_atoms) */
    AtomBack *atom_back;
    int       cap_atom_back;
    int       n_sat_atoms;    /* mirror of sat allocation count */

    /* bound trail */
    BoundUndo *btrail;
    int        n_btrail, cap_btrail;

    /* SAT trail mirror: for every theory_assign, record the bound-trail
     * height just before the assign happened. theory_unassign rewinds. */
    int *pre_btrail;
    int *trail_lits;            /* lits in order of assignment */
    int  n_trail, cap_trail;

    /* Propagation queue. Each entry holds a derived literal and an
     * explanation clause. theory_propagate drains this. */
    struct PropEntry {
        int  lit;
        int *reason;
        int  n_reason;
        int  reason_cap;
    } *propq;
    int n_propq, cap_propq, propq_head;

    /* Conflict scratch space. */
    int *conf_buf;
    int  cap_conf;

    /* Solved? (one-shot like the SMT solver) */
    int solved;

    /* Theory hooks owned by this context (so multiple LcgCtx coexist). */
    TheoryHooks hooks;

    /* Stats. */
    long n_explanations;

    /* Buffered clauses, held until lcg_solve creates the SAT engine.
     * We can't create the SAT engine eagerly because we need to know
     * the final atom count first -- pre-allocating a generous budget
     * leaves unused atoms on the VSIDS heap that the solver still
     * decides on, badly inflating the decision count. Each buffered
     * clause is owned (free on lcg_free / after replay). */
    struct ClauseBuf {
        int *lits;
        int  n_lits;
    } *cl_buf;
    int cl_buf_n, cl_buf_cap;

    /* Cached verbose level (applied when the SAT engine is created). */
    int verbose;

    /* ----- Phase 2 ----- */

    /* Registered propagators. Each propagator is a black-box; theory_check
     * iterates the list. Phase 1's linear_le constraints live here too,
     * wrapped as LcgPropagator instances. */
    LcgPropagator *props;
    int            n_props, cap_props;

    /* Scratch explanation buffer. Populated by lcg_expl_* during a
     * propagator's run; consumed by lcg_emit_*. Only one propagator
     * uses it at a time, and the LCG core clears it before each call.
     *
     * Contents: negated supports, ready to drop into a SAT clause.
     * If the propagator emits a derived literal, the final clause is
     * [derived_lit, neg_support_1, neg_support_2, ...].
     * For a conflict, the clause is just the neg_supports. */
    int *expl_neg;
    int  expl_n, expl_cap;
};

/* ====== bookkeeping helpers ====== */
static void grow_vars(LcgCtx *ctx, int needed) {
    if (needed <= ctx->cap_vars) return;
    int nc = ctx->cap_vars ? ctx->cap_vars * 2 : 8;
    while (nc < needed) nc *= 2;
    ctx->vars = realloc(ctx->vars, (size_t)nc * sizeof(VarRec));
    ctx->cap_vars = nc;
}

static void grow_atom_back(LcgCtx *ctx, int max_atom) {
    if (max_atom < ctx->cap_atom_back) return;
    int nc = ctx->cap_atom_back ? ctx->cap_atom_back * 2 : 32;
    while (nc <= max_atom) nc *= 2;
    ctx->atom_back = realloc(ctx->atom_back, (size_t)nc * sizeof(AtomBack));
    for (int i = ctx->cap_atom_back; i < nc; i++) {
        ctx->atom_back[i].var = -1;
        ctx->atom_back[i].threshold_idx = -1;
    }
    ctx->cap_atom_back = nc;
}

static void grow_btrail(LcgCtx *ctx, int needed) {
    if (needed <= ctx->cap_btrail) return;
    int nc = ctx->cap_btrail ? ctx->cap_btrail * 2 : 32;
    while (nc < needed) nc *= 2;
    ctx->btrail = realloc(ctx->btrail, (size_t)nc * sizeof(BoundUndo));
    ctx->cap_btrail = nc;
}

static void grow_trail(LcgCtx *ctx, int needed) {
    if (needed <= ctx->cap_trail) return;
    int nc = ctx->cap_trail ? ctx->cap_trail * 2 : 32;
    while (nc < needed) nc *= 2;
    ctx->pre_btrail = realloc(ctx->pre_btrail, (size_t)nc * sizeof(int));
    ctx->trail_lits = realloc(ctx->trail_lits, (size_t)nc * sizeof(int));
    ctx->cap_trail = nc;
}

static void grow_propq(LcgCtx *ctx) {
    if (ctx->n_propq < ctx->cap_propq) return;
    int oc = ctx->cap_propq;
    int nc = oc ? oc * 2 : 8;
    ctx->propq = realloc(ctx->propq, (size_t)nc * sizeof(*ctx->propq));
    for (int i = oc; i < nc; i++) {
        ctx->propq[i].lit = 0;
        ctx->propq[i].reason = NULL;
        ctx->propq[i].n_reason = 0;
        ctx->propq[i].reason_cap = 0;
    }
    ctx->cap_propq = nc;
}

/* ====== allocating SAT atoms on demand ======
 *
 * Atom ids are just an integer counter; they don't refer to anything in
 * the SAT engine until lcg_solve creates it. The SAT engine will be
 * sized to exactly ctx->n_sat_atoms at solve time, so unused atoms do
 * not appear on its VSIDS heap and do not inflate the decision count.
 */

static int alloc_sat_atom(LcgCtx *ctx) {
    int atom = ++ctx->n_sat_atoms;
    grow_atom_back(ctx, atom);
    return atom;
}

/* Buffer a clause for later replay into the SAT engine. */
static void buffer_clause(LcgCtx *ctx, const int *lits, int n) {
    if (ctx->cl_buf_n == ctx->cl_buf_cap) {
        int nc = ctx->cl_buf_cap ? ctx->cl_buf_cap * 2 : 32;
        ctx->cl_buf = realloc(ctx->cl_buf, (size_t)nc * sizeof(*ctx->cl_buf));
        ctx->cl_buf_cap = nc;
    }
    struct ClauseBuf *b = &ctx->cl_buf[ctx->cl_buf_n++];
    b->n_lits = n;
    if (n > 0) {
        b->lits   = malloc((size_t)n * sizeof(int));
        memcpy(b->lits, lits, (size_t)n * sizeof(int));
    } else {
        b->lits = NULL;     /* empty clause; will produce immediate UNSAT */
    }
}

/* ====== lifecycle ====== */
static int   theory_check    (void *data, int **conflict_out, int *n_conflict_out);
static int   theory_propagate(void *data, int *prop_out,
                              int **reason_out, int *n_reason_out);
static void  theory_assign   (void *data, int lit);
static void  theory_unassign (void *data, int lit);

LcgCtx *lcg_new(void) {
    LcgCtx *ctx = calloc(1, sizeof(LcgCtx));
    /* No SAT engine yet -- created on demand in lcg_solve so we can
     * size it exactly to n_sat_atoms. */
    ctx->hooks.assign    = theory_assign;
    ctx->hooks.unassign  = theory_unassign;
    ctx->hooks.check     = theory_check;
    ctx->hooks.propagate = theory_propagate;
    ctx->hooks.data      = ctx;
    return ctx;
}

void lcg_free(LcgCtx *ctx) {
    if (!ctx) return;
    if (ctx->sat) sat_free(ctx->sat);
    for (int i = 0; i < ctx->n_vars; i++) {
        free(ctx->vars[i].name);
        free(ctx->vars[i].atoms);
    }
    free(ctx->vars);
    free(ctx->atom_back);
    for (int i = 0; i < ctx->n_props; i++) {
        if (ctx->props[i].destroy) ctx->props[i].destroy(ctx->props[i].data);
    }
    free(ctx->props);
    free(ctx->expl_neg);
    free(ctx->btrail);
    free(ctx->pre_btrail);
    free(ctx->trail_lits);
    for (int i = 0; i < ctx->cap_propq; i++) free(ctx->propq[i].reason);
    free(ctx->propq);
    free(ctx->conf_buf);
    /* Any clauses not yet replayed (only happens if user freed without
     * solving). */
    for (int i = 0; i < ctx->cl_buf_n; i++) free(ctx->cl_buf[i].lits);
    free(ctx->cl_buf);
    free(ctx);
}

/* ====== creating an integer variable ====== */
IntVar lcg_new_int_var(LcgCtx *ctx, int lo, int hi, const char *name) {
    assert(hi >= lo);
    grow_vars(ctx, ctx->n_vars + 1);
    int v = ctx->n_vars++;
    VarRec *r = &ctx->vars[v];
    r->name    = name ? strdup(name) : NULL;
    r->dom_lo  = lo;
    r->dom_hi  = hi;
    r->cur_lo  = lo;
    r->cur_hi  = hi;
    r->n_atoms = hi - lo;
    r->atoms   = (r->n_atoms > 0) ? calloc((size_t)r->n_atoms, sizeof(int))
                                  : NULL;
    /* Order-encoding atoms and monotonicity clauses. */
    for (int k = 0; k < r->n_atoms; k++) {
        int a = alloc_sat_atom(ctx);
        r->atoms[k] = a;
        ctx->atom_back[a].var = v;
        ctx->atom_back[a].threshold_idx = k;
    }
    /* Monotonicity: a[k] -> a[k+1], i.e. -a[k] v a[k+1]. */
    for (int k = 0; k + 1 < r->n_atoms; k++) {
        int cl[2] = { -r->atoms[k], r->atoms[k+1] };
        buffer_clause(ctx, cl, 2);
    }
    return v;
}

/* ====== literal accessors ====== */
int lcg_lit_le(LcgCtx *ctx, IntVar x, int k) {
    VarRec *r = &ctx->vars[x];
    if (k >= r->dom_hi)        return LCG_LIT_TRUE;
    if (k <  r->dom_lo)        return LCG_LIT_FALSE;
    int idx = k - r->dom_lo;
    return r->atoms[idx];
}
int lcg_lit_ge(LcgCtx *ctx, IntVar x, int k) {
    /* x >= k  iff  not (x <= k-1) */
    VarRec *r = &ctx->vars[x];
    if (k <= r->dom_lo)        return LCG_LIT_TRUE;
    if (k >  r->dom_hi)        return LCG_LIT_FALSE;
    int idx = (k - 1) - r->dom_lo;
    return -r->atoms[idx];
}

/* ====== posting constraints ====== */
void lcg_post_clause(LcgCtx *ctx, const int *lits, int n) {
    /* Filter trivial-true / trivial-false sentinels before passing on. */
    int *buf = malloc((size_t)n * sizeof(int));
    int  nb = 0;
    int  trivially_true = 0;
    for (int i = 0; i < n; i++) {
        int L = lits[i];
        if (L ==  LCG_LIT_TRUE)  { trivially_true = 1; break; }
        if (L == -LCG_LIT_TRUE)  continue;            /* false: drop */
        if (L ==  LCG_LIT_FALSE) continue;
        if (L == -LCG_LIT_FALSE) { trivially_true = 1; break; }
        buf[nb++] = L;
    }
    if (!trivially_true) {
        if (nb == 0) {
            /* empty clause: trivially UNSAT */
            buffer_clause(ctx, NULL, 0);
        } else {
            buffer_clause(ctx, buf, nb);
        }
    }
    free(buf);
}

/* ====== bound trail operations ====== */
static void push_bound_undo(LcgCtx *ctx, int var, int side, int old_value) {
    grow_btrail(ctx, ctx->n_btrail + 1);
    BoundUndo *u = &ctx->btrail[ctx->n_btrail++];
    u->var = var; u->side = side; u->old_value = old_value;
}

static void undo_one_bound(LcgCtx *ctx) {
    assert(ctx->n_btrail > 0);
    BoundUndo u = ctx->btrail[--ctx->n_btrail];
    VarRec *r = &ctx->vars[u.var];
    if (u.side == 0) r->cur_lo = u.old_value;
    else             r->cur_hi = u.old_value;
}

/* ====== theory_assign / theory_unassign ====== */
static void theory_assign(void *data, int lit) {
    LcgCtx *ctx = (LcgCtx *)data;
    grow_trail(ctx, ctx->n_trail + 1);
    ctx->pre_btrail[ctx->n_trail] = ctx->n_btrail;
    ctx->trail_lits[ctx->n_trail++] = lit;

    int atom = lit > 0 ? lit : -lit;
    if (atom > ctx->n_sat_atoms) return;
    AtomBack *ab = &ctx->atom_back[atom];
    if (ab->var < 0) return;     /* not an order-encoding atom */

    VarRec *r = &ctx->vars[ab->var];
    int threshold = r->dom_lo + ab->threshold_idx;
    if (lit > 0) {
        /* x <= threshold  is true: tighten hi if useful. */
        if (r->cur_hi > threshold) {
            push_bound_undo(ctx, ab->var, 1, r->cur_hi);
            r->cur_hi = threshold;
        }
    } else {
        /* x <= threshold  is false, i.e. x > threshold: tighten lo. */
        if (r->cur_lo < threshold + 1) {
            push_bound_undo(ctx, ab->var, 0, r->cur_lo);
            r->cur_lo = threshold + 1;
        }
    }
}

static void theory_unassign(void *data, int lit) {
    LcgCtx *ctx = (LcgCtx *)data;
    assert(ctx->n_trail > 0);
    assert(ctx->trail_lits[ctx->n_trail - 1] == lit);
    ctx->n_trail--;
    int target = ctx->pre_btrail[ctx->n_trail];
    while (ctx->n_btrail > target) undo_one_bound(ctx);

    /* Drop any pending propagation entries that the SAT engine has not
     * yet consumed. They were derived under bounds that this unassign
     * is rolling back; their reason clauses may no longer reflect the
     * trail and would feed malformed clauses into the SAT engine if
     * replayed by theory_propagate. Consumed entries (already at
     * indices < propq_head) are fine: SAT has installed their reason
     * clauses as learnts and will manage them via watched literals. */
    if (ctx->n_propq > ctx->propq_head) {
        ctx->n_propq = ctx->propq_head;
    }
}

/* ====== propagation queue helpers ====== */
static void enqueue_prop(LcgCtx *ctx, int lit,
                         const int *reason, int n_reason)
{
    grow_propq(ctx);
    struct PropEntry *p = &ctx->propq[ctx->n_propq++];
    if (p->reason_cap < n_reason) {
        p->reason_cap = n_reason;
        p->reason = realloc(p->reason, (size_t)n_reason * sizeof(int));
    }
    memcpy(p->reason, reason, (size_t)n_reason * sizeof(int));
    p->n_reason = n_reason;
    p->lit      = lit;
    ctx->n_explanations++;
}

/* Reason format mirrors the EUF/LRA convention:
 *     reason[0]      = the propagated literal
 *     reason[1..n-1] = negations of the trail literals that imply it.
 */

/* ====== Phase 2: explanation builder ====== */

/* Grow the explanation buffer to fit at least `needed` literals. */
static void expl_grow(LcgCtx *ctx, int needed) {
    if (needed <= ctx->expl_cap) return;
    int nc = ctx->expl_cap ? ctx->expl_cap * 2 : 16;
    while (nc < needed) nc *= 2;
    ctx->expl_neg = realloc(ctx->expl_neg, (size_t)nc * sizeof(int));
    ctx->expl_cap = nc;
}

void lcg_expl_clear(LcgCtx *ctx) {
    ctx->expl_n = 0;
}

/* Support: "v >= lo". The asserted literal on the trail is -lit_le(v, lo-1).
 * In a conflict/reason clause we want its negation: +lit_le(v, lo-1).
 *
 * Edge cases (sentinels):
 *   lo <= dom_lo(v): support is trivially true; nothing to record.
 *   lo >  dom_hi(v): support can never hold; the caller is asking for
 *     an impossible support, which usually indicates a propagator bug,
 *     but we just drop it.
 */
void lcg_expl_lo(LcgCtx *ctx, IntVar v, int lo) {
    int L = lcg_lit_le(ctx, v, lo - 1);
    if (L == LCG_LIT_TRUE || L == LCG_LIT_FALSE) return;
    expl_grow(ctx, ctx->expl_n + 1);
    ctx->expl_neg[ctx->expl_n++] = L;
}

/* Support: "v <= hi". Asserted: +lit_le(v, hi). Negation: -lit_le(v, hi). */
void lcg_expl_hi(LcgCtx *ctx, IntVar v, int hi) {
    int L = lcg_lit_le(ctx, v, hi);
    if (L == LCG_LIT_TRUE || L == LCG_LIT_FALSE) return;
    expl_grow(ctx, ctx->expl_n + 1);
    ctx->expl_neg[ctx->expl_n++] = -L;
}

/* Emit "v <= k" with the current explanation as justification.
 * Returns LCG_PROP_OK if the literal is already known, _TIGHTENED if
 * we enqueued a propagation, _CONFLICT if asserting it conflicts with
 * existing bounds. */
int lcg_emit_le(LcgCtx *ctx, IntVar v, int k) {
    int target = lcg_lit_le(ctx, v, k);
    if (target == LCG_LIT_TRUE)  return LCG_PROP_OK;  /* already implied */
    if (target == LCG_LIT_FALSE) return LCG_PROP_CONFLICT;
    if (ctx->sat && sat_value(ctx->sat, target) == 1) return LCG_PROP_OK;

    int n = ctx->expl_n;
    int *reason = malloc((size_t)(n + 1) * sizeof(int));
    reason[0] = target;
    memcpy(reason + 1, ctx->expl_neg, (size_t)n * sizeof(int));
    enqueue_prop(ctx, target, reason, n + 1);
    free(reason);
    return LCG_PROP_TIGHTENED;
}

/* Emit "v >= k". Implemented as the negation of "v <= k-1". */
int lcg_emit_ge(LcgCtx *ctx, IntVar v, int k) {
    int target_le = lcg_lit_le(ctx, v, k - 1);
    if (target_le == LCG_LIT_FALSE) return LCG_PROP_OK;       /* v > k-1 always */
    if (target_le == LCG_LIT_TRUE)  return LCG_PROP_CONFLICT; /* v <= k-1 forced */
    int target = -target_le;
    if (ctx->sat && sat_value(ctx->sat, target_le) == 0) return LCG_PROP_OK;

    int n = ctx->expl_n;
    int *reason = malloc((size_t)(n + 1) * sizeof(int));
    reason[0] = target;
    memcpy(reason + 1, ctx->expl_neg, (size_t)n * sizeof(int));
    enqueue_prop(ctx, target, reason, n + 1);
    free(reason);
    return LCG_PROP_TIGHTENED;
}

int lcg_emit_conflict(LcgCtx *ctx) {
    (void)ctx;
    return LCG_PROP_CONFLICT;
}

/* Bounds accessors for use inside propagator code. */
int lcg_lb(const LcgCtx *ctx, IntVar v) { return ctx->vars[v].cur_lo; }
int lcg_ub(const LcgCtx *ctx, IntVar v) { return ctx->vars[v].cur_hi; }

/* ====== propagator registration ====== */
void lcg_register_propagator(LcgCtx *ctx, LcgPropagator p) {
    if (ctx->n_props == ctx->cap_props) {
        int nc = ctx->cap_props ? ctx->cap_props * 2 : 8;
        ctx->props = realloc(ctx->props, (size_t)nc * sizeof(LcgPropagator));
        ctx->cap_props = nc;
    }
    ctx->props[ctx->n_props++] = p;
}

/* ====== linear_le propagator ======
 *
 * Constraint:  sum a_i x_i  <=  K.
 *
 * For each variable x_j with coefficient a_j, isolate it:
 *
 *     a_j x_j  <=  K - sum_{i != j} a_i x_i
 *
 * Now bound the right-hand side. For each i != j with a_i > 0 use lo[i]
 * (since x_i is at least lo[i], the term a_i x_i is at least a_i lo[i],
 * so dropping it shrinks the RHS least when x_i is at its lower bound).
 * For each i != j with a_i < 0 use hi[i] (symmetric argument). Then:
 *
 *   If a_j > 0:    x_j <= (K - residual) / a_j         (upper bound tighten)
 *   If a_j < 0:    x_j >= (K - residual) / a_j         (lower bound tighten)
 *
 * The explanation for each tightening is the set of bound atoms used to
 * derive `residual`: for positive-coefficient i, the literal x_i >= lo[i];
 * for negative-coefficient i, the literal x_i <= hi[i].
 *
 * Phase 1: integer coefficients, integer rhs. No rational arithmetic
 * needed; we use floor/ceil division on the rounded bound. Empty domain
 * (lo > hi after tightening) is a conflict; we build a clause from the
 * supporting literals.
 */
/* Phase 2 rewrite of the propagator: same math, declarative explanations.
 *
 * For each linear_le constraint  sum a_i x_i <= K  we run two passes:
 *   1. If sum_min > K, the bounds are jointly infeasible. Build the
 *      explanation from the supporting bounds and emit a conflict.
 *   2. Otherwise, for each variable x_j try to tighten its upper bound
 *      (if a_j > 0) or lower bound (if a_j < 0). The explanation is
 *      the bounds of the OTHER variables that fixed the slack.
 *
 * The explanation logic --- which other-variable bounds to cite --- is
 * the same as in Phase 1 but is now expressed declaratively through
 * lcg_expl_lo / lcg_expl_hi. The arithmetic for the new bound is also
 * unchanged (floor/ceil division on the rounded slack).
 */

typedef struct {
    int   n_terms;
    int  *coeff;   /* heap-owned */
    int  *vars;    /* heap-owned */
    int   rhs;
} LinLePropData;

static void linear_le_destroy(void *data) {
    LinLePropData *d = (LinLePropData *)data;
    free(d->coeff);
    free(d->vars);
    free(d);
}

/* Add the supporting bound for OTHER terms (i != j) to ctx's explanation
 * buffer. For positive coefficients, the support is the variable's
 * current lower bound; for negative, its current upper bound. */
static void linear_le_add_supports(LcgCtx *ctx, LinLePropData *d, int skip) {
    for (int i = 0; i < d->n_terms; i++) {
        if (i == skip) continue;
        if (d->coeff[i] >= 0) lcg_expl_lo(ctx, d->vars[i], lcg_lb(ctx, d->vars[i]));
        else                  lcg_expl_hi(ctx, d->vars[i], lcg_ub(ctx, d->vars[i]));
    }
}

static int linear_le_propagate(LcgCtx *ctx, void *data) {
    LinLePropData *d = (LinLePropData *)data;

    long long sum_min = 0;
    for (int i = 0; i < d->n_terms; i++) {
        long long a  = d->coeff[i];
        long long lo = lcg_lb(ctx, d->vars[i]);
        long long hi = lcg_ub(ctx, d->vars[i]);
        sum_min += (a >= 0) ? a * lo : a * hi;
    }

    /* Pass 1: infeasibility. */
    if (sum_min > d->rhs) {
        lcg_expl_clear(ctx);
        linear_le_add_supports(ctx, d, /*skip=*/-1);   /* every term contributes */
        return lcg_emit_conflict(ctx);
    }

    /* Pass 2: per-variable tightening. */
    int tightened = 0;
    for (int j = 0; j < d->n_terms; j++) {
        long long a = d->coeff[j];
        if (a == 0) continue;
        long long residual_min = (a >= 0)
            ? sum_min - a * (long long)lcg_lb(ctx, d->vars[j])
            : sum_min - a * (long long)lcg_ub(ctx, d->vars[j]);
        long long room = (long long)d->rhs - residual_min;   /* a * x_j <= room */

        if (a > 0) {
            /* x_j <= floor(room / a) */
            long long bound = (room >= 0) ? room / a
                                          : -(((-room) + a - 1) / a);
            if (bound < lcg_ub(ctx, d->vars[j])) {
                lcg_expl_clear(ctx);
                linear_le_add_supports(ctx, d, j);
                int r = lcg_emit_le(ctx, d->vars[j], (int)bound);
                if (r == LCG_PROP_CONFLICT) return LCG_PROP_CONFLICT;
                if (r == LCG_PROP_TIGHTENED) tightened = 1;
            }
        } else {
            /* x_j >= ceil(room / a); a is negative so division direction flips */
            long long na = -a;
            long long bound;
            if (room >= 0) bound = -(room / na);
            else           bound = ((-room) + na - 1) / na;
            if (bound > lcg_lb(ctx, d->vars[j])) {
                lcg_expl_clear(ctx);
                linear_le_add_supports(ctx, d, j);
                int r = lcg_emit_ge(ctx, d->vars[j], (int)bound);
                if (r == LCG_PROP_CONFLICT) return LCG_PROP_CONFLICT;
                if (r == LCG_PROP_TIGHTENED) tightened = 1;
            }
        }
    }
    return tightened ? LCG_PROP_TIGHTENED : LCG_PROP_OK;
}

void lcg_post_linear_le(LcgCtx *ctx,
                        const int *coeff, const IntVar *vars, int n_terms,
                        int k)
{
    LinLePropData *d = malloc(sizeof(LinLePropData));
    d->n_terms = n_terms;
    d->coeff   = malloc((size_t)n_terms * sizeof(int));
    d->vars    = malloc((size_t)n_terms * sizeof(int));
    memcpy(d->coeff, coeff, (size_t)n_terms * sizeof(int));
    memcpy(d->vars,  vars,  (size_t)n_terms * sizeof(int));
    d->rhs     = k;
    LcgPropagator p = { .data = d,
                        .propagate = linear_le_propagate,
                        .destroy   = linear_le_destroy };
    lcg_register_propagator(ctx, p);
}

/* ====== theory_check: run all registered propagators once ======
 *
 * Each propagator may emit propagations (via lcg_emit_le / lcg_emit_ge)
 * or a conflict. We do NOT loop until fixpoint inside theory_check:
 * a propagator's emit_le enqueues the literal but does not update
 * cur_lo / cur_hi (that only happens when SAT calls theory_assign
 * with the literal). So a second pass would re-derive the same bound
 * and re-enqueue it forever. SAT drains the queue between calls and
 * re-invokes theory_check; that is the real fixpoint loop.
 */
static int theory_check(void *data, int **conflict_out, int *n_conflict_out) {
    LcgCtx *ctx = (LcgCtx *)data;
    for (int i = 0; i < ctx->n_props; i++) {
        int r = ctx->props[i].propagate(ctx, ctx->props[i].data);
        if (r == LCG_PROP_CONFLICT) {
            /* The propagator already populated ctx->expl_neg via
             * lcg_expl_lo / _hi; that array IS the conflict clause. */
            *conflict_out   = ctx->expl_neg;
            *n_conflict_out = ctx->expl_n;
            return 1;
        }
    }
    return 0;
}

static int theory_propagate(void *data, int *prop_out,
                            int **reason_out, int *n_reason_out)
{
    LcgCtx *ctx = (LcgCtx *)data;
    while (ctx->propq_head < ctx->n_propq) {
        struct PropEntry *p = &ctx->propq[ctx->propq_head++];
        if (p->lit == 0) continue;
        /* Skip already-assigned lits (their effect is harmless but redundant). */
        int atom = p->lit > 0 ? p->lit : -p->lit;
        if (sat_value(ctx->sat, atom) != -1) continue;
        *prop_out      = p->lit;
        *reason_out    = p->reason;
        *n_reason_out  = p->n_reason;
        return 1;
    }
    return 0;
}

/* ====== solve ====== */
int lcg_solve(LcgCtx *ctx) {
    if (ctx->solved) return LCG_UNKNOWN;
    ctx->solved = 1;

    /* Edge case: no SAT atoms were ever allocated (every variable is
     * pinned to a single value, and no user clauses were posted). The
     * SAT engine has nothing to decide, but propagators may still detect
     * infeasibility on the initial bounds --- run them once. */
    if (ctx->n_sat_atoms < 1) {
        for (int i = 0; i < ctx->n_props; i++) {
            if (ctx->props[i].propagate(ctx, ctx->props[i].data)
                    == LCG_PROP_CONFLICT) {
                return LCG_UNSAT;
            }
        }
        return LCG_SAT;
    }

    ctx->sat = sat_new(ctx->n_sat_atoms);
    sat_set_verbose(ctx->sat, ctx->verbose);
    sat_set_theory (ctx->sat, &ctx->hooks);

    /* Replay buffered clauses. The empty clause (n_lits==0) makes the
     * formula trivially UNSAT; sat_add_clause handles this when passed
     * a 0-length array. We pass a 1-element dummy for that case (sat's
     * API requires a non-NULL pointer). */
    for (int i = 0; i < ctx->cl_buf_n; i++) {
        struct ClauseBuf *b = &ctx->cl_buf[i];
        if (b->n_lits == 0) {
            int dummy = 0;
            sat_add_clause(ctx->sat, &dummy, 0);
        } else {
            sat_add_clause(ctx->sat, b->lits, b->n_lits);
        }
        free(b->lits);
        b->lits = NULL;
    }
    free(ctx->cl_buf);
    ctx->cl_buf = NULL;
    ctx->cl_buf_n = ctx->cl_buf_cap = 0;

    int r = sat_solve(ctx->sat);
    if (r == SAT_SATISFIABLE)   return LCG_SAT;
    if (r == SAT_UNSATISFIABLE) return LCG_UNSAT;
    return LCG_UNKNOWN;
}

/* ====== model query ====== */
int lcg_value(const LcgCtx *ctx, IntVar x) {
    const VarRec *r = &ctx->vars[x];
    /* The current bounds may have collapsed to a singleton; if so,
     * return it. Otherwise look at the order-encoding atoms to find
     * the smallest k such that x <= k is true. */
    if (r->cur_lo == r->cur_hi) return r->cur_lo;
    if (!ctx->sat) return r->dom_lo;   /* solve hasn't run */
    for (int k = 0; k < r->n_atoms; k++) {
        int a = r->atoms[k];
        int v = sat_value(ctx->sat, a);
        if (v == 1) return r->dom_lo + k;
    }
    /* Last threshold is "x <= hi - 1"; if all of them are false, x == hi. */
    return r->dom_hi;
}

/* ====== stats ====== */
long lcg_decisions    (const LcgCtx *ctx) {
    return ctx->sat ? sat_decisions(ctx->sat) : 0;
}
long lcg_conflicts    (const LcgCtx *ctx) {
    return ctx->sat ? sat_conflicts(ctx->sat) : 0;
}
long lcg_propagations (const LcgCtx *ctx) {
    return ctx->sat ? sat_propagations(ctx->sat) : 0;
}
long lcg_explanations (const LcgCtx *ctx) { return ctx->n_explanations; }

void lcg_set_verbose(LcgCtx *ctx, int level) {
    ctx->verbose = level;
    if (ctx->sat) sat_set_verbose(ctx->sat, level);
}
